/**************************************************** * Calculate the determinant of a real square matrix * * by the LU decomposition method * * ------------------------------------------------- * * SAMPLE RUN: * * (Calculate the determinant of real square matrix: * * | 1 0.5 0 0 0 0 0 0 0 | * * |-1 0.5 1 0 0 0 0 0 0 | * * | 0 -1 0 2.5 0 0 0 0 0 | * * | 0 0 -1 -1.5 4 0 0 0 0 | * * | 0 0 0 -1 -3 5.5 0 0 0 | * * | 0 0 0 0 -1 -4.5 7 0 0 | * * | 0 0 0 0 0 -1 -6 8.5 0 | * * | 0 0 0 0 0 0 -1 -7.5 10 | * * | 0 0 0 0 0 0 0 -1 -9 | ) * * * * * * Determinant = 1.000000 * * * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***************************************************** Exact value is: 1 ---------------------------------------------------*/ #include #include #include "lu.h" int main() { REAL **A; // input matrix of size n x n int *INDX; // dummy integer vector void *vmblock = NULL; // NOTA: index zero not used here. int i, id, j, n, rc; REAL det; n=9; vmblock = vminit(); A = (REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); INDX = (int *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; //not enough memory } for (i=0; i<=n; i++) for (j=0; j<=n; j++) A[i][j] = 0.0; // define matrix A column by column A[1][1]=1.0; A[2][1]=-1.0; A[1][2]=0.5; A[2][2]=0.5; A[3][2]=-1.0; A[2][3]=1.0; A[4][3]=-1.0; A[3][4]=2.5; A[4][4]=-1.5; A[5][4]=-1.0; A[4][5]=4.0; A[5][5]=-3.0; A[6][5]=-1.0; A[5][6]=5.5; A[6][6]=-4.5; A[7][6]=-1.0; A[6][7]=7.0; A[7][7]=-6.0; A[8][7]=-1.0; A[7][8]=8.5; A[8][8]=-7.5; A[9][8]=-1.0; A[8][9]=10.0; A[9][9]=-9.0; //call LU decomposition rc = LUDCMP(A,n,INDX,&id); //calculate determinant and display det = id; if (rc==0) { for (i=1; i<=n; i++) det *= A[i][i]; printf("\n Determinant = %f\n\n", det); return 0; } else { printf("\n Error in LU decomposition.\n\n"); return 2; } } // main // end of file deter1.cpp