/*----------------------------------------------------- ! Conjugate Gradient Method for a Sparse Linear System !------------------------------------------------------ ! Ref.: "Numerical Recipes, Cambridge University Press ! 1986 (chapter 2.10)" [BIBLI 08]. ! ! C++ demo1 version by J-P Moreau, Paris. ! (www.jpmoreau.fr) !------------------------------------------------------- ! SAMPLE RUN: ! ! Solve sparse linear system (size=10): ! ! 2x1 -x10 = 0 ! 2x2 -x9 -x10 = 0 ! 2x3 -x8 -x9 = 0 ! 2x4 -x7 -x6 = 0 ! 2x5 -x6 -x7 = 0 ! -x5 +2x6 = 11 ! -x4 -x5 +2x7 = 0 ! -x3 -x4 +2x8 = 0 ! -x2 -x3 +2x9 = 0 ! -x1 -x2 +2x10 = 0 ! ! ! Structure of sparse matrix A: ! ! X X ! X XX ! X XX ! X XX ! XXX ! XX ! XX X ! XX X ! XX X ! XX X ! ! Any key to continue... ! ! Number of iterations: 10 ! ! Solution vector: ! ! 1.0000 3.0000 5.0000 7.0000 9.0000 ! 10.0000 8.0000 6.0000 4.0000 2.0000 ! ! RSQ = 1.386462e-018 ! ! ! Product A.X: ! ! 0.0000 0.0000 0.0000 0.0000 0.0000 ! 11.0000 0.0000 0.0000 0.0000 0.0000 ! !-----------------------------------------------------*/ #include #include #define SIZE 11 #define EPS 1e-6 //global variables double A[SIZE][SIZE]; double B[SIZE], B1[SIZE], X[SIZE]; double RSQ; int i,j,N; void ASUB(double *,double *); void ATSUB(double *,double *); void SPARSE(double *B,int M,double *X,double *RSQ) { /*------------------------------------------------------------------ !Solves the linear system A.x = b for the vector X of length N, !given the right-hand vector B, and given two subroutines, !ASUB(XIN,XOUT) and ATSUB(XIN,XOUT), which respectively calculate !A.x and AT.x (AT for Transpose of A) for x given as first argument, !returning the result in their second argument. !These subroutines should take every advantage of the sparseness !of the matrix A. On input, X should be set to a first guess of the !desired solution (all zero components is fine). On output, X is !the solution vector, and RSQ is the sum of the squares of the !components of the residual vector A.x - b. If this is not small, !then the matrix is numerically singular and the solution represents !a least-squares approximation. !-----------------------------------------------------------------*/ //LABEL: L1 double G[SIZE],H[SIZE],XI[SIZE],XJ[SIZE]; double EPS2,BSQ,RP,ANUM,ADEN,GAM,GG,DGG; int IRST,ITER,J; EPS2=N*EPS*EPS; IRST=0; L1:IRST++; ASUB(X,XI); RP=0.0; BSQ=0.0; for (J=1; J<=M; J++) { BSQ=BSQ+B[J]*B[J]; XI[J]=XI[J]-B[J]; RP=RP+XI[J]*XI[J]; } ATSUB(XI,G); for (J=1; J<=M; J++) { G[J]=-G[J]; H[J]=G[J]; } for (ITER=1; ITER<=10*M; ITER++) { //main loop ASUB(H,XI); ANUM=0.0; ADEN=0.0; for (J=1; J<=M; J++) { ANUM=ANUM+G[J]*H[J]; ADEN=ADEN+XI[J]*XI[J]; } if (ADEN==0.0) printf(" Very singular matrix.\n"); ANUM=ANUM/ADEN; for (J=1; J<=M; J++) { XI[J]=X[J]; X[J]+=ANUM*H[J]; } ASUB(X,XJ); *RSQ=0.0; for (J=1; J<=M; J++) { XJ[J]=XJ[J]-B[J]; *RSQ=*RSQ+XJ[J]*XJ[J]; } if (*RSQ==RP || *RSQ<=BSQ*EPS2) { printf(" Number of iterations: %d\n",ITER); return; //normal return } if (*RSQ > RP) { for (J=1; J<=M; J++) X[J]=XI[J]; if (IRST>=3) return; //return if too many restarts goto L1; } RP=*RSQ; ATSUB(XJ,XI); //compute gradient for next iteration GG=0.0; DGG=0.0; for (J=1; J<=M; J++) { GG=GG+G[J]*G[J]; DGG=DGG+(XI[J]+G[J])*XI[J]; } if (GG==0.0) { printf(" Number of iterations: %d\n", ITER); return; //rare but normal return } GAM=DGG/GG; for (J=1; J<=M; J++) { G[J]=-XI[J]; H[J]=G[J]+GAM*H[J]; } } //main loop printf(" Too many iterations.\n\n"); } /*In these versions of ASUB1 and ATSUB1, we do not take any advantage of sparseness! On the contrary, ASUB and ATSUB take into account the sparseness of the matrix A. N.B. For big values of N, the difference in computation time may be important. */ void ASUB1(double *X,double *V) { int I,J; for (I=1; I<=N; I++) { V[I]=0.0; for (J=1; J<=N; J++) V[I]=V[I]+A[I][J]*X[J]; } } void ASUB(double *X,double *V) { int I; V[1]=A[1][1]*X[1]+A[1][10]*X[10]; for (I=2; I<6; I++) V[I]=A[I][I]*X[I]+A[I][N-I+1]*X[N-I+1]+A[I][N-I+2]*X[N-I+2]; V[6]=A[6][6]*X[6]+A[6][5]*X[5]; for (I=7; I<11; I++) V[I]=A[I][I]*X[I]+A[I][N-I+1]*X[N-I+1]+A[I][N-I+2]*X[N-I+2]; } void ATSUB1(double *X,double *V) { int I,J; for (I=1; I<=N; I++) { V[I]=0.0; for (J=1; J<=N; J++) V[I]=V[I]+A[J][I]*X[J]; } } void ATSUB(double *X,double *V) { int I; V[1]=A[1][1]*X[1]+A[10][1]*X[10]; for (I=2; I<6; I++) V[I]=A[I][I]*X[I]+A[N-I+1][I]*X[N-I+1]+A[N-I+2][I]*X[N-I+2]; V[6]=A[6][6]*X[6]+A[5][6]*X[5]; for (I=7; I<11; I++) V[I]=A[I][I]*X[I]+A[N-I+1][I]*X[N-I+1]+A[N-I+2][I]*X[N-I+2]; } // show structure of sparse matrix A void ShowMat() { int i,j; printf("\n"); for (i=1; i<=N; i++) { printf(" "); for (j=1; j<=N; j++) if (A[i][j]!=0.0) printf("X"); else printf(" "); printf("\n"); } printf("\n"); } void main() { // matrix A N=10; for (i=1; i<=N; i++) for (j=1; j<=N; j++) A[i][j]=0.0; A[1][1]=2.0; A[1][10]=-1.0; A[2][2]=2.0; A[2][9]=-1.0; A[2][10]=-1.0; A[3][3]=2.0; A[3][8]=-1.0; A[3][9]=-1.0; A[4][4]=2.0; A[4][7]=-1.0; A[4][8]=-1.0; A[5][5]=2.0; A[5][6]=-1.0; A[5][7]=-1.0; A[6][5]=-1.0; A[6][6]=2.0; A[7][4]=-1.0; A[7][5]=-1.0; A[7][7]=2.0; A[8][3]=-1.0; A[8][4]=-1.0; A[8][8]=2.0; A[9][2]=-1.0; A[9][3]=-1.0; A[9][9]=2.0; A[10][1]=-1.0; A[10][2]=-1.0; A[10][10]=2.0; // vector B for (i=1; i<=N; i++) B[i]=0.0; B[6]=11.0; // initial guess for X for (i=1; i<=N; i++) X[i]=0.0; printf("\n Structure of sparse matrix A:\n"); ShowMat(); printf(" Any key to continue... "); getchar(); printf("\n"); SPARSE(B,N,X,&RSQ); printf("\n Solution vector:\n\n"); for (i=1; i<6; i++) printf("%8.4f",X[i]); printf("\n"); for (i=6; i<=N; i++) printf("%8.4f",X[i]); printf("\n"); printf("\n RSQ = %e\n\n", RSQ); // verification A.x = B ASUB(X,B1); printf("\n Product A.X:\n\n"); for (i=1; i<6; i++) printf("%8.4f",B1[i]); printf("\n"); for (i=6; i<=N; i++) printf("%8.4f",B1[i]); printf("\n\n"); } // end of file tsparse1.cpp