/*************************************************************************** * Solve a symmetric linear system using function SYMSOL * * ------------------------------------------------------------------------ * * SAMPLE RUN: * * ( Solve symmetric linear system: * * X1 -2X2 +3X3 - 4X4 + 5X5 - 6X6 = -3 * * -2X1 +3X2 -5X3 + 7X4 - 4X5 + 8X6 = 7 * * 3X1 -5X2 +5X3 + 6X4 - 2X5 + X6 = 8 * * -4X1 +7X2 +6X3 - 9X4 + 20X5 + 0.5X6 = 20.5 * * 5X1 -4X2 -2X3 + 20X4 + 11X5 +0.3333X6 = 30.3333 * * -6X1 +8X2 + X3 +0.5X4 +0.3333X5 + 13X6 = 16.8333 ) * * * * SYMMETRIC SYSTEM TO SOLVE: * * 1.000000 -2.000000 3.000000 -4.000000 5.000000 -6.000000 -3.000000 * * -2.000000 3.000000 -5.000000 7.000000 -4.000000 8.000000 7.000000 * * 3.000000 -5.000000 5.000000 6.000000 -2.000000 1.000000 8.000000 * * -4.000000 7.000000 6.000000 -9.000000 20.000000 0.500000 20.500000 * * 5.000000 -4.000000 -2.000000 20.000000 11.000000 0.333300 30.333300 * * -6.000000 8.000000 1.000000 0.500000 0.333300 13.000000 16.833300 * * * * SOLUTIONS: * * 1.00000000 * * 1.00000000 * * 1.00000000 * * 1.00000000 * * 1.00000000 * * 1.00000000 * * * * ------------------------------------------------------------------------ * * Reference: From Numath Library By Tuan Dang Trong in Fortran 77 * * [BIBLI 18]. * * * * C++ Release 1.0 By J-P Moreau, Paris * * (www.jpmoreau.fr) * **************************************************************************** To link with basis_r.cpp and vmblock.cpp. ----------------------------------------- */ #include #include REAL **A; REAL *B,*X,*D; int I,J,N; void *vmblock = NULL; void SYMSOL(int N, REAL **A, REAL *B, REAL *X, REAL *D) { /*--------------------------------------------------------------------- ! SOLVE A SYMMETRIC LINEAR SYSTEM A.X = B ! A IS SUPPOSED WELL BALANCED AND CAN BE NOT POSITIVE DEFINITE, ! THE COEFFICIENTS OF A ARE STORED IN THE UPPER TRIANGULAR HALF ! OF THE MATRIX AND ARE LOST DURING THE PROCESS. ! THE USED METHOD IS BASED UPON THE DECOMPOSITION OF THE INITIAL ! MATRIX INTO A PRODUCT OF 3 MATRICES: A = (U TRANSPOSE).D.U ! INPUTS: ! N NUMBER OF LINES OF THE LINEAR SYSTEM ! A TABLE OF DIMENSION N X N STORING THE COEFFICIENTS ! OF THE SYSTEM ! B VECTOR OF DIMENSION N (SECOND MEMBER) ! OUTPUT: ! X SYSTEM SOLUTION ! WORKING ZONE: ! D DIAGONAL D OF DECOMPOSITION !---------------------------------------------------------------------*/ // Labels: e2, e4, e6, e7, e8, e10 REAL S; int I,K,L,M; for (I=1; I<=N; I++) { L = I-1; M = I+1; S = A[I][I]; if (I==1) goto e2; for (K=1; K<=L; K++) S = S - A[K][I]*D[K]*A[K][I]; e2: D[I] = S; if (M > N) goto e6; for (J=M; J<=N; J++) { S = A[I][J]; if (I==1) goto e4; for (K=1; K<=L; K++) S = S - A[K][I]*D[K]*A[K][J]; e4: A[I][J] = S/D[I]; } } e6: for (I=1; I<=N; I++) { S = B[I]; L = I-1; if (I==1) goto e8; for (K=1; K<=L; K++) S = S - A[K][I]*X[K]; e8: X[I] = S; } for (L=1; L<=N; L++) { I = N-L+1; S = X[I]/D[I]; M = I+1; if (M>N) goto e10; for (K=M; K<=N; K++) S = S - A[I][K]*X[K]; e10: X[I] = S; } } //SYMSOL() void main() { N = 6; // size of linear system // allocate memory for matrix A and vectors B, X, D (index 0 not used) vmblock = vminit(); A = (REAL **)vmalloc(vmblock, MATRIX, N+1, N+1); X = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); B = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); D = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); if (! vmcomplete(vmblock)) LogError ("No Memory", 0, __FILE__, __LINE__); //define upper half part of matrix A A[1][1]= 1.0; A[1][2]=-2.0; A[1][3]=3.0; A[1][4]=-4.0; A[1][5]=5.0; A[1][6]=-6.0; A[2][2]= 3.0; A[2][3]=-5.0; A[2][4]=7.0; A[2][5]=-4.0; A[2][6]=8.0; A[3][3]= 5.0; A[3][4]= 6.0; A[3][5]=-2.0; A[3][6]=1.0; A[4][4]=-9.0; A[4][5]=20.0; A[4][6]=0.50; A[5][5]=11.0; A[5][6]= 0.3333; A[6][6]=13.0; //define 2nd member B[1]=-3.0; B[2]=7.0; B[3]=8.0; B[4]=20.5; B[5]=30.3333; B[6]=16.8333; // optional section to print full data for (I=2; I<=N; I++) for (J=1; J