gmres.f90

Go to the documentation of this file.

 
      MODULE gmres
      USE v3_utilities, ONLY: assert
      USE stel_kinds
      USE stel_constants, ONLY: zero, one
      USE descriptor_mod, ONLY: siesta_comm, iam, nprocs
      USE shared_data, etak=>etak_tol
      USE shared_functions, ONLY: funct_island, linesearch
      USE siesta_state, ONLY: update_state, update_state
      USE nscalingtools, ONLY: parfunctisl, mpi_err, startglobrow,      &
                               endglobrow, parsolver
      USE timer_mod
      USE mpi_inc
 
      IMPLICIT NONE
 
      CONTAINS
 
      SUBROUTINE gmres_fun
      USE hessian, ONLY: levmarq_param, dmin=>mindiag, apply_precond,   &
                         levmarq_param0, l_diagonal_only, gather_array, &
                         mupar
      USE siesta_namelist, ONLY: lresistive, ftol
      USE blocktridiagonalsolver_s, ONLY: refactorhessian
 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: j, iter, istat
      REAL(dp), PARAMETER :: levscan(1) = (/16._dp/)           !!007C
!      REAL(dp), PARAMETER :: levscan(2) = (/4._dp, 16._dp/)   !!007B
      REAL(dp) :: wmhd, fmhd, skston, skstoff, eps, lm0, fsq_min,       &
                  xmax_lin, mupars, scale_fac
      REAL(dp), ALLOCATABLE :: xc1(:), xcmin(:), brhs(:)
      LOGICAL               :: bTerm
!-----------------------------------------------
 
!     SPH: ONLY works with 0 restart each time
!     EXTERNALLY, l_init_state = .TRUE.   
 
      xc = 0
      l_linearize = .false.
      l_applyprecon = .false.    !LET GMRES APPLY PRECONDITIONER ITSELF
      mupars = mupar
      
!
!     STORE INITIAL UNPRECONDITIONED RESIDUE (BRHS=-GC AT INIT PT)
!     SO DEVIATIONS FROM THEM CAN BE COMPUTED REPEATEDLY
!
      l_getfsq =.true.
      l_init_state=.true.
      fmhd=fsq_total1               !!initial fsq BEFORE application of resistive diffusion
 
      ALLOCATE (brhs(neqs), stat=istat)
      CALL assert(istat.eq.0,'ALLOCATION ERROR IN GMRES_FUN')
      l_getwmhd=.true.
 
 
      CALL second0(skston)
      CALL funct_island
      brhs = -gc 
      IF (parfunctisl) THEN
         CALL gather_array(brhs)
      END IF
      l_getwmhd=.false.
      CALL second0(skstoff)
      gmres_funct_island_time=gmres_funct_island_time+(skstoff-skston)
 
      wmhd = wtotal
      fsq_min = fsq_total1
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FAST CONVERGENCE - SPH 083016 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ALLOCATE(xc1(neqs), xcmin(neqs), stat=istat)
      CALL assert(istat.eq.0,'Allocation error in GMRES')
      xcmin = 0
 
      iter = 0
      lm0 = levmarq_param
      scale_fac = 1
 
920   CONTINUE
!
!     FAST CONVERGENCE GUESS
!     ITERATIVELY SOLVE (A* - eps*I) x = b, WHERE A* = A + eps*I 
!     eps = lev_marq_param*Dmin
!     THEN x = x0 + x1 + ..., 
!     x0 = P*b   where P* = (A*)^-1
!     xn = eps*P*(x[n-1])
!
      lfast: IF (.NOT.l_diagonal_only .AND. lcolscale) THEN
 
      IF (iam.EQ.0 .AND. lverbose) print 925, dmin
925   FORMAT(' DMIN: ',1pe12.3,/,1x,'FAST CONVERGENCE SUMMARY',         &
             /,' -------------',/,1x,'ITER',7x,'FSQ_NL',                &
             10x,'||X||',9x,'MAX|X|')
      CALL assert(all(xc.EQ.zero),' XC != 0')
      l_getfsq=.true.
      l_linearize = .false.
      CALL funct_island
      xc = -gc
      eps = abs(levmarq_param)*dmin
      IF (.not.l_applyprecon) THEN
         CALL apply_precond(xc)
      END IF
      xc1 = xc
      DO j = 1, 10
         IF (j .GT. 1) THEN
            IF (eps .EQ. zero) GOTO 927
            gc = xc1
            CALL apply_precond(gc)
            xc1 = eps*gc            !nth partial sum: store for next j iteration
            xc = xc + xc1
         END IF
         CALL funct_island
 
         IF (fsq_total1 .GT. fsq_min) THEN
            bterm = (j .GT. 1) 
         ELSE
            xcmin = xc
            bterm = (fsq_total1 .GE. 0.95_dp*fsq_min)
            fsq_min = fsq_total1
         END IF
         IF (iam.EQ.0 .AND. lverbose) print 912, j, fsq_total1,         &
            sqrt(sum(xc*xc)), maxval(abs(xc))
         IF (fsq_min.LE.ftol .OR. bterm .OR. eps.EQ.zero) EXIT
      END DO
      
!      ELSE IF (l_backslv) THEN
!         IF (.NOT.l_ApplyPrecon) CALL Apply_Precond(gc)
      ENDIF lfast
 
 927  CONTINUE
 
      lgmres: IF (fsq_min.GT.ftol .AND. fsq_min.GT.fmhd/10) THEN
 
!!!START GMRES CONVERGENCE
!
!     THIS IS STANDARD GMRES CALL
!
      CALL second0(skston)
      xc = xcmin
      CALL gmres_wrap (xc, brhs)
      CALL second0(skstoff)
      gmres_wrap_time=gmres_wrap_time+(skstoff-skston)
 
      IF (fsq_total1.LT.fsq_min .OR. all(xcmin.EQ.zero)) THEN
         xcmin = xc
         fsq_min = fsq_total1
         IF (lm0 .GT. zero) scale_fac = levmarq_param/lm0
      END IF
!!!END GMRES CONVERGENCE
!FOR NOW, THIS IS ONLY IMPLEMENTED FOR PARSOLVER=TRUE: MUST IMPLEMENT
!RefactorHessian FOR LSCALAPACK=T
      IF (.false.) THEN
!      IF (.NOT.l_Diagonal_Only .AND. PARSOLVER) THEN
         IF (fsq_min.GE.0.95_dp*fmhd .AND. iter.LT.SIZE(levscan)) THEN
            iter = iter+1
            levmarq_param = lm0*levscan(iter)
            IF (levmarq_param .LT. levm_ped) levmarq_param = 100*levm_ped
            IF (levmarq_param .GE. levmarq_param0) levmarq_param = levmarq_param0/10._dp**iter
            IF (iam.EQ.0 .AND. lverbose) print 930,' Refactoring Hessian: LM = ',levmarq_param
            CALL refactorhessian(levmarq_param)
            GOTO 920
         END IF
      END IF
      END IF lgmres
 
 912  FORMAT(i5,3(3x,1pe12.3))
 930  FORMAT (/,a,1pe12.3)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
      DEALLOCATE(xc1)
      l_getfsq = .true.
 
 1100 CONTINUE
 
      xmax_lin = maxval(abs(xcmin))
      IF (iam .EQ. 0) THEN
!         PRINT 1020,' ||X||-GMRES = ', SQRT(SUM(xcmin*xcmin)), ' MAX(|X|) = ',xmax_lin
         WRITE (unit_out, 1020) ' ||X||-GMRES = ', sqrt(sum(xcmin*xcmin)), ' MAX(|X|) = ',xmax_lin
      END IF
 1020 FORMAT(2(a,1p,e11.3))
 
      fsq_lin = fsq_total
      l_linearize = .false.
 
!UPDATE EVOLUTION FIRST
      IF (fsq_min .GT. ftol) THEN
 
!     TURN OFF NEXT 2 LINES FOR TESTING f1=levm IN EVOLUTION
         levm_scale = levm_scale*scale_fac  
         levm_scale = min(100._dp, levm_scale)      
!
!     SIMPLE LINE SEARCH SCAN OF |F| FOR VARIOUS |X| 
         mupar = 0
         CALL linesearch(xcmin, fsq_min)
 
      END IF
 
      xc = xcmin
!
!     RECOMPUTE PERTURBED FIELDS, PRESSURE FOR MINIMUM STATE
!     THEN UPDATE NONLINEAR STATE BY ADDING XC (PERTURBATION)
!
 
      l_printoriginforces=.true.
      l_init_state=.true.
      l_getwmhd= .true.
      CALL second0(skston)
      CALL funct_island
!      CALL update_state(.FALSE., zero, zero)                           !Get this in evolution, line 455
      CALL second0(skstoff)
      gmres_funct_island_time=gmres_funct_island_time+(skstoff-skston)
      l_printoriginforces=.false.
 
      mupar = mupars
 
      DEALLOCATE (xcmin, brhs)
      
      END SUBROUTINE gmres_fun
 
      SUBROUTINE init_gmres0 (icntl, cntl, etak, m, n)
      INTEGER, INTENT(OUT) :: icntl(9), m, n
      REAL(dp),INTENT(OUT) :: cntl(5)
      REAL(dp), INTENT(IN) :: etak
      REAL(dp)             :: skston, skstoff
      INTEGER, PARAMETER   :: noPrec=0, leftprec=1, rightprec=2, dbleprec=3
 
      CALL second0(skston)
      CALL init_dgmres(icntl ,cntl)
      CALL second0(skstoff)
      gmres_init_dgmres_time=gmres_init_dgmres_time+(skstoff-skston)
 
!*************************
!*  Tune some parameters
!*************************
! Tolerance
       cntl(1) = etak
!       cntl(1) = 1.E-5_dp
! Write errors to fort.21
!      icntl(1) = 21   !21
! Write warnings to fort.21
      icntl(2) = 21
! Save the convergence history in file fort.20
      IF (nprecon .gt. 1) icntl(3) = 20
! Preconditioning INPUT flags (note: different than revcom flags in driver)
!      icntl(4) = leftPrec
      icntl(4) = rightprec 
!      icntl(4) = dblePrec
!! ICGS orthogonalization
      icntl(5) = 3
! Initial guess
      icntl(6) = 0
!      icntl(6) = 1
! Maximum number of iterations at each step (~ns/5)
      icntl(7) = ngmres_steps
 
      icntl(8) = 1             !Default
      icntl(9) = 1             !Steps for peek at progress during rev com loop
!*********************************
! Choose the restart parameter
!*********************************
!      write(*,*) 'Restart  <', ldstrt
!      read(*,*) m
!
!     m <= n
!
      n   = neqs
      m = 200
 
      END SUBROUTINE init_gmres0
 
      SUBROUTINE gmres_wrap (x0, b)
      USE siesta_namelist, ONLY: ftol
      USE nscalingtools, ONLY: pargmres, rcounts, disp
      USE hessian, ONLY: apply_precond
      USE gmres_lib, ONLY: gmres_info, gmres_par, gmres_ser
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), INTENT(IN)    :: b(:)
      REAL(dp), INTENT(INOUT) :: x0(:)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      TYPE(gmres_info)   :: gi
      INTEGER :: n
!-----------------------------------------------
      l_linearize = .true.
      l_getfsq = .false.         !So call to funct in matvec does not include gc0
 
!******************************************************
!*  Initialize GMRES control parameters and Load GMRES_INFO structure
!*******************************************************
      CALL init_gmres0(gi%icntl, gi%cntl, etak, gi%m, n)
      gi%ftol=ftol
      gi%lverbose = lverbose
      gi%icntl(6) = 1
 
!
!     EASY-TO-USE WRAPPER FOR GMRES DRIVER CALL
!
!     X0: on input, initial guess if icntl(6) == 1
!         on output, solution of Ax = b
!         NOTE: it is not overwritten UNTIL the end of this routine
      gi%iam=iam; gi%nprocs=nprocs; gi%ngmres_type=ngmres_type
      IF (parsolver .AND. pargmres) THEN
         gi%endglobrow=endglobrow; gi%startglobrow=startglobrow
         gi%mblk_size=mblk_size; gi%rcounts=>rcounts; gi%disp=>disp
 
         gi%my_comm = siesta_comm
         gi%my_comm_world = siesta_comm
!         gi%lactive = lactive
         CALL gmres_par (n, gi, matvec_par, apply_precond, getnlforce,  &
                         x0, b)
      ELSE
         CALL gmres_ser (n, gi, matvec, apply_precond, getnlforce,      &
                         x0, b)
      END IF
 
      fsq_total1 = gi%ftol
      CALL assert(all(xc.EQ.x0),'XC != X0 IN GMRES_WRAP')
 
      END SUBROUTINE gmres_wrap
 
      SUBROUTINE matvec_par (ploc, Ap, nloc)
      USE island_params, ONLY: ns=>ns_i
      USE hessian, ONLY: eps_factor, mupar, gather_array
      USE blocktridiagonalsolver_s, ONLY: parmatvec, padsides 
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN)                    :: nloc
      REAL(dp), INTENT(IN), DIMENSION(nloc)  :: ploc
      REAL(dp), INTENT(OUT), DIMENSION(nloc) :: Ap
!-----------------------------------------------
      LOGICAL, PARAMETER    :: LPARMATVEC=.false.
      INTEGER               :: myrowstart, myrowend, istat
      REAL(dp), ALLOCATABLE :: p(:)
      REAL(dp)              :: delta, skston, skstoff
!-----------------------------------------------
!SPH NOTE 032713: EVENTUALLY REMOVE CALL TO ParMatVec IF THE FUNCT_ISLAND IS FASTER
!
!     NOTE: DO NOT CALL ParMatVec WHEN nprecon_type==PREDIAG
!           SINCE IN HESSIAN, WE SET ALL THE OFF-DIAGONAL COMPONENTS OF 
!           THE A-MATRIX  (USED IN ParMatVec) TO ZERO
!
!           AT PRESENT, ParMatVec WILL NOT WORK IF mupar != 0
!
      CALL assert(.NOT.l_init_state,'l_init_state = T in matvec_par')
      istat = (endglobrow-startglobrow+1)*mblk_size
      CALL assert(istat.EQ.nloc, 'nloc wrong in matvec_par')
 
      myrowstart=(startglobrow-1)*mblk_size+1
      myrowend=myrowstart+nloc-1
 
      CALL second0(skston)
      ALLOCATE (p(ns*mblk_size), stat=istat)
      CALL assert(istat.eq.0,'Allocation error in matvec_par')
      p = 0
 
      p(myrowstart:myrowend) = ploc
      CALL padsides(p, mblk_size, 1, 1)
 
      CALL second0(skstoff)
      gmres_wrap_allgather_time=gmres_wrap_allgather_time+(skstoff-skston)
 
      CALL second0(skston)
 
      IF (.NOT.lparmatvec .OR. nprecon_type.EQ.prediag                  &
          .OR. mupar.NE.zero) THEN
         delta = sqrt(epsilon(delta))*eps_factor
         xc = delta*p
 
         l_linearize=.true.
         CALL funct_island
         ap = gc(myrowstart:myrowend)/delta
 
      ELSE 
         CALL parmatvec(p,ap,nloc)
      END IF
 
      DEALLOCATE (p)
 
      WHERE (abs(ap) .LT. 1.e-10_dp) ap = 0
 
      CALL second0(skstoff)
 
      paryax_time=paryax_time+(skstoff-skston)
 
      END SUBROUTINE matvec_par
 
      SUBROUTINE matvec (p, Ap, ndim)
      USE stel_kinds
      USE hessian, ONLY: eps_factor, gather_array
#if defined(MPI_OPT)
      USE mpi_inc
#endif
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN)   :: ndim
      REAL(dp), INTENT(IN)  :: p(ndim)
      REAL(dp), INTENT(OUT) :: Ap(ndim)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      REAL(dp), PARAMETER   :: zero=0
      REAL(dp)              :: delta
!-----------------------------------------------
      CALL assert(.not.l_init_state,'l_init_state = T in matvec')
 
      delta = sqrt(epsilon(delta))*eps_factor
 
!     Note: GMRES without Preconditioner has pnorm = 1 (checked!)
!      pnorm = SUM(p*p)
!      delta = delta/SQRT(pnorm)
!      delta = delta*(pnorm + SUM(p*xc0))/pnorm
 
!IF CALLED IN PAR MODE, FIRST GATHER THE xc's
      CALL assert(SIZE(gc).EQ.SIZE(ap),'gc and Ap wrong sizes')
      xc = delta*p
 
      l_getwmhd=.true.
      CALL funct_island
      l_getwmhd=.false.
 
      IF (l_linearize) THEN
         ap = gc/delta
      ELSE
         ap = (gc-gc0)/delta
      END IF
 
      WHERE (abs(ap) .LT. 1.e-10_dp) ap = 0
 
      END SUBROUTINE matvec
 
 
      SUBROUTINE get_etak(fk, fkm, ftol, etak)
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp), INTENT(in)    :: fk, fkm, ftol
      REAL(dp), INTENT(inout) :: etak
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      REAL(dp), PARAMETER :: gamma1=0.9_dp, alpha=1.5_dp, eta0 = 1.e-02_dp
      REAL(dp)            :: etakm
!-----------------------------------------------
!     ROUTINE PROVIDED L. CHACON (11/09/06)
      etakm = etak
 
!     Superlinear convergence
      etak = gamma1*(fk/fkm)**alpha
 
!     Safeguard avoid sharp decrease in etak
      etak = min(eta0, max(etak, gamma1*etakm**alpha))
 
!     Safeguard avoid "oversolving"
      etak = min(eta0, max(etak, gamma1*ftol/fk))
 
      END SUBROUTINE get_etak
 
 
      SUBROUTINE getnlforce(xcstate, fsq_nl, bnorm)
#if defined(MPI_OPT)
      USE blocktridiagonalsolver_s, ONLY: padsides 
#endif
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      REAL(dp),INTENT(IN)  :: xcstate(neqs), bnorm
      REAL(dp),INTENT(OUT) :: fsq_nl
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER              :: nloc, myrowstart, myrowend
      LOGICAL              :: l_linloc, l_Apploc, l_getloc
!-----------------------------------------------
      CALL assert(.not.l_init_state,'l_init_state=T in GetNLForce!')
!Store state 
      l_linloc = l_linearize
      l_apploc = l_applyprecon
      l_getloc = l_getfsq
      
      nloc=(endglobrow-startglobrow+1)*mblk_size
      myrowstart=(startglobrow-1)*mblk_size+1
      myrowend=myrowstart+nloc-1
 
 
!Set state variables
      xc(myrowstart:myrowend) = bnorm*xcstate(myrowstart:myrowend)       !undo internal gmres normalization
#if defined(MPI_OPT)
      IF (parsolver) CALL padsides(xc, mblk_size, 1, 1)
#endif
      l_init_state = .false.
      l_linearize = .false.
      l_getfsq = .true.
      l_applyprecon=.false.
 
      CALL funct_island
      fsq_nl = fsq_total1
 
!Restore state variables
      l_linearize = l_linloc
      l_applyprecon = l_apploc
      l_getfsq = l_getloc
 
      END SUBROUTINE getnlforce
 
      SUBROUTINE qmr_fun
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: ndim, nlen, nlim, ierr, info(4), j
      INTEGER :: revcom, colx, colb
      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: vecs
      REAL(dp) :: tol = 1.e-4_dp, gnorm
!-----------------------------------------------
!
!     STORE INITIAL POINT AND INITIAL RESIDUE (AT INIT PT)
!     SO DEVIATIONS FROM THEM CAN BE COMPUTED REPEATEDLY
!
      ALLOCATE(xc0(neqs), gc0(neqs), vecs(neqs,9), stat=j)
      CALL assert(j.eq.0,'Allocation error in qmr_fun')
 
      ndim = SIZE(vecs,1)
      nlen = ndim
      nlim = 100
 
!     SPH: ONLY works with 0 restart each time??!!!     
      l_linearize = .false.
      l_applyprecon = .true.
 
      xc = 0
      CALL funct_island
      xc0 = xc
      gc0 = gc
 
      l_linearize = .true.
!
!     INITIALIZE vecs
!
      gnorm = sum(gc*gc)
      gnorm = sqrt(gnorm)
      vecs(:ndim,2) = -gc(:ndim)/gnorm
      vecs(:ndim,3) =  gc(:ndim)/gnorm
    
      ierr = 100000
      info = 0
      info(1) = ierr
 
 10   CALL dutfx (ndim,nlen,nlim,vecs,tol,info)
      revcom = info(2)
      colx   = info(3)
      colb   = info(4)
      IF (revcom .eq. 1) THEN
         CALL matvec (vecs(1,colx), vecs(1,colb), ndim)
         GO TO 10
      END IF
 
      xc(1:ndim) = xc0(1:ndim) + gnorm*vecs(:,1)
 
!
!     RECOMPUTE delta current,pressure prior to update_state call
!
      CALL funct_island
      fsq_lin = fsq_total
 
      l_linearize = .false.
      CALL funct_island
 
      DEALLOCATE (xc0, gc0, vecs)
 
      END SUBROUTINE qmr_fun
 
      END MODULE gmres