minres.f90

Go to the documentation of this file.

! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16  Time: 11:06:43



SUBROUTINE minres( n, b, r1, r2, v, w, w1, w2, x, y,  &
aprod, msolve, checka, precon, shift, nout , itnlim, rtol,  &
istop, itn, anorm, acond, rnorm, ynorm )

    IMPLICIT NONE
    EXTERNAL aprod, msolve
    INTEGER :: n, nout, itnlim, istop, itn
    LOGICAL :: checka, precon
    DOUBLE PRECISION :: shift, rtol, anorm, acond, rnorm, ynorm,  &
    b(n), r1(n), r2(n), v(n), w(n), w1(n), w2(n), x(n), y(n)

    EXTERNAL ddot  , dnrm2
    DOUBLE PRECISION :: ddot  , dnrm2

    !     Local variables

    DOUBLE PRECISION :: alfa  , beta  , beta1 , cs    ,  &
    dbar  , delta , denom , diag  , eps   , epsa  , epsln , epsr  , epsx  ,  &
    agamma, gbar  , gmax  , gmin  , oldb  , oldeps, qrnorm, phi   , phibar,  &
    rhs1  , rhs2  , s     , sn    , t     , tnorm2, ynorm2, z
    INTEGER :: i
    LOGICAL :: debug, prnt

    DOUBLE PRECISION :: zero, one, two, ten
    parameter( zero = 0.0d+0,  one =  1.0d+0,  &
    two  = 2.0d+0,  ten = 10.0d+0 )

    CHARACTER (LEN=16) :: enter, exit
    CHARACTER (LEN=52) :: msg(-1:8)

    DATA               enter /' Enter MINRES.  '/, exit  /' Exit  MINRES.  '/

    DATA               msg  &
    / 'beta2 = 0.  If M = I, b and x are eigenvectors of A',  &
    'beta1 = 0.  The exact solution is  x = 0',  &
    'Requested accuracy achieved, as determined by rtol',  &
    'Reasonable accuracy achieved, given eps',  &
    'x has converged to an eigenvector', 'Acond has exceeded 0.1/eps',  &
    'The iteration limit was reached',  &
    'Aprod  does not define a symmetric matrix',  &
    'Msolve does not define a symmetric matrix',  &
    'Msolve does not define a pos-def preconditioner' /
    !     ------------------------------------------------------------------

    debug = .false.

    !     ------------------------------------------------------------------
    !     Compute eps, the machine precision.  The call to daxpy is
    !     intended to fool compilers that use extra-length registers.
    !     31 May 1999: Hardwire eps so the debugger can step thru easily.
    !     ------------------------------------------------------------------
    eps    = 2.22d-16    ! Set eps = zero here if you want it computed.

    eps  = 0.0d0                             !!!!!!!!!!!!!!!
    gmin = 0.0d0
    gmax = 0.0d0

    IF (eps <= zero) THEN
        eps    = two**(-12)
        DO
            eps    = eps / two
            x(1)   = eps
            y(1)   = one
            CALL daxpy( 1, one, x, 1, y, 1 )
            IF (y(1) <= one) exit
        END DO
        eps    = eps * two
    END IF
    !     WRITE(*,*) 'computed epsilon is ',eps   !!!!!!!!!!!!!!!

    !     ------------------------------------------------------------------
    !     Print heading and initialize.
    !     ------------------------------------------------------------------
    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
    CALL dload2( n, zero, x )

    !     ------------------------------------------------------------------
    !     Set up y and v for the first Lanczos vector v1.
    !     y  =  beta1 P' v1,  where  P = C**(-1).
    !     v is really P' v1.
    !     ------------------------------------------------------------------
    CALL dcopy( n, b, 1, y , 1 )         ! y  = b
    CALL dcopy( n, b, 1, r1, 1 )         ! r1 = b
    IF ( precon ) CALL msolve( n, b, y )
    beta1  = ddot( n, b, 1, y, 1 )

    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      = ddot( n, y, 1, y, 1 )
        t      = ddot( n,r1, 1,r2, 1 )
        z      = abs( s - t )
        epsa   = (s + eps) * eps**0.33333d+0
        IF (z > epsa) THEN
            istop = 7
            go to 900
        END IF
    END IF

    !     ------------------------------------------------------------------
    !     See if Aprod  is symmetric.
    !     ------------------------------------------------------------------
    IF (checka) THEN
        CALL aprod( n, y, w )
        CALL aprod( n, w, r2 )
        s      = ddot( n, w, 1, w, 1 )
        t      = ddot( n, y, 1,r2, 1 )
        z      = abs( s - t )
        epsa   = (s + eps) * eps**0.33333d+0
        IF (z > epsa) THEN
            istop = 6
            go to 900
        END IF
    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
    CALL dload2( n, zero, w  )        ! w  = 0
    CALL dload2( n, zero, w2 )        ! w2 = 0
    CALL dcopy( n, r1, 1, r2, 1 )    ! r2 = 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
        IF (istop /= 0) exit

        !-----------------------------------------------------------------
        ! 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).
        CALL dscal2( n, s, y, v )      ! v = vk if P = I

        CALL aprod( n, v, y )
        CALL daxpy( n, (- shift), v, 1, y, 1 )
        IF (itn >= 2) THEN
            CALL daxpy( n, (- beta/oldb), r1, 1, y, 1 )
        END IF

        alfa   = ddot( n, v, 1, y, 1 )     ! alphak

        CALL daxpy( n, (- alfa/beta), r2, 1, y, 1 )
        CALL dcopy( n, r2, 1, r1, 1 )
        CALL dcopy( n,  y, 1, r2, 1 )
        IF ( precon ) CALL msolve( n, r2, y )

        oldb   = beta                        ! oldb = betak
        beta   = ddot( n, r2, 1, y, 1 )    ! beta = betak+1^2
        IF (beta < zero) THEN
            istop = 6
            WRITE(*,*) 'unfortunately beta < 0 ',beta,zero
            exit
        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

        agamma = sqrt( gbar**2 + beta**2 )   ! gammak
        cs     = gbar / agamma               ! ck
        sn     = beta / agamma               ! 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', agamma
            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/agamma

        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, agamma)
        gmin   = min( gmin, agamma)
        z      = rhs1 / agamma
        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
        rnorm  = qrnorm

        ! 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
            IF (qrnorm <= epsx      ) istop = 2
            IF (qrnorm <= epsr      ) 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 (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

    END DO
    !     ------------------------------------------------------------------
    !     End of main iteration loop.
    !     ------------------------------------------------------------------

    ! Display final status.

900 IF (nout  > 0) THEN
        WRITE(nout, 2000) exit, istop, itn, exit, anorm, acond,  &
        exit, rnorm, ynorm
        WRITE(nout, 3000) exit, 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)
3000 format(      a, 5x, a )

END SUBROUTINE minres