dpcgs Subroutine

subroutine dpcgs(Neq, Tn, Y, Savf, R, Wght, N, Maxl, Delta, Hl0, Jpre, Mnewt, f, psol, Npsl, X, P, W, Z, Lpcg, Wp, Iwp, Wk, Iflag)

This routine computes the solution to the system A*x = b using a scaled preconditioned version of the Conjugate Gradient algorithm.

It is assumed here that the scaled matrix D-1 * A * D and the scaled preconditioner D-1 * M * D are close to being symmetric positive definite.

On entry

NEQ

problem size, passed to F and PSOL (NEQ(1) = N).

TN

current value of t.

Y

array containing current dependent variable vector.

SAVF

array containing current value of f(t,y).

R

the right hand side of the system A*x = b.

WGHT

array of length N containing scale factors. 1/WGHT(i) are the diagonal elements of the diagonal scaling matrix D.

N

the order of the matrix A, and the lengths of the vectors Y, SAVF, R, WGHT, P, W, Z, WK, and X.

MAXL

the maximum allowable number of iterates.

DELTA

tolerance on residuals b - A*x in weighted RMS-norm.

HL0

current value of (step size h) * (coefficient l0).

JPRE

preconditioner type flag.

MNEWT

Newton iteration counter (.ge. 0).

WK

real work array used by routine DATP.

WP

real work array used by preconditioner PSOL.

IWP

integer work array used by preconditioner PSOL.

On return

X

the final computed approximation to the solution of the system A*x = b.

LPCG

the number of iterations performed, and current order of the upper Hessenberg matrix HES.

NPSL

the number of calls to PSOL.

IFLAG

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 scaled matrix or scaled preconditioner is not
          sufficiently close to being symmetric pos. definite.
       -1 means there was a nonrecoverable error in PSOL.

Arguments

Type IntentOptional Attributes Name
integer, dimension(*) :: Neq
real(kind=dp) :: Tn
real(kind=dp), dimension(*) :: Y
real(kind=dp), dimension(*) :: Savf
real(kind=dp), dimension(*) :: R
real(kind=dp), dimension(*) :: 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), dimension(*) :: X
real(kind=dp), intent(inout), dimension(*) :: P
real(kind=dp), intent(inout), dimension(*) :: W
real(kind=dp), intent(inout), dimension(*) :: Z
integer, intent(out) :: Lpcg
real(kind=dp), intent(inout), dimension(*) :: Wp
integer, dimension(*) :: Iwp
real(kind=dp), intent(inout), dimension(*) :: Wk
integer, intent(out) :: Iflag

Calls

proc~~dpcgs~~CallsGraph proc~dpcgs dpcgs.inc::dpcgs datp datp proc~dpcgs->datp daxpy daxpy proc~dpcgs->daxpy dvnorm dvnorm proc~dpcgs->dvnorm

Variables

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

Source Code

subroutine dpcgs(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 , dimension(*) :: Neq
real(kind=dp) :: Tn
real(kind=dp) , dimension(*) :: Y
real(kind=dp) , dimension(*) :: Savf
real(kind=dp) , dimension(*) :: R
real(kind=dp) , dimension(*) :: 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
!
real(kind=dp) :: alpha , beta , bnrm , ptw , rnrm , ztr , ztr0
integer :: i , ier
integer , intent(out) :: Iflag , Lpcg , Npsl
integer , dimension(*) :: Iwp
real(kind=dp) , intent(inout) , dimension(*) :: P , W , Wk , Wp , X , Z
!
external psol
   !
   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!
      X(1:N) = R(1:N)
      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 = 0.0D0
      do i = 1 , N
         ztr = ztr + Z(i)*R(i)*Wght(i)**2
      enddo
      if ( Lpcg/=1 ) then
         if ( ztr0==0.0D0 ) then
            !-----------------------------------------------------------------------
            !  This block handles division by zero errors.
            !-----------------------------------------------------------------------
            Iflag = 4
            return
         endif
         beta = ztr/ztr0
         do i = 1 , N
            P(i) = Z(i) + beta*P(i)
         enddo
      else
         !X!call dcopy(N,Z,1,P,1)
         P(1:N) = Z(1:N) !X!
      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 = 0.0D0
      do i = 1 , N
         ptw = ptw + P(i)*W(i)*Wght(i)**2
      enddo
      if ( ptw==0.0D0 ) then
         !-----------------------------------------------------------------------
         !  This block handles division by zero errors.
         !-----------------------------------------------------------------------
         Iflag = 4
         return
      endif
      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 ) exit INFINITE
   enddo INFINITE
   Iflag = 2
   if ( rnrm<=1.0D0 ) Iflag = 1
   if ( rnrm<=bnrm .and. Mnewt==0 ) Iflag = 1
end subroutine dpcgs