This routine computes the solution to the system A*x = b using a preconditioned version of the Conjugate Gradient algorithm. It is assumed here that the matrix A and the preconditioner matrix M are symmetric positive definite or nearly so.
problem size, passed to F and PSOL (NEQ(1) = N).
current value of t.
array containing current dependent variable vector.
array containing current value of f(t,y).
the right hand side of the system A*x = b.
array of length N containing scale factors. 1/WGHT(i) are the diagonal elements of the diagonal scaling matrix D.
the order of the matrix A, and the lengths of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.
the maximum allowable number of iterates.
tolerance on residuals b - A*x in weighted RMS-norm.
current value of (step size h) * (coefficient l0).
preconditioner type flag.
Newton iteration counter (.ge. 0).
real work array used by routine DATP.
real work array used by preconditioner PSOL.
integer work array used by preconditioner PSOL.
the final computed approximation to the solution of the system A*x = b.
the number of iterations performed, and current order of the upper Hessenberg matrix HES.
the number of calls to PSOL.
integer error flag:
0 means convergence in LPCG iterations, LPCG .le. MAXL.
1 means the convergence test did not pass in MAXL
iterations, but the residual norm is .lt. 1,
or .lt. norm(b) if MNEWT = 0, and so X is computed.
2 means the convergence test did not pass in MAXL
iterations, residual .gt. 1, and X is undefined.
3 means there was a recoverable error in PSOL
caused by the preconditioner being out of date.
4 means there was a zero denominator in the algorithm.
The system matrix or preconditioner matrix is not
sufficiently close to being symmetric pos. definite.
-1 means there was a nonrecoverable error in PSOL.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | Neq(*) | ||||
| real(kind=dp) | :: | Tn | ||||
| real(kind=dp) | :: | Y(*) | ||||
| real(kind=dp) | :: | Savf(*) | ||||
| real(kind=dp) | :: | R(*) | ||||
| real(kind=dp) | :: | Wght(*) | ||||
| integer | :: | N | ||||
| integer, | intent(in) | :: | Maxl | |||
| real(kind=dp), | intent(in) | :: | Delta | |||
| real(kind=dp) | :: | Hl0 | ||||
| integer, | intent(in) | :: | Jpre | |||
| integer, | intent(in) | :: | Mnewt | |||
| real | :: | f | ||||
| real | :: | psol | ||||
| integer, | intent(out) | :: | Npsl | |||
| real(kind=dp), | intent(inout) | :: | X(*) | |||
| real(kind=dp), | intent(inout) | :: | P(*) | |||
| real(kind=dp), | intent(inout) | :: | W(*) | |||
| real(kind=dp), | intent(inout) | :: | Z(*) | |||
| integer, | intent(out) | :: | Lpcg | |||
| real(kind=dp), | intent(inout) | :: | Wp(*) | |||
| integer | :: | Iwp(*) | ||||
| real(kind=dp), | intent(inout) | :: | Wk(*) | |||
| integer, | intent(out) | :: | Iflag |
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| real(kind=dp), | public | :: | alpha | ||||
| real(kind=dp), | public | :: | beta | ||||
| real(kind=dp), | public | :: | bnrm | ||||
| integer, | public | :: | i | ||||
| integer, | public | :: | ier | ||||
| real(kind=dp), | public | :: | ptw | ||||
| real(kind=dp), | public | :: | rnrm | ||||
| real(kind=dp), | public | :: | ztr | ||||
| real(kind=dp), | public | :: | ztr0 |
subroutine dpcg(Neq,Tn,Y,Savf,R,Wght,N,Maxl,Delta,Hl0,Jpre,Mnewt,f,psol,Npsl,X,P,W,Z,Lpcg,Wp,Iwp,Wk,Iflag) ! integer :: Neq(*) real(kind=dp) :: Tn real(kind=dp) :: Y(*) real(kind=dp) :: Savf(*) real(kind=dp) :: R(*) real(kind=dp) :: Wght(*) integer :: N integer,intent(in) :: Maxl real(kind=dp),intent(in) :: Delta real(kind=dp) :: Hl0 integer,intent(in) :: Jpre integer,intent(in) :: Mnewt external :: f external :: psol integer,intent(out) :: Npsl real(kind=dp),intent(inout) :: X(*) real(kind=dp),intent(inout) :: P(*) real(kind=dp),intent(inout) :: W(*) real(kind=dp),intent(inout) :: Z(*) integer,intent(out) :: Lpcg real(kind=dp),intent(inout) :: Wp(*) integer :: Iwp(*) real(kind=dp),intent(inout) :: Wk(*) integer,intent(out) :: Iflag ! real(kind=dp) :: alpha , beta , bnrm , ptw , rnrm , ztr , ztr0 integer :: i , ier ! Iflag = 0 Npsl = 0 Lpcg = 0 do i = 1 , N X(i) = 0.0D0 enddo bnrm = dvnorm(N,R,Wght) ! Test for immediate return with X = 0 or X = b. ----------------------- if ( bnrm>Delta ) then ! ztr = 0.0D0 else if ( Mnewt>0 ) return !X!call dcopy(N,R,1,X,1) X(1:N) = R(1:N)!X! return endif ! Loop point for PCG iterations. --------------------------------------- INFINITE: do Lpcg = Lpcg + 1 !X!call dcopy(N,R,1,Z,1) Z(1:N) = R(1:N)!X! ier = 0 if ( Jpre/=0 ) then call psol(Neq,Tn,Y,Savf,Wk,Hl0,Wp,Iwp,Z,3,ier) Npsl = Npsl + 1 if ( ier/=0 ) then !----------------------------------------------------------------------- ! This block handles error returns from PSOL. !----------------------------------------------------------------------- if ( ier<0 ) Iflag = -1 if ( ier>0 ) Iflag = 3 return endif endif ztr0 = ztr ztr = ddot(N,Z,1,R,1) if ( Lpcg==1 ) then !X!call dcopy(N,Z,1,P,1) P(1:N) = Z(1:N)!X! elseif ( ztr0==0.0D0 ) then !----------------------------------------------------------------------- ! This block handles division by zero errors. !----------------------------------------------------------------------- Iflag = 4 return else beta = ztr/ztr0 do i = 1 , N P(i) = Z(i) + beta*P(i) enddo endif !----------------------------------------------------------------------- ! Call DATP to compute A*p and return the answer in W. !----------------------------------------------------------------------- call datp(Neq,Y,Savf,P,Wght,Hl0,Wk,f,W) ! ptw = ddot(N,P,1,W,1) if ( ptw==0.0D0 ) then Iflag = 4 exit INFINITE else alpha = ztr/ptw call daxpy(N,alpha,P,1,X,1) alpha = -alpha call daxpy(N,alpha,W,1,R,1) rnrm = dvnorm(N,R,Wght) if ( rnrm<=Delta ) return if ( Lpcg<Maxl ) cycle INFINITE Iflag = 2 if ( rnrm<=1.0D0 ) Iflag = 1 if ( rnrm<=bnrm .and. Mnewt==0 ) Iflag = 1 exit INFINITE endif enddo INFINITE end subroutine dpcg