/************************************************************* * SOLVE A LINEAR SYSTEM BY DIRECT FACTORIZATION * * ---------------------------------------------------------- * * SAMPLE RUN: * * (Solve linear system A X = B, where: * * * * 1 0 0 0 0 1 1 * * 1 1 0 0 0 -1 0 * * A = -1 1 1 0 0 1 B = 1 * * 1 -1 1 1 0 -1 0 * * -1 1 -1 1 1 1 1 * * 1 -1 1 -1 1 -1 0 ) * * * * LINEAR SYSTEM AX = B: * * * * 1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 1.0000 * * 1.0000 1.0000 0.0000 0.0000 0.0000 -1.0000 0.0000 * * -1.0000 1.0000 1.0000 0.0000 0.0000 1.0000 1.0000 * * 1.0000 -1.0000 1.0000 1.0000 0.0000 -1.0000 0.0000 * * -1.0000 1.0000 -1.0000 1.0000 1.0000 1.0000 1.0000 * * 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 0.0000 * * * * SOLUTION: * * 0.343750 * * 0.312500 * * 0.375000 * * 0.250000 * * 0.500000 * * 0.656250 * * * * ---------------------------------------------------------- * * Ref.: From Numath Library By Tuan Dang Trong in fortran 77 * * [BIBLI 18]. * * * * C++ Release By J-P Moreau, Paris * * (www.jpmoreau.fr) * *************************************************************/ #include #include #define NMAX 50 typedef double MAT[NMAX][NMAX]; MAT A, W; double B[NMAX], X[NMAX], Z[NMAX]; int I,J,N; void DLITTL(MAT A, double *B, int N, double *X, MAT W, double *Z) { /*---------------------------------------------------------------------- ! SOLVE A LINEAR SYSTEM BY DIRECT FACTORIZATION ! A*X = L*U*X = B ! L(LOWER TRIANGULAR MATRIX WITH L(I,I)=1) ! U(UPPER TRIANGULAR MATRIX),THE SUB-DIAGONAL TERMS OF L AND THE ! TERMS U(I,J) ARE STORED IN THE SAME MATRIX, W. THE RESOLUTION ! IS OBTAINED BY SOLVING TWO TRIANGULAR SYSTEMS: ! L*Z = B AND U*X = Z. ! THE PARTIAL PIVOTING IS MADE BY CHOOSING THE BIGGEST ELEMENT OF ! EACH COLUMN OF THE TRANSfor (MATION MATRIX. ! INPUTS: ! A = TABLE OF SIZE (N,N) R*8 ! B = SECOND MEMBER VECTOR (N) R*8 ! N = ORDER OF LINEAR SYSTEM I*4 ! OUTPUTS: ! X = SOLUTION VECTOR (N) R*8 ! WORKING ZONE: ! W = TABLE OF SIZE (LDA,N) R*8 ! Z = AUXILIARY VECTOR (N) R*8 ! NOTE: ! MESSAGE '** DLITTL ** NO UNIQUE SOLUTION' IS GIVEN WHEN A ! NULL PIVOT IS FOUND. !----------------------------------------------------------------------*/ // Labels: e2,e10,e32 double SUM,AMAX,T; int I,J,K,P; P=1; AMAX=fabs(A[1][1]); // SEEK MAX. ELEMENT OF FIRST COLUMN for (J=2; J<=N; J++) { if (fabs(A[J][1]) < AMAX) goto e2; AMAX=A[J][1]; P=J; e2: ;} if (AMAX == 0.0) goto e32; if (P != 1) { // EXCHANGE LINES OF MATRIX A AND ORDER OF UNKNOWNS // LINKED TO B for (J=1; J<=N; J++) { T=A[1][J]; A[1][J]=A[P][J]; A[P][J]=T; } T=B[1]; B[1]=B[P]; B[P]=T; } // FIRST LINE OF U AND FIRST COLUMN OF L W[1][1]=A[1][1]; for (J=2; J<=N; J++) { W[1][J]=A[1][J]; W[J][1]=A[J][1]/W[1][1]; } // SECOND LINE OF U AND SECOND COLUMN OF L // (N-1)TH LINE OF U AND (N-1)TH COLUMN OF L for (I=2; I0; I--) { SUM=0.0; for (K=I+1; K<=N; K++) SUM += W[I][K]*X[K]; X[I]=(Z[I]-SUM)/W[I][I]; } // END OF RESOLUTION OF A*X = L*(U*X) = B return; e32: printf("\n ** DLITTL ** NO UNIQUE SOLUTION.\n\n"); } void main() { N=6; //size of system for (I=1; I<=N; I++) { for (J=1; J<=N; J++) A[I][J]=0.0; B[I]=0.0; } // define elements of A <> 0 A[1][1]= 1.0; A[1][6]= 1.0; A[2][1]= 1.0; A[2][2]= 1.0; A[2][6]=-1.0; A[3][1]=-1.0; A[3][2]= 1.0; A[3][3]= 1.0; A[3][6]= 1.0; A[4][1]= 1.0; A[4][2]=-1.0; A[4][3]= 1.0; A[4][4]= 1.0; A[4][6]=-1.0; A[5][1]=-1.0; A[5][2]= 1.0; A[5][3]=-1.0; A[5][4]= 1.0; A[5][5]= 1.0; A[5][6]= 1.0; A[6][1]= 1.0; A[6][2]=-1.0; A[6][3]= 1.0; A[6][4]=-1.0; A[6][5]= 1.0; A[6][6]=-1.0; //define elements of B <> 0 B[1]=1.0; B[3]=1.0; B[5]=1.0; // print data printf("\n LINEAR SYSTEM AX=B:\n\n"); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) printf("%8.4f", A[I][J]); printf(" %8.4f\n", B[I]); } // call main function DLITTL(A,B,N,X,W,Z); //print results printf("\n SOLUTION:\n\n"); for (I=1; I<=N; I++) printf(" %f\n", X[I]); printf("\n"); } // End of file tdlittl.cpp