'************************************************************************ '* This is a test program for the subroutine 1000 to solve a linear * '* system * '* A * X = Y * '* for a symmetric positive definite matrix A using the conjugate * '* gradient method. * '* -------------------------------------------------------------------- * '* Ref.: "Numerical algorithms with C, by Gisela Engeln-Muellges and * '* Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* -------------------------------------------------------------------- * * '* SAMPLE RUN: * '* * '* Solve the following linear system: * '* * '* 4 size N of matrix * '* 5.0 -1.0 -1.0 -1.0 the upper triangle (with diagonal) * '* 5.0 -1.0 -1.0 of the positive definite matrix A * '* 5.0 -1.0 * '* 5.0 * '* 1.0 right hand side of the linear system * '* 1.0 * '* 1.0 * '* 1.0 * '* * '* ----------------------------------------------------- * '* Symmetric Linear System (Conjugate Gradient Method) * '* ----------------------------------------------------- * '* Test data (Matrix and right hand side): * '* 5.00000 -1.00000 -1.00000 -1.00000 1.00000 * '* -1.00000 5.00000 -1.00000 -1.00000 1.00000 * '* -1.00000 -1.00000 5.00000 -1.00000 1.00000 * '* -1.00000 -1.00000 -1.00000 5.00000 1.00000 * '* * '* Solution for the linear system: * '* 0.50000 0.50000 0.50000 0.50000 * '* ----------------------------------------------------- * '* * '* Basic version without dynamic allocations by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '************************************************************************ defdbl a-h,o-z defint i-n ISIZE=25 ZERO=0 dim a(ISIZE,ISIZE) 'the upper triangle of a positive definite real 'matrix. dim y(ISIZE) 'right hand side of the linear system dim x(ISIZE) 'solution vector of the system 'integer n size of matrix A 'integer fehler% return value of cg_method() ' ----------------------------- print header --------------------- cls print print "-----------------------------------------------------" print " Symmetric Linear System (Conjugate Gradient Method) " print "-----------------------------------------------------" ' ----------------------------- Set inputs ----------------------- n=4 'size of system to solve 'define a matrix a(0,0)=5: a(0,1)=-1: a(0,2)=-1: a(0,3)=-1 a(1,1)= 5: a(1,2)=-1: a(1,3)=-1 a(2,2)= 5: a(2,3)=-1 a(3,3)= 5 'define b vector y(0)=1: y(1)=1: y(2)=1: y(3)=1 ' -------------- Print input data as a safeguard ---------------- F$="###.#####" print " Test data (Matrix and right hand side):" for i = 0 to n-1 for j=0 to i-1 print using F$; a(j,i); next for j=i to n-1 print using F$; a(i,j); next print using F$; y(i) next ' ---------------- solve linear system --------------------------- 'call cg_method(n, a, y, x, fehler) perform CG method gosub 1000 ' -------------------- Print output results ---------------------- if fehler%<>0 then 'Error in CG method? print " Error in executing CG method." print " Error code: "; fehler% else print " Solution for the linear system:" for i=0 to n-1 print using F$; x(i); next print end if print "-----------------------------------------------------" print print end 'Subroutine cg_method ( Conjugate Gradient Method ' n, Size of the linear system ' a, System matrix ' y, right hand side ' x, solution vector ' fehler% error code ' ) ' original name: cg_verfahren() 1000 MACHEPS=2e-16 '************************************************************************ '* CG solves the linear system * '* A * X = Y * '* for a symmetric, positive definite matrix A via the conjugate * '* gradient method. * '* * '* Input parameters: * '* ================= * '* n Size of the linear system * '* a [0..n-1,0..n-1] system matrix A. Only the upper triangle of A is * '* used. * '* y [0..n-1] vector of the right hand side * '* * '* Output parameters: * '* ================== * '* x [0..n-1] vector giving the solution * '* * '* Return value: * '* ============= * '* = 0: all is ok * '* = 1: n < 2 or other disallowed input parameters * '* = 2: memory exceeded * '* * '************************************************************************ dim d(n) ' (0..n-1) auxiliary vectors d and g dim g(n) dim AmalD(n) ' (0..n-1) auxiliary vector A * d ' ' main variables: ' ' real alpha, coefficient ' beta, coefficient ' dividend, numerator and denominator of a fraction ' divisor, respectively, used to compute alpha, beta ' hilf, auxiliary variables ' hilf2, ' abstand, distance of two successive approximations ' for the solution vector x (taken in the ' euclidean norm) ' xnorm euklidean norm of x ' ' integer k, i, j loop variables if n < 2 then ' invalid parameter? fehler%=1 return end if '------------------------------------------------------------------ ' start with x at the origin '------------------------------------------------------------------ for i = n - 1 to 0 step -1 x(i) = ZERO next '------------------------------------------------------------------ ' initialize d and g : ' d = -g = -(a*x - y) = y (since x = 0) '------------------------------------------------------------------ for i = n - 1 to 0 step -1 hilf = y(i) d(i) = hilf g(i) = -hilf next '------------------------------------------------------------------ ' perform at most n steps of the CG Method '------------------------------------------------------------------ for k = n to 0 step -1 '---------------------------------------------------------------- ' compute new alpha: ' alpha = -(d(transp) * g) / (d(transp) * (a * d)) '---------------------------------------------------------------- dividend = ZERO divisor = ZERO for i = n - 1 to 0 step -1 dividend = dividend + d(i) * g(i) hilf = ZERO for j = 0 to i-1 hilf = hilf + a(j,i) * d(j) next j for j = i to n-1 hilf = hilf + a(i,j) * d(j) next j AmalD(i) = hilf divisor = divisor + d(i) * hilf next i if divisor = ZERO then fehler%=0 return end if alpha = -dividend / divisor '---------------------------------------------------------------- ' compute the norm of x und alpha * d and find a new x: ' x = x + alpha * d, then check whether x is close enough, ' in order to stop the process before n complete steps '---------------------------------------------------------------- xnorm = ZERO abstand = ZERO for i = n - 1 to 0 step -1 hilf = x(i) xnorm = xnorm + hilf*hilf hilf2 = alpha * d(i) abstand = abstand + hilf2*hilf2 x(i) = hilf + hilf2 next i if abstand < MACHEPS * xnorm then fehler%=0 return end if '---------------------------------------------------------------- ' compute new g: g = g + alpha * (a * d) '---------------------------------------------------------------- for i = n - 1 to 0 step -1 g(i) = g(i) + alpha * AmalD(i) next '---------------------------------------------------------------- ' compute new beta: ' beta = (g(transp) * (a * d)) / (d(transp) * (a * d)) '---------------------------------------------------------------- dividend = ZERO for i = n - 1 to 0 step -1 dividend = dividend + g(i) * AmalD(i) next beta = dividend / divisor '---------------------------------------------------------------- ' compute new d : d = - g + beta * d '---------------------------------------------------------------- for i = n - 1 to 0 step -1 d(i) = -g(i) + beta * d(i) next i next k 'end k loop fehler%=0 return ' end of file cgtst1.bas