/****************************************************** * Solving a linear system AX = B by LU decomposition * * * * C++ version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * --------------------------------------------------- * * SAMPLE RUN: * * * * Input file (test_lu.dat): * * * * 4 * * 8 2 3 12 25.0 * * 2 4 7 0.25 13.25 * * 3 7 3 5 18.0 * * 12 0.25 5 2 19.25 * * * * Output file (test_lu.lst): * * * * --------------------------------------------------- * * LINEAR SYSTEM TO BE SOLVED: * * --------------------------------------------------- * * N=4 * * * * 8.000000 2.000000 3.000000 12.00000 25.00000 * * 2.000000 4.000000 7.000000 0.250000 13.25000 * * 3.000000 7.000000 3.000000 5.000000 18.00000 * * 12.00000 0.250000 5.000000 2.000000 19.25000 * * * * System solution: * * * * X1= 1.000000 * * X2= 1.000000 * * X3= 1.000000 * * X4= 1.000000 * * --------------------------------------------------- * * Uses: basis_r.cpp, lu.cpp, wmblock.cpp * ******************************************************/ #include #include //Functions of lu.cpp called by main int LUDCMP(REAL **, int, int *, int *); void LUBKSB(REAL **, int, int *, REAL *); int main() { REAL **A; // matrice n+1 x n+2 REAL *B; // vecteur n+1 REAL *temp; // vecteur n+1 int *INDX; // vecteur n+1 void *vmblock = NULL; // NOTA: index zero not used here. int d=0, i, j, n, rc=0; char input[20], output[20], s[20]; FILE *fp1, *fp2; strcpy(s,"Test_lu"); strcpy(input,s); strcat(input,".dat"); strcpy(output,s); strcat(output,".lst"); fp1 = fopen(input,"r"); fp2 = fopen(output,"w"); fscanf(fp1,"%d",&n); // size of given system AX=B //dynamic allocations of real matrix A, real vector B, //temporary vector temp and integer vector INDX. vmblock = vminit(); A = (REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); B = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); temp = (REAL *) vmalloc(vmblock, VEKTOR, n+2, 0); INDX = (int *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; } WriteHead(fp2," LINEAR SYSTEM TO BE SOLVED:"); fprintf(fp2," N=%d\n\n",n); for (i=1; i