Nonstationary.F
#include "parmdefines.h"
#include "iadefines.h"
c***********************************************************************
#include "author.inc"
c* $Id: Nonstationary.F,v 1.18 1996/04/24 19:34:42 turner Exp $
c*
c* Wrapper for routines that compute the solution to a system of
c* linear equations of the form Ax=b by one of several nonstationary
c* iterative methods.
c*
c* <PARAMETER LIST> (not in same order as arg list)
c*
c* Input:
c* mode - operational mode
c* 0 ==> black-box mode
c* Control parameters (iparm and rparm values) are set
c* by the routine to reasonable defaults, depending on
c* method. One parameter, however, must still be set.
c* That is the maximum number of iterations to allow:
c* itmax - iparm(_JT_itmax_)
c* 1 ==> hybrid, or partial black-box, mode
c* Control parameters related to output, and overall
c* convergence criterion, are set by user. Remainder
c* are set automatically by the routine. Parameters
c* that can be set are:
c* iparm(_JT_out_)
c* iparm(_JT_luout_)
c* iparm(_JT_luerr_)
c* iparm(_JT_itmax_)
c* iparm(_JT_stop_) (see Note)
c* iparm(_JT_pre_)
c* rparm(_JT_eps_) (see Note)
c* See the discussion of the elements of iparm and rparm
c* below for more info.
c*
c* Note: Defaults can be used for iparm(_JT_stop_) and
c* rparm(_JT_eps_) by by setting them to zero.
c* 2 ==> full manual mode
c* Iparm and rparm are assumed to be set.
c*
c* b - source vector
c* a - coefficient matrix
c* ia - integer vector containing info about how "a" is stored
c* NOTE: see description of ia below
c* ja - integer map for coefficient
c*
c* In/Out:
c* method - indicates method to use
c* 0 ==> unspecified (black-box)
c* 1 ==> standard conjugate gradients (CG)
c* 2 ==> transpose-free quasi-minimal residuals method (TFQMR)
c* 3 ==> generalized minimal residual (GMRES)
c*
c* NOTES:
c* - If a method is not specified (method=0), the routine tries
c* a sequence of methods and options, from fastest to most
c* likely to succeed. The sequence is:
c* TFQMR -> GMRES(20) -> GMRES(40) -> GMRES(80)
c* - Some control parameters are modified in an attempt to achieve
c* convergence if methods start failing.
c* - CG is for symmetric systems only.
c* - TFQMR is a fast, fairly stable method. Residual norms
c* decrease nearly monotonically, but breakdowns can occur.
c* - GMRES is the safest of these nonsymmetric methods, but also
c* typically the slowest. Breakdowns should not occur, and
c* although the method can fail to converge, increasing iparm(_JT_nold_)
c* (see below) increases the likelihood of convergence (and
c* also cost, both in CPU time and memory required).
c* - Routine automatically enters black-box mode if chosen solver
c* fails (even if chosen solver was CG).
c* - On return the value of method will indicate the final method
c* used. Hence it should be reset prior to a subsequent call.
c*
c* iparm - array of integer parameters
c* rparm - array of floating point parameters
c* NOTE: see description of iparm and rparm below
c* x - solution vector (whatever is in x on entry is used as an
c* initial guess)
c*
c* Output:
c* rnormt - norm of true residual
c* errt - error estimate based on true residual
c* rnorm - norm of residual at each iteration
c* err - error estimate at each iteration
c* ap - preconditioner
c* iap - integer vector containing info about how "ap" is stored
c* NOTE: see description of ia below
c* jap - same as ja, but for preconditioner
c* status - return status:
c* -5 ==> iteration failed due to breakdown
c* -4 ==> itmax exceeded
c* -3 ==> internal error
c* -2 ==> memory allocation failure
c* -1 ==> invalid argument(s)
c* 0 ==> convergence was achieved
c* >0 ==> preconditioning needed
c*
#include "parmdesc.inc"
c*
#include "iadesc.inc"
c*
c* Variables in iap are analogous to those in ia.
c*
c* <SUBROUTINES REQUIRED>
c*
c* JT_BCGS
c* JT_CG
c* JT_GMRES
c* JT_SetDefaults
c* JT_StoppingTest
c*
c* <FUNCTIONS REQUIRED>
c*
c* JT_VectorNorm
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Nonstationary (method, mode,
& b, a, ia, ja, iparm, rparm, x, cpu,
& rnormt, errt, rnorm, err, ap, iap, jap, status)
implicit none
c
c ... Input:
integer mode
integer ia(_JT_no_of_storage_parameters_), ja(*)
real a(*)
#ifdef strict_f77
real b(*)
#else
real b(ia(_JT_nrows_))
#endif
c
c ... In/Out:
integer method
integer iap(_JT_no_of_storage_parameters_), jap(*)
integer iparm(_JT_no_of_iparms_)
real rparm(_JT_no_of_rparms_)
#ifdef strict_f77
real x(*)
#else
real x(ia(_JT_nrows_))
#endif
real ap(*)
c
c ... Output:
integer status
#ifdef strict_f77
real cpu(0:*), rnorm(0:*), err(0:*), rnormt(0:*), errt(0:*)
#else
real cpu(0:iparm(_JT_itmax_))
real rnorm(0:iparm(_JT_itmax_)), err(0:iparm(_JT_itmax_))
real rnormt(0:iparm(_JT_itmax_)), errt(0:iparm(_JT_itmax_))
#endif
c
c ... Local:
real JT_VectorNorm, zero
real bnorm, err_end, rnormt_end, errt_end
real user_out, user_luout, user_luerr, user_itmax,
& user_stop, user_pre, user_eps
c
parameter (zero=0.0d0)
c
c ... Ensure that all variables are saved, in case reverse communication
c is being used.
save
c
c ... If re-entering after preconditioning, assume everything is set
c correctly and skip straight to solver.
if (status .gt. 0) goto 50
c
c ... Set control parameters as necessary.
if (mode .eq. 2) then
c
c Full manual mode, so skip straight to solver.
goto 40
else
c
c Need to save iparm(_JT_itmax_) for either black-box or hybrid mode.
user_itmax = iparm(_JT_itmax_)
if (mode .eq. 1) then
c
c Hybrid mode. Save other user-specified control parameters.
user_out = iparm(_JT_out_)
user_luout = iparm(_JT_luout_)
user_luerr = iparm(_JT_luerr_)
user_stop = iparm(_JT_stop_)
user_pre = iparm(_JT_pre_)
user_eps = rparm(_JT_eps_)
elseif (mode .ne. 0) then
if (iparm(_JT_out_) .ge. 2) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '***********************************'
write(iparm(_JT_luerr_),901) 'JTpack77: WARNING'
write(iparm(_JT_luerr_),907) ' Invalid value for mode: ', mode
write(iparm(_JT_luerr_),901) ' Continuing in black-box mode...'
write(iparm(_JT_luerr_),901) '***********************************'
mode = 0
endif
endif
endif
c
c ... Set defaults.
call JT_SetDefaults (iparm, rparm, status)
c
c Use || r || / || b || stopping test.
iparm(_JT_stop_) = 2
#ifndef unicos
c
c Use SSOR preconditioning.
iparm(_JT_pre_) = 2
#endif
c
c ... Restore user setting of iparm(_JT_itmax_).
iparm(_JT_itmax_) = user_itmax
c
c ... For hybrid mode, restore user settings of control parameters as
c appropriate.
if (mode .eq. 1) then
iparm(_JT_out_) = user_out
iparm(_JT_luout_) = user_luout
iparm(_JT_luerr_) = user_luerr
if (user_stop .ne. 0) iparm(_JT_stop_) = user_stop
iparm(_JT_pre_) = user_pre
if (user_eps .ne. zero) rparm(_JT_eps_) = user_eps
endif
c
40 continue
c
c ... Check to make sure iparm(_JT_itmax_) was set reasonably.
if (iparm(_JT_itmax_) .le. 0) then
if (iparm(_JT_out_) .ge. 2) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901)
& '**********************************************************'
write(iparm(_JT_luerr_),901) 'JTpack77: ERROR'
write(iparm(_JT_luerr_),907) ' Invalid value for itmax: ', iparm(_JT_itmax_)
write(iparm(_JT_luerr_),901) ' (itmax must be specified, and must'//
& ' be less than or equal to the upper'
write(iparm(_JT_luerr_),901) ' dimension of arrays cpu, rnormt,'//
& ' errt, rnorm, and err)'
write(iparm(_JT_luerr_),901)
& '**********************************************************'
endif
status = -1
return
endif
c
c ... Use black-box mode if no method is specified.
if (method.lt.0 .or. method.gt.3) then
if (iparm(_JT_out_) .ge. 2) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '***********************************'
write(iparm(_JT_luerr_),901) 'JTpack77: WARNING'
write(iparm(_JT_luerr_),907) ' Invalid value for method: ', method
write(iparm(_JT_luerr_),901) ' Continuing in black-box mode...'
write(iparm(_JT_luerr_),901) '***********************************'
method = 0
endif
endif
if (method .eq. 0) method = 2
c
50 continue
c
c ... Call appropriate solver routine.
if (method .eq. 1) then
c
c Standard conjugate gradient method.
iparm(_JT_resid_) = 2 ! replace r with b-Ax every MAX(10,SQRT(n))
! iterations
call JT_CG (b, a, ia, ja, iparm, rparm, x,
& cpu, rnormt, errt, rnorm, err, ap, iap, jap, status)
c
elseif (method .eq. 2) then
c
c Transpose-free QMR.
iparm(_JT_method_) = 1
call JT_BCGS (b, a, ia, ja, iparm, rparm, x,
& cpu, rnormt, errt, rnorm, err, ap, iap, jap, status)
c
c Protect against false convergence unless stopping test based on
c the true residual was used. Can't do this if stopping test 0
c was used either, but that test is not recommended for these
c methods anyway.
if (status.eq.0 .and. iparm(_JT_stop_).gt.0) then
c
c Negate iparm(_JT_stop_) to force JT_StoppingTest to calculate || b-Ax ||
c and the error estimate based on || b-Ax ||.
iparm(_JT_stop_) = -iparm(_JT_stop_)
c
c Calculate norm of true residual (b-Ax), and error estimate based
c on it.
bnorm = JT_VectorNorm(iparm(_JT_norm_), ia(_JT_nrows_), b, status)
rnormt_end = zero
call JT_StoppingTest (bnorm, rnorm(0), rnorm(iparm(_JT_iter_)), x, x, b,
& a, ia, ja, iparm, rparm, err_end, rnormt_end, errt_end,
& status)
if (errt_end .gt. rparm(_JT_eps_)) status = -4
c
c Restore iparm(_JT_stop_).
iparm(_JT_stop_) = -iparm(_JT_stop_)
endif
c
elseif (method .eq. 3) then
c
c Generalized minimal residual method.
call JT_GMRES (b, a, ia, ja, iparm, rparm, x,
& cpu, rnormt, errt, rnorm, err, ap, iap, jap, status)
endif
c
c ... Check return status.
if (status .lt. 0) then
if (status .eq. -1) then
c
c Invalid arguments to solver.
if (iparm(_JT_out_) .ge. 1) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '*****************************'
write(iparm(_JT_luerr_),901) 'JTpack77: ERROR'
write(iparm(_JT_luerr_),901) ' invalid arguments to solver'
write(iparm(_JT_luerr_),901) '*****************************'
endif
elseif (status .eq. -2) then
c
c Memory allocation failure.
if (iparm(_JT_out_) .ge. 1) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '*************************************'
write(iparm(_JT_luerr_),901) 'JTpack77: ERROR'
write(iparm(_JT_luerr_),901) ' memory allocation failure in solver'
write(iparm(_JT_luerr_),901) '*************************************'
endif
elseif (status .eq. -3) then
c
c Internal error.
if (iparm(_JT_out_) .ge. 1) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '**************************'
write(iparm(_JT_luerr_),901) 'JTpack77: ERROR'
write(iparm(_JT_luerr_),901) ' internal error in solver'
write(iparm(_JT_luerr_),901) '**************************'
endif
else
c
c Failure to converge.
if (method .eq. 1) then
c
c CG was used. Maybe system isn't really symmetric, so try
c TFQMR, whether failure was due to breakdown or itmax exceeded.
if (iparm(_JT_out_) .ge. 2) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '************************************'
write(iparm(_JT_luerr_),901) 'JTpack77: WARNING'
write(iparm(_JT_luerr_),901) ' CG failed'
write(iparm(_JT_luerr_),901) ' maybe system is really asymmetric?'
write(iparm(_JT_luerr_),901) ' trying TFQMR...'
write(iparm(_JT_luerr_),901) '************************************'
endif
method = 2
status = 0
goto 50
elseif (method .eq. 2) then
c
c TFQMR was used, so try GMRES, whether failure was due to
c breakdown, false convergence, or itmax exceeded.
if (iparm(_JT_out_) .ge. 2) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '*********************'
write(iparm(_JT_luerr_),901) 'JTpack77: WARNING'
write(iparm(_JT_luerr_),901) ' TFQMR failed'
write(iparm(_JT_luerr_),901) ' trying GMRES(20)...'
write(iparm(_JT_luerr_),901) '*********************'
endif
method = 3
status = 0
goto 50
else
c
c GMRES failed to converge.
if (iparm(_JT_nold_) .le. 40) then
c
c Try again with more vectors, but don't use more than 80.
if (iparm(_JT_out_) .ge. 2) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '*******************************'
write(iparm(_JT_luerr_),901) 'JTpack77: WARNING'
write(iparm(_JT_luerr_),903) ' GMRES(', iparm(_JT_nold_), ') failed'
write(iparm(_JT_luerr_),903) ' trying GMRES(', 2*iparm(_JT_nold_), ')...'
write(iparm(_JT_luerr_),901) '*******************************'
endif
iparm(_JT_nold_) = 2*iparm(_JT_nold_)
c
c Also increase number of steps if using Jacobi or SSOR
c preconditioning.
if (iparm(_JT_pre_).eq.1 .or. iparm(_JT_pre_).eq.2) then
iparm(_JT_steps_) = 2*iparm(_JT_steps_)
endif
c
c Use SSOR preconditioning if had been using none.
if (iparm(_JT_pre_) .eq. 0) iparm(_JT_pre_) = 2
c
status = 0
goto 50
else
c
c Give up.
if (iparm(_JT_out_) .ge. 1) then
write(iparm(_JT_luerr_),*)
write(iparm(_JT_luerr_),901) '**************************'
write(iparm(_JT_luerr_),901) 'JTpack77: ERROR'
write(iparm(_JT_luerr_),903) ' GMRES(', iparm(_JT_nold_), ') failed'
write(iparm(_JT_luerr_),901) ' giving up'
write(iparm(_JT_luerr_),901) '**************************'
endif
endif
endif
endif
endif
c
9999 continue
c
#include "formats.inc"
return
end