subroutine nitstb (n, xcur, fcur, fcnrm, step, eta, f, jacv, rpar, $ ipar, ijacv, irpre, iksmax, ifdord, nfe, njve, $ nrpre, nli, r, rtil, p, phat, v, t, rwork1, rwork2, $ rsnrm, dinpr, dnorm, itrmks) implicit none integer n, ipar(*), ijacv, irpre, iksmax, ifdord, nfe, njve, $ nrpre, nli, itrmks double precision xcur(n), fcur(n), fcnrm, step(n), eta, rpar(*), $ r(n), rtil(n), p(n), phat(n), v(n), t(n), rwork1(n), $ rwork2(n), rsnrm, dinpr, dnorm external f, jacv, dinpr, dnorm c ------------------------------------------------------------------------ c c This is nitstb v0.3, the BiCGSTAB routine for determining (trial) inexact c Newton steps. The original reference is H. van der Vorst, "Bi-CGSTAB: c A fast and smoothly converging variant of Bi-CG for the soluton of c nonsymmetric linear systems," SIAM J. Sci. Statist. Comput., 13 (1992), c pp. 631--644. c c ------------------------------------------------------------------------ c c Explanation: c c n = dimension of the problem. c c xcur = vector of length n, current approximate solution. c c fcur = vector of length n, value of f at xcur. c c fcnrm = norm of fcur. c c step = vector of length n, (trial) step. c c eta = relative residual reduction factor. c c f = name of user-supplied subroutine for evaluating the function c the zero of which is sought; this routine has the form c c subroutine f(n, xcur, fcur, rpar, ipar, itrmf) c c where xcur is the array containing the current x value, fcur c is f(xcur) on output, rpar and ipar are, respectively, real c and integer parameter/work arrays for use by the subroutine, c and itrmf is an integer termination flag. The meaning of c itrmf is as follows: c 0 => normal termination; desired function value calculated. c 1 => failure to produce f(xcur). c c jacv = name of user-supplied subroutine for evaluating J*v or c P(inverse)*v, where J is the Jacobian of f and P is a c right preconditioning operator. If neither analytic J*v c evaluations nor right preconditioning is used, this can c be a dummy subroutine; if right preconditioning is used but c not analytic J*v evaluations, this need only evaluate c P(inverse)*v. The form is c c subroutine jacv(n, xcur, fcur, ijob, v, z, rpar, ipar, itrmjv) c c where xcur and fcur are vectors of length n containing the c current x and f values, ijob is an integer flag indicating c which product is desired, v is a vector of length n to be c multiplied, z is a vector of length n containing the desired c product on output, rpar and ipar are, respectively, real c and integer parameter/work arrays for use by the subroutine, c and itrmjv is an integer termination c flag. The meaning of ijob is as follows: c 0 => z = J*v c 1 => z = P(inverse)*v c The meaning of itrmjv is as follows: c 0 => normal termination; desired product evaluated. c 1 => failure to produce J*v. c 2 => failure to produce P(inverse)*v. c This subroutine is called only from nitjv, and is always c called with v .ne. z. c c rpar = real parameter/work array passed to the f and jacv routines. c c ipar = integer parameter/work array passed to the f and jacv routines. c c ijacv = flag for determining method of J*v evaluation. c 0 => finite-difference evaluation (default). c 1 => analytic evaluation. c c irpre = flag for right preconditioning. c 0 => no right preconditioning c 1 => right preconditioning c c iksmax = maximum allowable number of BiCGSTAB iterations. c c ifdord = order of the finite-difference formula used in BiCGSTAB c when J*v products are evaluated using finite-differences. c When ijacv = 0 on input to nitsol, ifdord is set to 1, 2, or c 4 in nitsol; otherwise, it is irrelevant. When ijacv = 0 on c input to this subroutine, ifdord determines the order of the c finite-difference formula used at each BiCGSTAB iteration c (default 1). In this case, ijacv is set to -1 below to c signal to nitjv that the order of the finite-difference c formula is to be determined by ifdord. The original value c ijacv = 0 is restored on return. c c nfe = number of function evaluations. c c njve = number of J*v evaluations. c c nrpre = number of P(inverse)*v evaluations. c c nli = number of linear iterations. c c r = residual vector c c rtil = "r-tilde" vector used in BiCGSTAB c c p = vector used in BiCGSTAB c c phat = vector used in BiCGSTAB c c v = vector used in BiCGSTAB c c t = vector used in BiCGSTAB c c rwork1 = work vector, passed on to nitjv c c rwork2 = work vector, passed on to nitjv c c rsnrm = BiCGSTAB residual norm on return. c c dinpr = inner-product routine, either user-supplied or blas ddot. c c dnorm = norm routine, either user-supplied or blas dnrm2. c c itrmks = termination flag; values have the following meanings: c 0 => normal termination: acceptable step found. c 1 => J*v failure in nitjv. c 2 => P(inverse)*v failure in nitjv. c 3 => acceptable step not found in iksmax BiCGSTAB iterations. c 4 => BiCGSTAB breakdown. c c Note: On return, nitsol terminates if itrmks is 1 or 2. c If itrmks is 3 or 4, nitsol may terminate or continue. c In this event, the step returned is a meaningful inexact c Newton step only if the residual norm has been reduced. c A decision on termination/continuation is made in nitdrv c according to whether there is sufficient residual norm c reduction, even though the desired inexact Newton condition c may not hold. c c ------------------------------------------------------------------------- c c Subroutines required by this and all called routines: c c user supplied: f, jacv c c nitsol routines: nitjv, nitfd c c blas routines: daxpy, dcopy, dscal c c lapack routine: dlamch c c user supplied or blas: dinpr, dnorm c c explanation: In nitsol, dinpr and dnorm are set to either the blas c ddot and dnrm2 routines or the user-supplied usrnpr and usrnrm c routines. c c This subroutine called by: nitdrv c c Subroutines called by this subroutine: daxpy, dcopy, dscal, dinpr, dlamch, c dnorm, nitjv c c Common block: c include 'nitprint.h' c c If diagnostic information is desired, include this common block in the c main program and set iplvl and ipunit according to the following: c c iplvl = 0 => no printout c = 1 => iteration numbers and F-norms c = 2 => ... + some stats, step norms, and linear model norms c = 3 => ... + some Krylov solver and backtrack information c = 4 => ... + more Krylov solver and backtrack information c c ipunit = printout unit number. c c ------------------------------------------------------------------------ c double precision abstol, alpha, beta, omega, $ rho, rhomns, tau, temp integer i, istb, itask, itrmjv double precision dlamch external dlamch c double precision sfmin data sfmin /0.0d0/ c ------------------------------------------------------------------------ c c Initialize sfmin only on first entry. c if ( sfmin .eq. 0.0d0 ) sfmin = dlamch( 's' ) c ------------------------------------------------------------------------ c If finite-differences are used to evaluate J*v products (ijacv= 0), then c ijacv is set to -1 within this subroutine to signal to nitjv that the c order of the finite-difference formula is to be determined by ifdord. c The original value ijacv= 0 is restored on return. c ------------------------------------------------------------------------ if (ijacv .eq. 0) ijacv = -1 c ------------------------------------------------------------------------ c Set the stopping tolerance, initialize the step, etc. c ------------------------------------------------------------------------ rsnrm = fcnrm abstol = eta*rsnrm do 10 i = 1, n step(i) = 0.d0 10 continue istb = 0 c ------------------------------------------------------------------------ c For printing: if (iplvl .ge. 3) then write(ipunit,*) write(ipunit,800) eta endif 800 format('nitstb: eta =', 1pd10.3) if (iplvl .ge. 4) then write(ipunit,810) write(ipunit,*) write(ipunit,820) istb, rsnrm endif 810 format('nitstb: BiCGSTAB iteration no. (parts a and b)', $ ' linear residual norm, ') 820 format(5x,i4,5x,1pd10.3) c ------------------------------------------------------------------------ c c ------------------------------------------------------------------------ c Set up r and rtil. c ------------------------------------------------------------------------ call dcopy(n,fcur,1,r,1) temp = -1.d0 call dscal(n,temp,r,1) call dcopy(n,r,1,rtil,1) c ------------------------------------------------------------------------ c Top of the iteration loop. c ------------------------------------------------------------------------ 100 continue istb = istb + 1 nli = nli + 1 c ------------------------------------------------------------------------ c Perform the first "half-iteration". c ------------------------------------------------------------------------ rho = dinpr(n,rtil,1,r,1) if (istb .eq. 1) then call dcopy(n,r,1,p,1) else if ( abs(rhomns) .lt. sfmin*abs(rho) ) then itrmks = 4 goto 900 else beta = (rho/rhomns)*(alpha/omega) call daxpy(n,-omega,v,1,p,1) call dscal(n,beta,p,1) call daxpy(n,1.d0,r,1,p,1) endif endif if (irpre .eq. 0) then call dcopy(n,p,1,phat,1) else itask = 2 call nitjv(n, xcur, fcur, f, jacv, rpar, ipar, ijacv, $ ifdord, itask, nfe, njve, nrpre, p, phat, $ rwork1, rwork2, dnorm, itrmjv) if (itrmjv .gt. 0) then itrmks = 2 go to 900 endif endif itask = 0 call nitjv(n, xcur, fcur, f, jacv, rpar, ipar, ijacv, $ ifdord, itask, nfe, njve, nrpre, phat, v, $ rwork1, rwork2, dnorm, itrmjv) if (itrmjv .gt. 0) then itrmks = 1 go to 900 endif tau = dinpr(n,rtil,1,v,1) if ( abs(tau) .lt. sfmin*abs(rho) ) then itrmks = 4 goto 900 else alpha = rho/tau endif call daxpy(n,-alpha,v,1,r,1) call daxpy(n,alpha,phat,1,step,1) rsnrm = dnorm(n,r,1) c ------------------------------------------------------------------------ c For printing: if (iplvl .ge. 4) then write(ipunit,830) istb, rsnrm 830 format(5x,i4,'.a',3x,1pd10.3) endif c ------------------------------------------------------------------------ c c ------------------------------------------------------------------------ c Test for termination. c ------------------------------------------------------------------------ if (rsnrm .le. abstol) then itrmks = 0 go to 900 endif c ------------------------------------------------------------------------ c Perform the second "half-iteration". c ------------------------------------------------------------------------ if (irpre .eq. 0) then call dcopy(n,r,1,phat,1) else itask = 2 call nitjv(n, xcur, fcur, f, jacv, rpar, ipar, ijacv, $ ifdord, itask, nfe, njve, nrpre, r, phat, $ rwork1, rwork2, dnorm, itrmjv) if (itrmjv .gt. 0) then itrmks = 2 go to 900 endif endif itask = 0 call nitjv(n, xcur, fcur, f, jacv, rpar, ipar, ijacv, $ ifdord, itask, nfe, njve, nrpre, phat, t, $ rwork1, rwork2, dnorm, itrmjv) if (itrmjv .gt. 0) then itrmks = 1 go to 900 endif tau = dnorm(n,t,1) tau = tau*tau temp = dinpr(n,t,1,r,1) if ( tau .le. sfmin*abs(temp) ) then itrmks = 4 goto 900 else omega = temp/tau endif if ( abs(omega) .lt. sfmin*abs(alpha) ) then itrmks = 4 go to 900 endif call daxpy(n,-omega,t,1,r,1) call daxpy(n,omega,phat,1,step,1) rsnrm = dnorm(n,r,1) c ------------------------------------------------------------------------ c For printing: if (iplvl .ge. 4) then write(ipunit,840) istb, rsnrm 840 format(5x,i4,'.b',3x,1pd10.3) endif c ------------------------------------------------------------------------ c c ------------------------------------------------------------------------ c Test for termination. c ------------------------------------------------------------------------ if (rsnrm .le. abstol) then itrmks = 0 go to 900 endif if (istb .ge. iksmax) then itrmks = 3 go to 900 endif c ------------------------------------------------------------------------ c If continuing, update and return to the top of the iteration loop. c ------------------------------------------------------------------------ rhomns = rho go to 100 c ------------------------------------------------------------------------ c Bottom of the iteration loop. c ------------------------------------------------------------------------ c c ------------------------------------------------------------------------ c All returns made here. c ------------------------------------------------------------------------ 900 continue c ------------------------------------------------------------------------ c For printing: if (iplvl .ge. 3) then write(ipunit,*) if (itrmks .ne. 1 .and. itrmks .ne. 2) then write(ipunit,850) itrmks, rsnrm 850 format('nitstb: itrmks =', i2, ' final lin. res. norm =', $ 1pd10.3) else write(ipunit,860) itrmks 860 format('nitstb: itrmks:', i4) endif endif c ------------------------------------------------------------------------ c c ------------------------------------------------------------------------ c If ijacv = -1, then restore it to the original value ijacv = 0. c ------------------------------------------------------------------------ if (ijacv .eq. -1) ijacv = 0 return end