/* WR-59 */ /* THE CONJUGATE GRADIENT METHOD */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 59-61. */ /* n = Order of the matrix A. p = Upper bound in the declaration of the arrays q, c[1:p]. For k = 1, 2, ..., p the values q[k], e[k] are stored. [see text]. b = array b[1:n] containing the constant terms of Eq. 1 [see text]. op= procedure op(n, v) res:(av); value n; integer n; arrayv, av; code computing for an arbitrary vector v[1:n] the product av[1:n] = A x v without changing v. x = In the normal case (n > 0) x is irrelevant at the beginning since the proces then starts with x = null vector. However, n may be negative which serves as a flag that an initial guess for the solution must be given. Of course, in this case abs(n) is used internally instead of n. nn = Maximum number of CG-steps to be performed. Output: x = array x[1:n] = approximate solution of (1) [see text]. nn = The number of CG-steps actually performed. If nn is equal to the originally given value, this indicates that more steps might have produced a more precise solution. q, e = array q, e[1:p] containing the values q[k], e[k] (k = 1, ..., p). If the final value of nn is below p, then these arrays are filled only up to q[nn], e[nn]. indef = Exit is used if for the first time a value q[k] becomes non-positive. In this case the matrix A is either not positive definite or so ill-conditioned that its smallest eigenvalue is in the roundoff-error level. In either case no results are given. */ CG:PROC(N,P,B,OP,X,NN,Q,E,INDEF) OPTIONS (REORDER); /*PAGE 59*/ /*00000010*/ DCL (N,P,NN,K,J,SQN,KSQN,LB,RB,LBN,UBN,LBP) BIN FIXED(31), /*00000020*/ (B,X) (*) BINARY FLOAT (21), (TR,RFE) BIN FLOAT (53), /*00000030*/ (Q,E)(*) BINARY FLOAT (21), /*00000040*/ INDEF LABEL, OP ENTRY, /*00000050*/ (RR,RAR,OLDRR, ALPHA,F,F2,EK,QK) BIN FLOAT (21), /*00000060*/ (((DR,DX,AR,R)(LBN:UBN)) CTL) BIN FLOAT(21); /*00000070*/ LBN=LBOUND(X,1); UBN =LBN+N-1; /*00000080*/ LBP=LBOUND(Q,1); /*00000090*/ ALLOCATE DR,DX,AR,R; /*00000100*/ IF N<0 THEN DO; /*00000110*/ N=-N ; /*00000120*/ CALL OP(N,X,AR); /*00000130*/ R=AR+B; /*00000140*/ END; /*00000150*/ ELSE DO; /*00000160*/ X=0E0;R=B; /*00000170*/ END; /*00000180*/ RR,RFE=0E0;KSQN,SQN=SQRT(N);ALPHA=1E0;F=EXP(1E0/N);F2=F*F; /*00000190*/ DX,DR=0E0; /*00000200*/ /* BIG LOOP FOR CG-STEP */ /*00000210*/ DO K=0 TO NN; /*00000220*/ /*TEST FOR TRUE RESIDUALS */ /*00000230*/ IF K=KSQN THEN DO; /*00000240*/ KSQN=KSQN+SQN;ALPHA=ALPHA*F;F=F*2; /*00000250*/ CALL OP(N,X,AR);RFE=0E0; /*00000260*/ DO J=LBN TO UBN; /*00000270*/ TR=(AR(J)+B(J)-R(J)); /*00000280*/ RFE=RFE+TR*TR; /*00000290*/ END; /*00000300*/ END; /*00000310*/ CALL OP(N,R,AR);OLDRR=RR; /*00000320*/ TR=0E0; /*00000330*/ DO J=LBN TO UBN; /*00000340*/ TR=TR+R(J)*R(J); /*00000350*/ END; /*00000360*/ RR=TR;TR=0E0; /*00000370*/ DO J=LBN TO UBN; /*00000380*/ TR=TR+R(J)*AR(J); /*00000390*/ END; /*00000400*/ RAR=TR; /*00000410*/ IF K=0 THEN EK=0; ELSE EK=QK*RR/OLDRR; /*00000420*/ IF K>0&K<=P THEN DO; /*00000430*/ Q(LBP+K-1)=QK;E(LBP+K-1)=EK; /*00000440*/ END; /*00000450*/ IF RR<=ALPHA*RFE THEN NN=K; /*00000460*/ IF K=NN THEN GO TO REG; /*00000470*/ QK=RAR/RR-EK; /*00000480*/ IF QK<=0E0 THEN GO TO INDEF; /*00000490*/ /*COMPUTATION OF NEW R AND X */ /*00000500*/ TR=1E0/QK; /*00000510*/ DO J=LBN TO UBN; /*00000520*/ DX(J)=(EK*DX(J)-R(J))*TR; /*00000530*/ X(J)=X(J)+DX(J); /*00000540*/ END; /*00000550*/ DO J=LBN TO UBN; /*00000560*/ DR(J)=(EK*DR(J)-AR(J))*TR; /*00000570*/ R(J)=R(J)+DR(J); /*00000580*/ END; /*00000590*/ END; /*00000600*/ REG: /*00000610*/ END CG; /*00000620*/