dspiom – odepack

subroutine dspiom(Neq, Tn, Y, Savf, B, Wght, N, Maxl, Kmp, Delta, Hl0, Jpre, Mnewt, f, psol, Npsl, X, V, Hes, Ipvt, Liom, Wp, Iwp, Wk, Iflag)

This routine solves the linear system A * x = b using a scaled preconditioned version of the Incomplete Orthogonalization Method. An initial guess of x = 0 is assumed.

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).
B: the right hand side of the system Ax = b. B is also used as work space when computing the final approximation. (B is the same as V(,MAXL+1) in the call to DSPIOM.)
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, B, WGHT, and X.
MAXL: the maximum allowable order of the matrix HES.
KMP: the number of previous vectors the new vector VNEW must be made orthogonal to. KMP .le. MAXL.
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 of length N used by DATV and PSOL.
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.

V

the N by (LIOM+1) array containing the LIOM orthogonal vectors V(,1) to V(,LIOM).

HES

the LU factorization of the LIOM by LIOM upper Hessenberg matrix whose entries are the scaled inner products of AV(,k) and V(*,i).

IPVT

an integer array containg pivoting information. It is loaded in DHEFA and used in DHESL.

LIOM

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 LIOM iterations, LIOM.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.
       -1 means there was a nonrecoverable error in PSOL.

Arguments

Type	Intent	Attributes		Name
integer,		dimension(*)	::	Neq
real(kind=dp)			::	Tn
real(kind=dp),		dimension(*)	::	Y
real(kind=dp),		dimension(*)	::	Savf
real(kind=dp),	intent(inout),	dimension(*)	::	B
real(kind=dp),		dimension(*)	::	Wght
integer			::	N
integer			::	Maxl
integer			::	Kmp
real(kind=dp),	intent(inout)		::	Delta
real(kind=dp)			::	Hl0
integer			::	Jpre
integer,	intent(in)		::	Mnewt
real			::	f
real			::	psol
integer,	intent(out)		::	Npsl
real(kind=dp),		dimension(*)	::	X
real(kind=dp),		dimension(N,*)	::	V
real(kind=dp),	intent(inout),	dimension(Maxl,Maxl)	::	Hes
integer,		dimension(*)	::	Ipvt
integer,	intent(out)		::	Liom
real(kind=dp),		dimension(*)	::	Wp
integer,		dimension(*)	::	Iwp
real(kind=dp),		dimension(*)	::	Wk
integer,	intent(out)		::	Iflag

Calls

Help

Variables

Type	Visibility		Name
real(kind=dp),	public	::	bnrm
real(kind=dp),	public	::	bnrm0
integer,	public	::	i
integer,	public	::	ier
integer,	public	::	info
integer,	public	::	j
integer,	public	::	k
integer,	public	::	ll
integer,	public	::	lm1
real(kind=dp),	public	::	prod
real(kind=dp),	public	::	rho
real(kind=dp),	public	::	snormw
real(kind=dp),	public	::	tem

Subroutines

subroutine approximate()

Arguments

None

Source Code

subroutine dspiom(Neq,Tn,Y,Savf,B,Wght,N,Maxl,Kmp,Delta,Hl0,Jpre,Mnewt,f,psol,Npsl,X,V,Hes,Ipvt,Liom,Wp,Iwp,Wk,Iflag)
!
integer , dimension(*) :: Neq
real(kind=dp) :: Tn
real(kind=dp) , dimension(*) :: Y
real(kind=dp) , dimension(*) :: Savf
real(kind=dp) , intent(inout) , dimension(*) :: B
real(kind=dp) , dimension(*) :: Wght
integer :: N
integer :: Maxl
integer :: Kmp
real(kind=dp) , intent(inout) :: Delta
real(kind=dp) :: Hl0
integer :: Jpre
integer , intent(in) :: Mnewt
external :: f
!
real(kind=dp) :: bnrm , bnrm0 , prod , rho , snormw , tem
real(kind=dp) , intent(inout) , dimension(Maxl,Maxl) :: Hes
integer :: i , ier , info , j , k , ll , lm1
integer , intent(out) :: Iflag , Liom , Npsl
integer , dimension(*) :: Ipvt , Iwp
real(kind=dp) , dimension(N,*) :: V
real(kind=dp) , dimension(*) :: Wk , Wp , X
external psol
   !
   Iflag = 0
   Liom = 0
   Npsl = 0
   !-----------------------------------------------------------------------
   !  The initial residual is the vector b.  Apply scaling to b, and test
   !  for an immediate return with X = 0 or X = b.
   !-----------------------------------------------------------------------
   do i = 1 , N
      V(i,1) = B(i)*Wght(i)
   enddo
   bnrm0 = dnrm2(N,V,1)
   bnrm = bnrm0
   IFDELTA: if ( bnrm0>Delta ) then
      !  Apply inverse of left preconditioner to vector b. --------------------
      ier = 0
      if ( Jpre/=0 .and. Jpre/=2 ) then
         call psol(Neq,Tn,Y,Savf,Wk,Hl0,Wp,Iwp,B,1,ier)
         Npsl = 1
         if ( ier/=0 )then
            !  This block handles error returns forced by routine PSOL.
            if ( ier<0 ) Iflag = -1
            if ( ier>0 ) Iflag = 3
         return
         endif
         !  Calculate norm of scaled vector V(*,1) and normalize it. -------------
         do i = 1 , N
            V(i,1) = B(i)*Wght(i)
         enddo
         bnrm = dnrm2(N,V,1)
         Delta = Delta*(bnrm/bnrm0)
      endif
      tem = 1.0D0/bnrm
      call dscal(N,tem,V(1,1),1)
      !  Zero out the HES array. ----------------------------------------------
      do j = 1 , Maxl
         do i = 1 , Maxl
            Hes(i,j) = 0.0D0
         enddo
      enddo
      !-----------------------------------------------------------------------
      !  Main loop on LL = l to compute the vectors V(*,2) to V(*,MAXL).
      !  The running product PROD is needed for the convergence test.
      !-----------------------------------------------------------------------
      prod = 1.0D0
      do ll = 1 , Maxl
         Liom = ll
         !-----------------------------------------------------------------------
         !  Call routine DATV to compute VNEW = Abar*v(l), where Abar is
         !  the matrix A with scaling and inverse preconditioner factors applied.
         !  Call routine DORTHOG to orthogonalize the new vector vnew = V(*,l+1).
         !  Call routine DHEFA to update the factors of HES.
         !-----------------------------------------------------------------------
         call datv(Neq,Y,Savf,V(1,ll),Wght,X,f,psol,V(1,ll+1),Wk,Wp,Iwp,Hl0,Jpre,ier,Npsl)
         if ( ier/=0 ) then
            !  This block handles error returns forced by routine PSOL.
            if ( ier<0 ) Iflag = -1
            if ( ier>0 ) Iflag = 3
         return
         endif
         call dorthog(V(1,ll+1),V,Hes,N,ll,Maxl,Kmp,snormw)
         call dhefa(Hes,Maxl,ll,Ipvt,info,ll)
         lm1 = ll - 1
         if ( ll>1) then
             if (Ipvt(lm1)==lm1 ) prod = prod*Hes(ll,lm1)
         endif
         if ( info/=ll ) then
            !-----------------------------------------------------------------------
            !  Update RHO, the estimate of the norm of the residual b - A*x(l).
            !  test for convergence.  If passed, compute approximation x(l).
            !  If failed and l .lt. MAXL, then continue iterating.
            !-----------------------------------------------------------------------
            rho = bnrm*snormw*abs(prod/Hes(ll,ll))
            if ( rho<=Delta ) then
               call approximate()
               return
            endif
            if ( ll==Maxl ) exit
         else
            !-----------------------------------------------------------------------
            !  The last pivot in HES was found to be zero.
            !  If vnew = 0 or l = MAXL, take an error return with IFLAG = 2.
            !  otherwise, continue the iteration without a convergence test.
            !-----------------------------------------------------------------------
            if ( snormw==0.0D0 ) exit IFDELTA
            if ( ll==Maxl ) exit IFDELTA
         endif
         !  If l .lt. MAXL, store HES(l+1,l) and normalize the vector v(*,l+1).
         Hes(ll+1,ll) = snormw
         tem = 1.0D0/snormw
         call dscal(N,tem,V(1,ll+1),1)
      enddo
      !-----------------------------------------------------------------------
      !  l has reached MAXL without passing the convergence test:
      !  If RHO is not too large, compute a solution anyway and return with
      !  IFLAG = 1.  Otherwise return with IFLAG = 2.
      !-----------------------------------------------------------------------
      if ( rho<=1.0D0 ) then
         Iflag = 1
         call approximate()
         return
      elseif ( rho<=bnrm .and. Mnewt==0 ) then
         Iflag = 1
         call approximate()
         return
      endif
   elseif ( Mnewt>0 ) then
      do i = 1 , N
         X(i) = 0.0D0
      enddo
      return
   else
      !X!call dcopy(N,B,1,X,1)
      X(1:N) = B(1:N)!X!
      return
   endif IFDELTA

   Iflag = 2

contains

subroutine approximate()
!-----------------------------------------------------------------------
!  Compute the approximation x(l) to the solution.
!  Since the vector X was used as work space, and the initial guess
!  of the Newton correction is zero, X must be reset to zero.
!-----------------------------------------------------------------------
   integer :: ll
   ll = Liom
   do k = 1 , ll
      B(k) = 0.0D0
   enddo
   B(1) = bnrm
   call dhesl(Hes,Maxl,ll,Ipvt,B)
   do k = 1 , N
      X(k) = 0.0D0
   enddo
   do i = 1 , ll
      call daxpy(N,B(i),V(1,i),1,X,1)
   enddo
   do i = 1 , N
      X(i) = X(i)/Wght(i)
   enddo
   if ( Jpre<=1 ) return
   call psol(Neq,Tn,Y,Savf,Wk,Hl0,Wp,Iwp,X,2,ier)
   Npsl = Npsl + 1
   if ( ier==0 ) return
   !-----------------------------------------------------------------------
   !  This block handles error returns forced by routine PSOL.
   !-----------------------------------------------------------------------
   if ( ier<0 ) Iflag = -1
   if ( ier>0 ) Iflag = 3

end subroutine approximate

end subroutine dspiom

dspiom Subroutine

Contents

Variables

Subroutines

Source Code

subroutine dspiom(Neq, Tn, Y, Savf, B, Wght, N, Maxl, Kmp, Delta, Hl0, Jpre, Mnewt, f, psol, Npsl, X, V, Hes, Ipvt, Liom, Wp, Iwp, Wk, Iflag)

On entry

On return

Arguments

Calls

Variables

Subroutines

subroutine approximate()

Arguments

Source Code