Skip to content

Commit b5db893

Browse files
fpalaciosfpalacios
committed
Better method to compute the inverse of a square matrix.
1 parent e9d3ca8 commit b5db893

File tree

4 files changed

+312
-21
lines changed

4 files changed

+312
-21
lines changed

‎Common/include/matrix_structure.hpp‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -239,6 +239,29 @@ class CSysMatrix {
239239
*/
240240
void DeleteValsRowi(unsigned long i);
241241

242+
/*!
243+
* \brief Recursive definition of determinate using expansion by minors. Written by Paul Bourke
244+
* \param[in] a - Matrix to compute the determinant.
245+
* \param[in] n - Size of the quare matrix.
246+
* \return Value of the determinant.
247+
*/
248+
su2double MatrixDeterminant(su2double **a, unsigned long n);
249+
250+
/*!
251+
* \brief Find the cofactor matrix of a square matrix. Written by Paul Bourke
252+
* \param[in] a - Matrix to compute the determinant.
253+
* \param[in] n - Size of the quare matrix.
254+
* \param[out] b - cofactor matrix
255+
*/
256+
void MatrixCoFactor(su2double **a, unsigned long n, su2double **b) ;
257+
258+
/*!
259+
* \brief Transpose of a square matrix, do it in place. Written by Paul Bourke
260+
* \param[in] a - Matrix to compute the determinant.
261+
* \param[in] n - Size of the quare matrix.
262+
*/
263+
void MatrixTranspose(su2double **a, unsigned long n) ;
264+
242265
/*!
243266
* \brief Performs the Gauss Elimination algorithm to solve the linear subsystem of the (i, i) subblock and rhs.
244267
* \param[in] block_i - Index of the (i, i) subblock in the matrix-by-blocks structure.

‎Common/src/grid_movement_structure.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1470,7 +1470,7 @@ void CVolumetricMovement::SetFEA_StiffMatrix2D(CGeometry *geometry, CConfig *con
14701470
switch (config->GetDeform_Stiffness_Type()) {
14711471
case INVERSE_VOLUME: E = 1.0 / ElemVolume; break;
14721472
case WALL_DISTANCE: E = 1.0 / ElemDistance; break;
1473-
case CONSTANT_STIFFNESS: E=1.0/EPS; break;
1473+
case CONSTANT_STIFFNESS: E = 1.0 / EPS; break;
14741474
}
14751475

14761476
Nu = config->GetDeform_Coeff();
@@ -1609,7 +1609,7 @@ void CVolumetricMovement::SetFEA_StiffMatrix3D(CGeometry *geometry, CConfig *con
16091609
switch (config->GetDeform_Stiffness_Type()) {
16101610
case INVERSE_VOLUME: E = 1.0 / ElemVolume; break;
16111611
case WALL_DISTANCE: E = 1.0 / ElemDistance; break;
1612-
case CONSTANT_STIFFNESS: E=1.0/EPS; break;
1612+
case CONSTANT_STIFFNESS: E = 1.0 / EPS; break;
16131613
}
16141614

16151615
Nu = config->GetDeform_Coeff();

‎Common/src/matrix_structure.cpp‎

Lines changed: 183 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -489,6 +489,102 @@ void CSysMatrix::DeleteValsRowi(unsigned long i) {
489489

490490
}
491491

492+
493+
su2double CSysMatrix::MatrixDeterminant(su2double **a, unsigned long n) {
494+
495+
unsigned long i, j, j1, j2;
496+
su2double det = 0;
497+
su2double **m = NULL;
498+
499+
if (n < 1) { }
500+
else if (n == 1) { det = a[0][0]; }
501+
else if (n == 2) { det = a[0][0] * a[1][1] - a[1][0] * a[0][1]; }
502+
else {
503+
det = 0.0;
504+
505+
for (j1=0;j1<n;j1++) {
506+
m = new su2double*[n-1];
507+
for (i=0;i<n-1;i++)
508+
m[i] = new su2double[n-1];
509+
510+
for (i=1;i<n;i++) {
511+
j2 = 0;
512+
for (j=0;j<n;j++) {
513+
if (j == j1)
514+
continue;
515+
m[i-1][j2] = a[i][j];
516+
j2++;
517+
}
518+
}
519+
520+
det += pow(-1.0,j1+2.0) * a[0][j1] * MatrixDeterminant(m,n-1);
521+
for (i=0;i<n-1;i++)
522+
delete [] m[i];
523+
delete [] m;
524+
}
525+
526+
}
527+
528+
return(det);
529+
530+
}
531+
532+
void CSysMatrix::MatrixCoFactor(su2double **a, unsigned long n, su2double **b) {
533+
534+
unsigned long i,j,ii,jj,i1,j1;
535+
su2double det;
536+
su2double **c;
537+
538+
c = new su2double*[n-1];
539+
for (i=0;i<n-1;i++)
540+
c[i] = new su2double[n-1];
541+
542+
for (j=0;j<n;j++) {
543+
for (i=0;i<n;i++) {
544+
545+
/*--- Form the adjoint a_ij ---*/
546+
i1 = 0;
547+
for (ii=0;ii<n;ii++) {
548+
if (ii == i)
549+
continue;
550+
j1 = 0;
551+
for (jj=0;jj<n;jj++) {
552+
if (jj == j)
553+
continue;
554+
c[i1][j1] = a[ii][jj];
555+
j1++;
556+
}
557+
i1++;
558+
}
559+
560+
/*--- Calculate the determinate ---*/
561+
det = MatrixDeterminant(c,n-1);
562+
563+
/*--- Fill in the elements of the cofactor ---*/
564+
b[i][j] = pow(-1.0,i+j+2.0) * det;
565+
}
566+
}
567+
for (i=0;i<n-1;i++)
568+
delete [] c[i];
569+
delete [] c;
570+
571+
}
572+
573+
void CSysMatrix::MatrixTranspose(su2double **a, unsigned long n) {
574+
575+
unsigned long i, j;
576+
su2double tmp;
577+
578+
for (i=1;i<n;i++) {
579+
for (j=0;j<i;j++) {
580+
tmp = a[i][j];
581+
a[i][j] = a[j][i];
582+
a[j][i] = tmp;
583+
}
584+
}
585+
586+
}
587+
492588
void CSysMatrix::Gauss_Elimination(unsigned long block_i, su2double* rhs, bool transposed) {
493589

494590
short iVar, jVar, kVar; // This is important, otherwise some compilers optimizations will fail
@@ -1040,13 +1136,57 @@ void CSysMatrix::InverseDiagonalBlock(unsigned long block_i, su2double *invBlock
10401136
invBlock[jVar*nVar+iVar] = aux_vector[jVar];
10411137
}
10421138

1139+
// cout <<"Gauss" <<endl;
1140+
// cout << invBlock[0] <<" "<< invBlock[1] <<" "<< invBlock[2] << endl;
1141+
// cout << invBlock[3] <<" "<< invBlock[4] <<" "<< invBlock[5] << endl;
1142+
// cout << invBlock[6] <<" "<< invBlock[7] <<" "<< invBlock[8] << endl;
1143+
//
1144+
// su2double *Block = GetBlock(block_i, block_i);
1145+
//
1146+
// Matrix = new su2double*[nVar];
1147+
// CoFactor = new su2double*[nVar];
1148+
// for (iVar=0;iVar<nVar;iVar++) {
1149+
// Matrix[iVar] = new su2double[nVar];
1150+
// CoFactor[iVar] = new su2double[nVar];
1151+
// }
1152+
//
1153+
// for (iVar = 0; iVar < nVar; iVar++) {
1154+
// for (jVar = 0; jVar < nVar; jVar++)
1155+
// Matrix[iVar][jVar] = Block[jVar*nVar+iVar];
1156+
// }
1157+
//
1158+
// Det = MatrixDeterminant(Matrix, nVar);
1159+
// MatrixCoFactor(Matrix, nVar, CoFactor);
1160+
// MatrixTranspose(CoFactor, nVar);
1161+
//
1162+
//
1163+
// for (iVar = 0; iVar < nVar; iVar++) {
1164+
// for (jVar = 0; jVar < nVar; jVar++)
1165+
// invBlock[jVar*nVar+iVar] = CoFactor[iVar][jVar]/Det;
1166+
// }
1167+
//
1168+
// for (iVar = 0; iVar < nVar; iVar++) {
1169+
// delete [] Matrix[iVar];
1170+
// delete [] CoFactor[iVar];
1171+
// }
1172+
// delete [] Matrix;
1173+
// delete [] CoFactor;
1174+
//
1175+
// cout <<"CoFactor" <<endl;
1176+
// cout << invBlock[0] <<" "<< invBlock[1] <<" "<< invBlock[2] << endl;
1177+
// cout << invBlock[3] <<" "<< invBlock[4] <<" "<< invBlock[5] << endl;
1178+
// cout << invBlock[6] <<" "<< invBlock[7] <<" "<< invBlock[8] << endl;
1179+
1180+
10431181
}
10441182

10451183

10461184
void CSysMatrix::InverseDiagonalBlock_ILUMatrix(unsigned long block_i, su2double *invBlock) {
10471185

10481186
unsigned long iVar, jVar;
1049-
1187+
su2double Det;
1188+
su2double **Matrix, **CoFactor;
1189+
10501190
for (iVar = 0; iVar < nVar; iVar++) {
10511191
for (jVar = 0; jVar < nVar; jVar++)
10521192
aux_vector[jVar] = 0.0;
@@ -1059,6 +1199,48 @@ void CSysMatrix::InverseDiagonalBlock_ILUMatrix(unsigned long block_i, su2double
10591199
invBlock[jVar*nVar+iVar] = aux_vector[jVar];
10601200
}
10611201

1202+
// cout <<"Gauss" <<endl;
1203+
// cout << invBlock[0] <<" "<< invBlock[1] <<" "<< invBlock[2] << endl;
1204+
// cout << invBlock[3] <<" "<< invBlock[4] <<" "<< invBlock[5] << endl;
1205+
// cout << invBlock[6] <<" "<< invBlock[7] <<" "<< invBlock[8] << endl;
1206+
//
1207+
// su2double *Block = GetBlock_ILUMatrix(block_i, block_i);
1208+
//
1209+
// Matrix = new su2double*[nVar];
1210+
// CoFactor = new su2double*[nVar];
1211+
// for (iVar=0;iVar<nVar;iVar++) {
1212+
// Matrix[iVar] = new su2double[nVar];
1213+
// CoFactor[iVar] = new su2double[nVar];
1214+
// }
1215+
//
1216+
// for (iVar = 0; iVar < nVar; iVar++) {
1217+
// for (jVar = 0; jVar < nVar; jVar++)
1218+
// Matrix[iVar][jVar] = Block[jVar*nVar+iVar];
1219+
// }
1220+
//
1221+
// Det = MatrixDeterminant(Matrix, nVar);
1222+
// MatrixCoFactor(Matrix, nVar, CoFactor);
1223+
// MatrixTranspose(CoFactor, nVar);
1224+
//
1225+
//
1226+
// for (iVar = 0; iVar < nVar; iVar++) {
1227+
// for (jVar = 0; jVar < nVar; jVar++)
1228+
// invBlock[jVar*nVar+iVar] = CoFactor[iVar][jVar]/Det;
1229+
// }
1230+
//
1231+
// for (iVar = 0; iVar < nVar; iVar++) {
1232+
// delete [] Matrix[iVar];
1233+
// delete [] CoFactor[iVar];
1234+
// }
1235+
// delete [] Matrix;
1236+
// delete [] CoFactor;
1237+
//
1238+
// cout <<"CoFactor" <<endl;
1239+
// cout << invBlock[0] <<" "<< invBlock[1] <<" "<< invBlock[2] << endl;
1240+
// cout << invBlock[3] <<" "<< invBlock[4] <<" "<< invBlock[5] << endl;
1241+
// cout << invBlock[6] <<" "<< invBlock[7] <<" "<< invBlock[8] << endl;
1242+
1243+
10621244
}
10631245

10641246
void CSysMatrix::BuildJacobiPreconditioner(bool transpose) {

0 commit comments

Comments
 (0)