'------------------------------------------------------ ' Conjugate Gradient Method for a Sparse Linear System '------------------------------------------------------ ' Ref.: "Numerical Recipes, Cambridge University Press ' 1986 (chapter 2.10)" [BIBLI 08]. ' ' Basic demo version by J-P Moreau, Paris. ' (www.jpmoreau.fr) '------------------------------------------------------- ' SAMPLE RUN: ' ' Solve sparse linear system (size=5): ' ' 4x1 + x2 +2x3 +0.500x4 + 2x5 = 7 ' x1 +0.5x2 = 3 ' 2x1 +3x3 = 7 ' 0.5x1 +0.625x4 = -4 ' 2x1 +16x5 = -4 ' ' ' Structure of matrix A: ' ' XXXXX ' XX ' X X ' X X ' X X ' ' Any key + Enter to continue... ' ' Number of iterations: 6 ' ' Solution vector: ' ' 2.0000 2.0000 1.0000 -8.0000 -0.5000 ' ' ' RSQ = 4.891844969899742E-010 ' ' ' Product A;X: ' ' 7.0000 3.0000 7.0000 -4.0000 -4.0000 ' '------------------------------------------------------- 'PROGRAM TSPARSE defdbl a-h,o-z defint i-n N=5 F\$="####.####" DIM A(N,N), B(N), B1(N), X(N) DIM XX(N), VV(N) 'matrix A for i=1 to N for j=1 to N A(i,j)=0# next j next i A(1,1)=4# : A(1,2)=1# : A(1,3)=2# : A(1,4)=0.5# : A(1,5)=2# A(2,1)=1# : A(2,2)=0.5# A(3,1)=2# : A(3,3)=3# A(4,1)=0.5# : A(4,4)=0.625# A(5,1)=2# : A(5,5)=16# 'vector B B(1)=7# : B(2)=3# : B(3)=7# : B(4)=-4# : B(5)=-4# 'initial guess for X for i=1 to N X(i)=0# next i cls 'clear screen print print " Structure of matrix A:" print gosub 500 'call SHOWMAT print print input " Any key to continue... ", r\$ cls : print gosub 1000 'call SPARSE(B,N,X,RSQ) print print " Solution vector:" print for i=1 to N print using F\$; X(I); next i print print print print " RSQ = "; RSQ print 'verification A.x = B for i=1 to N XX(i)=X(i) next i gosub 1100 'call ASUB(X,B1) for i=1 to N B1(i)=VV(i) next i print print " Product A.X:" print for i=1 to N print using F\$; B1(I); next i print END 'of main 'show structure of matrix (non zero elements) 500 'SUBROUTINE SHOWMAT for i=1 to N print " "; for j=1 to N if A(i,j)<>0# then print "X"; else print " "; end if next j print next i RETURN 1000 'SUBROUTINE SPARSE(B,N,X,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. '------------------------------------------------------------------- EPS=1E-6 DIM G(N),H(N),XI(N),XJ(N) EPS2=N*EPS*EPS IRST=0 1 IRST=IRST+1 for i=1 to N XX(i)=X(i) next i gosub 1100 'call ASUB(X,XI) for i=1 to N XI(i)=VV(i) next i RP=0# BSQ=0# for J=1 to N BSQ=BSQ+B(J)*B(J) XI(J)=XI(J)-B(J) RP=RP+XI(J)*XI(J) next J for i=1 to N XX(i)=XI(i) next i gosub 1200 'call ATSUB(XI,G) for i=1 to N G(i)=VV(i) next i for J=1 to N G(J)=-G(J) H(J)=G(J) next J for ITER=1 to 10*N 'main loop for i=1 to N XX(i)=H(i) next i gosub 1100 'call ASUB(H,XI) for i=1 to N XI(i)=VV(i) next i ANUM=0# ADEN=0# for J=1 to N ANUM=ANUM+G(J)*H(J) ADEN=ADEN+XI(J)*XI(J) next J if ADEN=0# then Print " Very singular matrix." end if ANUM=ANUM/ADEN for J=1 to N XI(J)=X(J) X(J)=X(J)+ANUM*H(J) next J for i=1 to N XX(i)=X(i) next i gosub 1100 'call ASUB(X,XJ) for i=1 to N XJ(i)=VV(i) next i RSQ=0# for J=1 to N XJ(J)=XJ(J)-B(J) RSQ=RSQ+XJ(J)*XJ(J) next J if RSQ=RP or RSQ<=BSQ*EPS2 then print " Number of iterations: "; ITER return 'normal return end if if RSQ > RP then for J=1 to N X(J)=XI(J) next J if IRST>=3 then return 'return if too many restarts goto 1 end if end if RP=RSQ for i=1 to N XX(i)=XJ(i) next i gosub 1200 'call ATSUB(XJ,XI) (compute gradient for next iteration) for i=1 to N XI(i)=VV(i) next i GG=0# DGG=0# for J=1 to N GG=GG+G(J)*G(J) DGG=DGG+(XI(J)+G(J))*XI(J) next J if GG=0# then print " Number of iterations: "; ITER return 'rare but normal return end if GAM=DGG/GG for J=1 to N G(J)=-XI(J) H(J)=G(J)+GAM*H(J) next J next ITER Print " Too many iterations." return 'In these versions of ASUB and ATSUB, 'we do not take any advantage of 'sparseness! (N is small). 1100 'SUBROUTINE ASUB(XX,VV) for I=1 to N VV(I)=0# for J=1 to N VV(I)=VV(I)+A(I,J)*XX(J) next J next I return 1200 'SUBROUTINE ATSUB(XX,VV) for I=1 to N VV(I)=0# for J=1 to N VV(I)=VV(I)+A(J,I)*XX(J) next J next I return 'end of file tsparse.bas 