Page MenuHomec4science
Files F63138011

Functions.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, May 18, 00:39

Functions.cpp
View Options

#include <iostream>
#include <cassert>
#include <cmath>
#include "Functions.hpp"
#include "Conjugate_gradient.hpp"
double** AllocateMemory(int rows, int cols){
/*Allocates memory for matrix A*/
double** A;
A=new double* [rows];
for (int i=0;i<rows;i++){
A[i]=new double[cols];
}
return A;
}
double* Substract(double* u_1, double* u_2, int size){
double* temp=AllocateMemory(size);
for (int i=0;i<size;i++) {
temp[i]=u_1[i]-u_2[i];
}
return temp;
}
double* Add(double* u_1, double* u_2,int size){
double* temp=AllocateMemory(size);
for (int i=0;i<size;i++) {
temp[i]=u_1[i]+u_2[i];
}
return temp;
}
double** Initialize(int size1,int size2){
double** A=AllocateMemory(size1,size2);
for (int i=0;i<size1;i++) {
for (int j = 0; j < size2; j++) {
A[i][j] = 0;
}
}
for (int i=0;i<size1;i++){
A[i][i]=1;
}
return A;
}
double* Initialize(int size){
double* x_0=AllocateMemory(size);
for (int i=0;i<size;i++){
x_0[i]=0;
}
return x_0;
}
double** Copy(double** A,int size){
double** B=AllocateMemory(size,size);
for (int i=0;i<size;i++){
for (int j=0;j<size;j++){
B[i][j]=A[i][j];
}
}
return B;
}
double* Copy(double* b,int size){
double* b_star=AllocateMemory(size);
for (int i=0;i<size;i++){
b_star[i]=b[i];
}
return b_star;
}
void SubstractRows (double** A,double scalar,int row_1, int row_2, int size){
for (int j=0;j<size;j++){
A[row_1][j]-=A[row_2][j]*scalar;
}
}
void SubstractRows (double* u,double scalar,int row_1, int row_2){
u[row_1]-=u[row_2]*scalar;
}
double Normalize (double** A,int row, int size){
double scalar=A[row][row];
for (int j=0;j<size;j++){
A[row][j]/=scalar;
}
return scalar;
}
LU_matrices Pivot(double** A,double* b,int size) {
double** L=Initialize(size,size);
double** U=Copy(A,size);
double* b_star=Copy(b,size);
LU_matrices output;
for (int i=0;i<size-1;i++){
for (int j=i+1;j<size;j++){
L[j][i]=U[j][i]/U[i][i];
SubstractRows(U,U[j][i]/U[i][i],j,i,size);
SubstractRows(b_star,L[j][i],j,i);
}
}
output.U_matrix=U;
output.L_matrix=L;
output.b_star_vector=b_star;
return output;
}
double* AllocateMemory(int rows){
/*Allocates memory for vector u*/
double* u=new double [rows];
return u;
}
void FreeMemory(double** A,int rows){
/*Frees memory of matrix A*/
for (int i=0;i<rows;i++){
delete[] A[i];
}
delete[] A;
}
void FreeMemory(double* u){
/*Frees memory of vector u*/
delete[] u;
}
int FindMax (double** A,int col,int size){
int index=0;
for (int i=0;i<size;i++){
if (A[i][col]>A[index][col]){
index=i;
}
}
return index;
}
void SwapRows (double** A,int last_index, int new_index, int size){
double* temp=AllocateMemory(size);
for (int i=0;i<size;i++){
temp[i]=A[last_index][i];
A[last_index][i]=A[new_index][i];
A[new_index][i]=temp[i];
}
FreeMemory(temp);
}
void TerminalPrint(double** A,int rows, int cols){
/*Prints to terminal values of a matrix A*/
for (int i=0;i<rows;i++){
for (int j=0;j<cols;j++){
std::cout<<A[i][j]<<" ";
}
std::cout<<std::endl;
}
std::cout<<std::endl;
}
void TerminalPrint(double* u,int size){
/*Prints to terminal values of a vector u*/
for (int i=0;i<size;i++) {
std::cout << u[i]<<std::endl;
}
std::cout<<std::endl;
}
double** Multiply(double** A,double** B,int rowsA, int colsA, int rowsB, int colsB){
/*Multiplies matrix A[rowsA][colsA] by matrix B[rowsB][colsB]*/
assert(rowsB==colsA);
double** C=AllocateMemory(rowsA,colsB);
for (int i=0;i<rowsA;i++){
for (int j=0;j<colsB;j++){
for (int k=0;k<rowsB;k++){
C[i][j]+=A[i][k]*B[k][j];
}
}
}
return C;
}
double* Multiply(double** A,double* u,int rowsA, int colsA, int rowsU){
/*Multiplies matrix A[rowsA][colsA] by vector U[rowsU]*/
assert(rowsU==colsA);
double* v=AllocateMemory(rowsA);
for (int i=0;i<rowsA;i++){
for (int j=0;j<colsA;j++){
v[i]+=A[i][j]*u[j];
}
}
return v;
}
double* Multiply(double* u,double** A,int rowsA, int colsA, int rowsU){
/*Multiplies vector U[rowsU] by matrix A[rowsA][colsA]*/
assert(rowsU==rowsA);
double* v=AllocateMemory(colsA);
for (int i=0;i<colsA;i++){
for (int j=0;j<rowsA;j++){
v[i]+=A[j][i]*u[j];
}
}
return v;
}
double** Multiply(double** A,double scalar, int rowsA, int colsA){
/*Multiplies matrix A[rowsA][colsA] by scalar*/
double** C=AllocateMemory(rowsA,colsA);
for (int i=0;i<rowsA;i++){
for (int j=0;j<colsA;j++){
C[i][j]=A[i][j]*scalar;
}
}
return C;
}
double Multiply(double* u_1,double* u_2,int size){
double temp=0;
for (int i=0;i<size;i++){
temp+=u_1[i]*u_2[i];
}
return temp;
}
double* Multiply(double* u,double scalar, int rowsU){
/*Multiplies vector u[rowsU] by scalar*/
double* v=AllocateMemory(rowsU);
for (int i=0;i<rowsU;i++){
v[i]=u[i]*scalar;
}
return v;
}
void ReduceMatrix(double** A,int size,int col, double** &A_hat){
/*Resizing of matrix A for determinant computation*/
for (int row = 1; row < size; ++row) {
for (int k = 0; k < col; ++k) {
A_hat[row - 1][k] = A[row][k];
}
for (int k = col + 1; k < size; ++k) {
A_hat[row - 1][k - 1] = A[row][k];
}
}
}
double Determinant(double** A,int size){
/*Computation of determinant : direct formula for 2x2, recursive otherwise*/
if (size==2){
return A[0][0]*A[1][1]-A[1][0]*A[0][1];
} else {
double det=0;
for (int j = 0; j < size; ++j) {
double** A_hat=AllocateMemory(size-1,size-1);
ReduceMatrix(A,size,j,A_hat);
det +=A[0][j]*std::pow(-1.0,j)*Determinant(A_hat,size-1);
FreeMemory(A_hat,size-1);
}
return det;
}
}
double* Up(double** U, double* b_star, int size){
double* x=AllocateMemory(size);
for (int i=2;i>-1;i--){
x[i]+=b_star[i];
for (int j=i+1;j<size;j++){
x[i]-=U[i][j]*x[j];
}
x[i]/=U[i][i];
}
return x;
}
double* Down(double** L, double* b_star, int size){
double* x=AllocateMemory(size);
for (int i=0;i<size;i++){
x[i]+=b_star[i];
for (int j=i-1;j>=0;j--){
x[i]-=L[i][j]*x[j];
}
x[i]/=L[i][i];
}
return x;
}
double** Transpose(double** A,int size){
double** temp=Copy(A,size);
for (int i=0;i<size;i++){
for (int j=0;j<size;j++){
temp[i][j]=A[j][i];
}
}
return temp;
}
Cholesky_matrices Cholesky_L(double** A,double* b,int size){
double** L=AllocateMemory(size,size);
double* b_star=Copy(b,size);
Cholesky_matrices output;
for (int i=0;i<size;i++){
for (int k=0;k<i+1;k++){
double sum=0;
for (int j=0;j<k;j++){
sum+=(L[i][j]*L[k][j]);
}
if (i==k){
L[i][k]=sqrt(A[i][i]-sum);
} else {
L[i][k]=(1.0/L[k][k])*(A[i][k]-sum);
}
}
}
output.L_matrix=L;
output.L_T_matrix=Transpose(L,size);
output.b_star_vector=Down(L,b,size);
return output;
}

Event Timeline

c4science · Help