# Untitled

unknown
plain_text
a year ago
6.4 kB
2
Indexable
Never
```// matrix_operations in 3D
#include <cmath>
#include "vector.hpp"
#include "matrix.hpp"

void product_scalar_by_matrix(float scalar, Matrix *M, int n, int m, Matrix *R)
{
for (int r = 0; r < n; r++)
for (int c = 0; c < m; c++)
R->set(scalar * M->get(r, c), r, c);
}

void product_matrix_by_vector(Matrix *M, Vector *V, int n, int m, Vector *R)
{
for (int r = 0; r < n; r++)
{
float acc_n = 0;
for (int c = 0; c < n; c++)
acc_n += M->get(r, c) * V->get(c);
R->set(acc_n, r);
}
}

// Funcion que se usara en MEF_Process
void product_matrix_by_matrix(Matrix *A, Matrix *B, Matrix *R)
{
int n = A->getNrows(), m = A->getNcols(), p = B->getNrows(), q = B->getNcols();
if (m == p)
{
R->setSize(n, q);
R->init();

for (int r = 0; r < n; r++)
for (int c = 0; c < q; c++)
for (int i = 0; i < m; i++)
R->add(A->get(r, i) * B->get(i, c), r, c);
}
else
{
cout << "*****Dimension incompatibility when multiplying matrices.\n\n";
exit(EXIT_FAILURE);
}
}

float determinant(Matrix *M);

float determinant_auxiliar(Matrix *M)
{
int n = M->getNcols();
float acc_n = 0;

for (int c = 0; c < n; c++)
{
Matrix clon(n, n);
M->clone(&clon);
clon.removeRow(0);
clon.removeColumn(c);
acc_n += pow(-1, c) * M->get(0, c) * determinant(&clon);
}
return acc_n;
}

float determinant(Matrix *M)
{
float ans;
switch (M->getNcols())
{
case 1:
ans = M->get(0, 0);
break;
case 2:
ans = M->get(0, 0) * M->get(1, 1) - M->get(0, 1) * M->get(1, 0);
break;
case 3:
ans = M->get(0, 0) * M->get(1, 1) * M->get(2, 2) - M->get(0, 0) * M->get(1, 2) * M->get(2, 1) - M->get(0, 1) * M->get(1, 0) * M->get(2, 2) + M->get(0, 1) * M->get(1, 2) * M->get(2, 0) + M->get(0, 2) * M->get(1, 0) * M->get(2, 1) - M->get(0, 2) * M->get(1, 1) * M->get(2, 0);
break;

default:
ans = determinant_auxiliar(M);
}
return ans;
}

/*float getMinor(Matrix *M, int n, int r, int c)
{
Matrix clon(n, n);
M->clone(&clon);

// M.show();
// clon.show();
clon.removeRow(r);
clon.removeColumn(c);

// clon.show();
return determinant(&clon);
}*/

/*void conjugate_matrix(Matrix *M, int n, Matrix *C)
{
for (int r = 0; r < n; r++)
for (int c = 0; c < n; c++)
C->set(pow(-1, r + c) * getMinor(M, n, r, c), r, c);
}*/

// Funcion transpose
void transpose_matrix(Matrix *M, int n, int m, Matrix *T)
{
for (int r = 0; r < n; r++)
for (int c = 0; c < m; c++)
T->set(M->get(r, c), c, r);
}

void calculateZ_matrix(Matrix *A, int n, Matrix *Z)
{
// utilizando archivo de Cholesky del repositorio de github
float acum;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
if (i == j)
{
acum = 0;
for (int k = 0; k < j; k++)
{
// cout << L->get(j,k) << " ";
acum += pow(X->get(k, j), 2);
}
cout << A->get(j, j) << " " << acum << "\n";
if (A->get(j, j) - acum >= 0)
{
Z->set(sqrt(A->get(j, j) - acum), j, j);
}
else
{
Z->set(0.000001, j, j);
}
}
else
{
if (i > j)
{
if (X->get(j, j) != 0)
{
acum = 0;
for (int k = 0; k < j; k++)
{
acum += Z->get(i, k) * Z->get(j, k);
}
// cout << L->get(j,j) << "\n";
Z->set((1 / Z->get(j, j)) * (A->get(i, j) - acum), i, j);
}
else
{
Z->set(0.000001, i, j);
}
}
else
{
Z->set(0, i, j);
}
}
}
}
}

//mirar sino hay error en sintaxis
void calculateY_matrix(Matrix *Z, int n, Matrix *Y)
{
// utilizando archivo de Cholesky del repositorio de github
float acum;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
if (j == i)
{
if (Z->get(i, j) != 0)
{
Y->set(1 /Z->get(i, i), i, i);
}
else
{
Z->set(0.000001, i, j);
}
}
else
{
if (i > j)
{
if (Z->get(i, i) != 0)
{
acum = 0;
for (int k = j; k < i; k++)
{
acum += Z->get(i, k) * Y->get(k, j);
}
Y->set(-1(1 / Z->get(i, i)) * acum, i, j);
}
else{
Z->set(0.000001, i, j);
}
}
}
else{
Y->set(0, i, j);
}
}
}
}
void calculateX_matrix(Matrix* Y, Matrix* Z, int n, Matrix*X){
float acum;
for(int i = n-1 ; i>=0; i--)
for(int j = 0; j<n; j++){
if(Z->get(i,i)!=0){
acum = 0;
for(int k = i+1; k<n; k++){
acum += Z->get(k,i)*X->get(k,j);
}
X->set((1/Z->get(i,i))*(Y->get(i,j)-acum),i,j);
}
else{
X->set(0.000001,i,j);
}
}

}

void calculate_inverse(Matrix *M, int n, Matrix *X)
{
//Aplicando las formulas de las matriz de X,Y y Z
Matrix Z(n, n), Y(n,n);
Z.init();
Y.init();

calculateX_matrix(A,n,&L);
calculateY_matrix(&Z,n,&Y);
calculateZ_matrix(&Y,&Z,n,X);

}```