minresModule.f90

Go to the documentation of this file.

 
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File minresModule.f90
!
 
module minresmodule
 
  use  minresdatamodule, only : dp
 
  implicit none
  public   :: minres
 
contains
 
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  subroutine minres( n, Aprod, Msolve, b, shift, checkA, precon, &
                     x, itnlim, nout, rtol,                      &
                     istop, itn, Anorm, Acond, rnorm, Arnorm, ynorm )
 
    integer,  intent(in)    :: n, itnlim, nout
    logical,  intent(in)    :: checka, precon
    real(dp), intent(in)    :: b(n)
    real(dp), intent(in)    :: shift, rtol
    real(dp), intent(out)   :: x(n)
    integer,  intent(out)   :: istop, itn
    real(dp), intent(out)   :: anorm, acond, rnorm, arnorm, ynorm
 
    interface
       subroutine aprod (n,x,y)                   ! y := A*x
         use       minresdatamodule
         integer,  intent(in)    :: n
         real(dp), intent(in)    :: x(n)
         real(dp), intent(out)   :: y(n)
       end subroutine aprod
 
       subroutine msolve(n,x,y)                   ! Solve M*y = x
         use       minresdatamodule
         integer,  intent(in)    :: n
         real(dp), intent(in)    :: x(n)
         real(dp), intent(out)   :: y(n)
       end subroutine msolve
    end interface
 
!     Local arrays and variables
      real(dp)  :: r1(n), r2(n), v(n), w(n), w1(n), w2(n), y(n)
      real(dp)  :: alfa  , beta  , beta1 , cs    ,          &
                   dbar  , delta , denom , diag  ,          &
                   eps   , epsa  , epsln , epsr  , epsx  ,  &
                   gamma , gbar  , gmax  , gmin  ,          &
                   oldb  , oldeps, qrnorm, phi   , phibar,  &
                   rhs1  , rhs2  , rnorml, rootl ,          &
                   arnorml,        relarnorml,              &
                   s     , sn    , t     , tnorm2, ynorm2, z
      integer   :: i
      logical   :: debug, prnt
 
    ! Local constants
      real(dp),         parameter :: zero =  0.0,  one = 1.0
      real(dp),         parameter :: ten  = 10.0
      character(len=*), parameter :: enter = ' Enter MINRES.  '
      character(len=*), parameter :: exitt = ' Exit  MINRES.  '
      character(len=*), parameter :: msg(-1:8) =                  &
        (/ 'beta2 = 0.  If M = I, b and x are eigenvectors of A', & ! -1
           'beta1 = 0.  The exact solution is  x = 0           ', & !  0
           'Requested accuracy achieved, as determined by rtol ', & !  1
           'Reasonable accuracy achieved, given eps            ', & !  2
           'x has converged to an eigenvector                  ', & !  3
           'Acond has exceeded 0.1/eps                         ', & !  4
           'The iteration limit was reached                    ', & !  5
           'Aprod  does not define a symmetric matrix          ', & !  6
           'Msolve does not define a symmetric matrix          ', & !  7
           'Msolve does not define a pos-def preconditioner    ' /) !  8
    !-------------------------------------------------------------------
 
    !obsolete intrinsic       :: abs, dot_product, epsilon, min, max, sqrt
 
    ! Print heading and initialize.
 
    debug = .false.
    eps   = epsilon(eps)
    if (nout > 0) then
       write(nout, 1000) enter, n, checka, precon, itnlim, rtol, shift
    end if
    istop  = 0
    itn    = 0
    anorm  = zero
    acond  = zero
    rnorm  = zero
    ynorm  = zero
    x(1:n) = zero
    arnorml = zero
    gmin = zero
    gmax = zero
 
    !-------------------------------------------------------------------
    ! Set up y and v for the first Lanczos vector v1.
    ! y = beta1 P' v1, where P = C**(-1).
    ! v is really P' v1.
    !-------------------------------------------------------------------
    y      = b
    r1     = b
    if ( precon ) call msolve( n, b, y )
    beta1  = dot_product(b,y)
 
    if (beta1 < zero) then     ! M must be indefinite.
       istop = 8
       go to 900
    end if
 
    if (beta1 == zero) then    ! b = 0 exactly.  Stop with x = 0.
       istop = 0
       go to 900
    end if
 
    beta1  = sqrt( beta1 )     ! Normalize y to get v1 later.
 
    !-------------------------------------------------------------------
    ! See if Msolve is symmetric.
    !-------------------------------------------------------------------
    if (checka  .and.  precon) then
       call msolve( n, y, r2 )
       s      = dot_product(y ,y )
       t      = dot_product(r1,r2)
       z      = abs(s - t)
       epsa   = (s + eps) * eps**0.33333
       if (z > epsa) then
          istop = 7
          go to 900
       end if
    end if
 
    !-------------------------------------------------------------------
    ! See if Aprod  is symmetric.  Initialize Arnorm.
    !-------------------------------------------------------------------
    if (checka) then
       call aprod ( n, y, w )
       call aprod ( n, w, r2 )
       s      = dot_product(w,w )
       t      = dot_product(y,r2)
       z      = abs(s - t)
       epsa   = (s + eps) * eps**0.33333
       if (z > epsa) then
          istop = 6
          go to 900
       end if
       arnorml = sqrt(s);
    else
       call aprod ( n, y, w )
       arnorml = sqrt( dot_product(w,w) )
    end if
 
    !-------------------------------------------------------------------
    ! Initialize other quantities.
    !-------------------------------------------------------------------
    oldb   = zero
    beta   = beta1
    dbar   = zero
    epsln  = zero
    qrnorm = beta1
    phibar = beta1
    rhs1   = beta1
    rhs2   = zero
    tnorm2 = zero
    ynorm2 = zero
    cs     = - one
    sn     = zero
    w(1:n) = zero
    w2(1:n)= zero
    r2(1:n)= r1
 
    if (debug) then
       write(*,*) ' '
       write(*,*) 'b    ', b
       write(*,*) 'beta ', beta
       write(*,*) ' '
    end if
 
    !===================================================================
    ! Main iteration loop.
    !===================================================================
    do
       itn = itn + 1               ! k = itn = 1 first time through
 
       !----------------------------------------------------------------
       ! Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
       ! The general iteration is similar to the case k = 1 with v0 = 0:
       !
       !   p1      = Operator * v1  -  beta1 * v0,
       !   alpha1  = v1'p1,
       !   q2      = p2  -  alpha1 * v1,
       !   beta2^2 = q2'q2,
       !   v2      = (1/beta2) q2.
       !
       ! Again, y = betak P vk,  where  P = C**(-1).
       ! .... more description needed.
       !----------------------------------------------------------------
       s      = one / beta            ! Normalize previous vector (in y).
       v      = s*y(1:n)              ! v = vk if P = I
 
       call aprod ( n, v, y )
       y      = y - shift*v           ! call daxpy ( n, (- shift), v, 1, y, 1 )
       if (itn >= 2) then
          y   = y - (beta/oldb)*r1    ! call daxpy ( n, (- beta/oldb), r1, 1, y, 1 )
       end if
 
       alfa   = dot_product(v,y)      ! alphak
       y      = y - (alfa/beta)*r2    ! call daxpy ( n, (- alfa/beta), r2, 1, y, 1 )
       r1     = r2
       r2     = y
       if ( precon ) call msolve( n, r2, y )
 
       oldb   = beta                  ! oldb = betak
       beta   = dot_product(r2,y)     ! beta = betak+1^2
       if (beta < zero) then
          istop = 6
          go to 900
       end if
 
       beta   = sqrt( beta )          ! beta = betak+1
       tnorm2 = tnorm2 + alfa**2 + oldb**2 + beta**2
 
       if (itn == 1) then                   ! Initialize a few things.
          if (beta/beta1 <= ten*eps) then   ! beta2 = 0 or ~ 0.
             istop = -1                     ! Terminate later.
          end if
         !tnorm2 = alfa**2
          gmax   = abs( alfa )              ! alpha1
          gmin   = gmax                     ! alpha1
       end if
 
       ! Apply previous rotation Qk-1 to get
       !   [deltak epslnk+1] = [cs  sn][dbark    0   ]
       !   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
 
       oldeps = epsln
       delta  = cs * dbar  +  sn * alfa ! delta1 = 0         deltak
       gbar   = sn * dbar  -  cs * alfa ! gbar 1 = alfa1     gbar k
       epsln  =               sn * beta ! epsln2 = 0         epslnk+1
       dbar   =            -  cs * beta ! dbar 2 = beta2     dbar k+1
 
       ! Compute the next plane rotation Qk
 
       gamma  = sqrt( gbar**2 + beta**2 )   ! gammak
       cs     = gbar / gamma                ! ck
       sn     = beta / gamma                ! sk
       phi    = cs * phibar                 ! phik
       phibar = sn * phibar                 ! phibark+1
 
       if (debug) then
          write(*,*) ' '
          write(*,*) 'v    ', v
          write(*,*) 'alfa ', alfa
          write(*,*) 'beta ', beta
          write(*,*) 'gamma', gamma
          write(*,*) 'delta', delta
          write(*,*) 'gbar ', gbar
          write(*,*) 'epsln', epsln
          write(*,*) 'dbar ', dbar
          write(*,*) 'phi  ', phi
          write(*,*) 'phiba', phibar
          write(*,*) ' '
       end if
 
       ! Update  x.
 
       denom = one/gamma
 
       do i = 1, n
          w1(i) = w2(i)
          w2(i) = w(i)
          w(i)  = ( v(i) - oldeps*w1(i) - delta*w2(i) ) * denom
          x(i)  =   x(i) +   phi * w(i)
       end do
 
       ! Go round again.
 
       gmax   = max( gmax, gamma )
       gmin   = min( gmin, gamma )
       z      = rhs1 / gamma
       ynorm2 = z**2  +  ynorm2
       rhs1   = rhs2  -  delta * z
       rhs2   =       -  epsln * z
 
       ! Estimate various norms and test for convergence.
 
       anorm  = sqrt( tnorm2 )
       ynorm  = sqrt( ynorm2 )
       epsa   = anorm * eps
       epsx   = anorm * ynorm * eps
       epsr   = anorm * ynorm * rtol
       diag   = gbar
       if (diag == zero) diag = epsa
 
       qrnorm = phibar
       rnorml = rnorm
       rnorm  = qrnorm
       rootl       = sqrt( gbar**2 +dbar**2  )  ! norm([gbar; dbar]);
       arnorml     = rnorml*rootl               ! ||A r_{k-1} ||
       relarnorml  = rootl  /  anorm;           ! ||Ar|| / (||A|| ||r||)     
       !relArnorml = Arnorml / Anorm;           ! ||Ar|| / ||A|| 
 
       ! Estimate  cond(A).
       ! In this version we look at the diagonals of  R  in the
       ! factorization of the lower Hessenberg matrix,  Q * H = R,
       ! where H is the tridiagonal matrix from Lanczos with one
       ! extra row, beta(k+1) e_k^T.
 
       acond  = gmax / gmin
 
       ! See if any of the stopping criteria are satisfied.
       ! In rare cases, istop is already -1 from above (Abar = const*I).
 
       if (istop == 0) then
          if (itn    >= itnlim    ) istop = 5
          if (acond  >= 0.1d+0/eps) istop = 4
          if (epsx   >= beta1     ) istop = 3
          ! original
          !if (qrnorm <= epsx  .or.  relArnorml <= epsx) istop = 2
          !if (qrnorm <= epsr  .or.  relArnorml <= epsr) istop = 1
          ! C. Kleinwort, DESY, 131002
          if (qrnorm <= epsx  .or.  relarnorml <= eps) istop = 2
          if (qrnorm <= epsr  .or.  relarnorml <= rtol) istop = 1
       end if
 
 
       ! See if it is time to print something.
 
       if (nout > 0) then
          prnt   = .false.
          if (n      <= 40         ) prnt = .true.
          if (itn    <= 10         ) prnt = .true.
          if (itn    >= itnlim - 10) prnt = .true.
          if (mod(itn,10)  ==     0) prnt = .true.
          if (qrnorm <=  ten * epsx) prnt = .true.
          if (qrnorm <=  ten * epsr) prnt = .true.
          if (relarnorml<= ten*epsx) prnt = .true.
          if (relarnorml<= ten*epsr) prnt = .true.
          if (acond  >= 1.0d-2/eps ) prnt = .true.
          if (istop  /=  0         ) prnt = .true.
 
          if ( prnt ) then
             if (    itn     == 1) write(nout, 1200)
             write(nout, 1300) itn, x(1), qrnorm, anorm, acond
             if (mod(itn,10) == 0) write(nout, 1500)
          end if
       end if
       if (istop /= 0) exit
 
    end do
    !===================================================================
    ! End of iteration loop.
    !===================================================================
 
    ! Display final status.
 
900 arnorm = arnorml
    if (nout  > 0) then
       write(nout, 2000) exitt, istop, itn,   &
                         exitt, anorm, acond, &
                         exitt, rnorm, ynorm, arnorm
       write(nout, 3000) exitt, msg(istop)
    end if
 
    return
 
 1000 format(// 1p,    a, 5x, 'Solution of symmetric   Ax = b'    &
              / ' n      =', i7, 5x, 'checkA =', l4, 12x,         &
                 'precon =', l4                                   &
              / ' itnlim =', i7, 5x, 'rtol   =', e11.2, 5x,       &
                 'shift  =', e23.14)
 1200 format(// 5x, 'itn', 8x, 'x(1)', 10x,                       &
                'norm(r)', 3x, 'norm(A)', 3x, 'cond(A)')
 1300 format(1p, i8, e19.10, 3e10.2)
 1500 format(1x)
 2000 format(/ 1p, a, 5x, 'istop =', i3,   14x, 'itn   =', i8     &
             /     a, 5x, 'Anorm =', e12.4, 5x, 'Acond =', e12.4  &
             /     a, 5x, 'rnorm =', e12.4, 5x, 'ynorm =', e12.4, 5x, 'Arnorml =', e12.4)
 3000 format(      a, 5x, a )
 
  end subroutine minres
 
end module minresmodule