minresqlpModule.f90

Go to the documentation of this file.

 
 
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File minresqlpModule.f90
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
module minresqlpmodule
 
  use  minresqlpdatamodule, only : dp, ip, one, zero, eps, realmin, prcsn, debug
  use  minresqlpblasmodule, only : dnrm2
 
  implicit none
 
  private                                       ! sets default for module
  public   :: minresqlp, symortho
 
contains
 
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
  subroutine minresqlp( n, Aprod, b, shift, Msolve, disable, nout,        &
                        itnlim, rtol, maxxnorm, trancond, Acondlim,       &
                        x, istop, itn, rnorm, Arnorm, xnorm, Anorm, Acond )
    ! Inputs
    integer(ip), intent(in)             :: n
    real(dp),    intent(in)             :: b(n)
    integer(ip), intent(in), optional   :: itnlim, nout
    logical,     intent(in), optional   :: disable
    real(dp),    intent(in), optional   :: shift
    real(dp),    intent(in), optional   :: rtol, maxxnorm, trancond, acondlim
 
    ! Outputs
    real(dp),    intent(out)            :: x(n)
    integer(ip), intent(out), optional  :: istop, itn
    real(dp),    intent(out), optional  :: rnorm, arnorm, xnorm, anorm, acond
 
    optional :: msolve
 
    interface
       subroutine aprod(n,x,y)                    ! y := A*x
         use minresqlpdatamodule
         integer(ip), 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 minresqlpdatamodule
         integer(ip), intent(in)    :: n
         real(dp),    intent(in)    :: x(n)
         real(dp),    intent(out)   :: y(n)
       end subroutine msolve
    end interface
 
    intrinsic :: abs, sqrt, present, floor, log10, dot_product
 
    ! Local arrays and variables
    real(dp)     :: shift_
    real(dp)     :: rtol_, maxxnorm_, trancond_, acondlim_
    real(dp)     :: rnorm_, arnorm_, xnorm_, anorm_, acond_
    logical      :: checka_, precon_, disable_
    integer(ip)  :: itnlim_, nout_, istop_, itn_
 
    real(dp)     :: r1(n), r2(n), v(n), w(n), wl(n), wl2(n),&
                    xl2(n), y(n), vec2(2), vec3(3)
    real(dp)     :: abs_gama, acondl, alfa, anorml, arnorml,&
                    axnorm, beta, beta1, betal, betan,      &
                    epsa, epsx, gminl2, ieps, pnorm,        &
                    relares, relaresl, relres, relresl,     &
                    rnorml, rootl, t1, t2, xl2norm,         &
                    xnorm_tmp, xnorml, z,                   &
                    cr1, cr2, cs, dbar, dlta, dlta_qlp,     &
                    dlta_tmp, dltan, epln, eplnn, eta, etal,&
                    etal2, gama, gama_qlp, gama_tmp, gamal, &
                    gamal_qlp, gamal_tmp, gamal2, gamal3,   &
                    gbar, gmin, gminl, phi, s, sn, sr1, sr2,&
                    t, tau, taul, taul2, u, u_qlp,          &
                    ul, ul_qlp, ul2, ul2_qlp, ul3, ul4,     &
                    vepln, vepln_qlp, veplnl, veplnl2, x1last
 
    integer(ip)  :: j, qlpiter, headlines, lines, nprint, flag0, ios
    logical      :: prnt, done, lastiter, connected, named_file, likels
    character(len=20) :: filename
    character(len=2)  :: qlpstr = ' '
 
    ! Local constants
    real(dp), parameter :: epsinv  = 10.0_dp**floor(log10(one/eps))
    real(dp)            :: normmax = 10.0_dp**floor(log10(one/eps)/2)
    character(len=*), parameter :: enter = ' Enter MINRES-QLP.  '
    character(len=*), parameter :: exitt = ' Exit  MINRES-QLP.  '
    character(len=*), parameter :: msg(1:15) =                                &
      (/ 'beta_{k+1} < eps.                                                ', & ! 1
         'beta2 = 0.  If M = I, b and x are eigenvectors of A.             ', & ! 2
         'beta1 = 0.  The exact solution is  x = 0.                        ', & ! 3
         'A solution to (poss. singular) Ax = b found, given rtol.         ', & ! 4
         'A solution to (poss. singular) Ax = b found, given eps.          ', & ! 5
         'Pseudoinverse solution for singular LS problem, given rtol.      ', & ! 6
         'Pseudoinverse solution for singular LS problem, given eps.       ', & ! 7
         'The iteration limit was reached.                                 ', & ! 8
         'The operator defined by Aprod appears to be unsymmetric.         ', & ! 9
         'The operator defined by Msolve appears to be unsymmetric.        ', & ! 10
         'The operator defined by Msolve appears to be indefinite.         ', & ! 11
         'xnorm has exceeded maxxnorm or will exceed it next iteration.    ', & ! 12
         'Acond has exceeded Acondlim or 0.1/eps.                          ', & ! 13
         'Least-squares problem but no converged solution yet.             ', & ! 14
         'A null vector obtained, given rtol.                              ' /) ! 15
 
     character(len=*), parameter :: ddebugstr1 = "(a, T5, i0, a, 5(e12.3))"
     character(len=*), parameter :: ddebugstr2 = "(5(a, i0, a, e12.3, a))"
 
     character(len=*), parameter :: headerstr =                    &
        "(// 1p,    a, 4x, 'Solution of symmetric   Ax = b'"    // &
        "  / ' n        =', i7, 6x,  '||b||    =', e11.2, 3x,"  // &
        "     'precon   =', l4                           "      // &
        "  / ' itnlim   =', i7, 6x, 'rtol     =', e11.2, 3x,"   // &
        "     'shift    =', e23.14                       "      // &
        "  / ' maxxnorm =', e11.2, 2x, 'Acondlim =', e11.2, 3x,"// &
        "     'trancond =', e11.2)"
     character(len=*), parameter :: tableheaderstr =                         &
      "(// '    iter   x(1)              xnorm     rnorm     Arnorm   '," // &
      " 'Compatible   LS      norm(A)   cond(A)')"
     character(len=*), parameter :: itnstr = "(1p, i8, e19.10, 7e10.2, a)"
     character(len=*), parameter :: finalstr1 =                    &
        "(/ 1p, a, 5x, a, i3,   14x, a, i8 "                    // &
        " /     a, 5x, a, e12.4, 5x, a, e12.4 "                 // &
        " /     a, 5x, a, e12.4, 5x, a, e12.4 "                 // &
        " /     a, 5x, a, e12.4 )"
     character(len=*), parameter :: finalstr2 = "(      a, 5x, a )"
 
    !------------------------------------------------------------------
    prnt = .false.
    t1 = zero
    t2 = zero
    gminl = zero
    ! Optional inputs
    if (present(shift)) then
       shift_ = shift
    else
       shift_ = zero
    end if
 
    checka_ = .true.
 
    if (present(disable)) then
       disable_ = disable
    else
       disable_ = .false.
    end if
 
    if (present(itnlim)) then
       itnlim_ = itnlim
    else
       itnlim_ = 4 * n
    end if
 
    connected = .false.
    filename = "MINRESQLP_tmp.txt"
    nout_ = 10
 
    if (present(nout)) then
       nout_ = nout
       inquire(unit=nout, opened=connected, named=named_file, name=filename)
       !write(*,*) connected, named_file, filename
       if (.not. connected) then
          write(*,*) "File unit 'nout' is not open."
          if (nout==5  .or.  nout == 6) then
             nout_ = 10
          end if
       end if
    end if
 
    if (.not. connected) then
        write(*,*) 'nout_ = ', nout_
        open(nout_, file=filename, status='unknown', iostat=ios)
        write(*,*) 'ios = ', ios
        if (ios /= 0) then
           write(*,*) "Error opening file '", filename, "'."
           stop
        end if
    end if
 
    if (present(rtol)) then
       rtol_ = rtol
    else
       rtol_ = eps
    end if
 
    if (prcsn == 6) then
       normmax = 1.0e4_dp
    end if
    if (present(maxxnorm)) then
       maxxnorm_ = min(maxxnorm, one/eps)
    else
       maxxnorm_ = normmax
    end if
 
    if (present(trancond)) then
       trancond_ = min(trancond, normmax)
    else
       trancond_ = normmax
    end if
 
    if (present(acondlim)) then
       acondlim_ = min(acondlim, epsinv)
    else
       acondlim_ = epsinv
    end if
 
    if (present(msolve)) then
       precon_ = .true.
    else
       precon_ = .false.
    end if
 
    !------------------------------------------------------------------
    ! Print heading and initialize.
    !------------------------------------------------------------------
    nprint   = min(n,20)
    !debug   = .true.
    lastiter = .false.
    flag0    = 0
    istop_   = flag0
    beta1    = dnrm2(n, b, 1)
    ieps     = 0.1_dp/eps
    itn_     = 0
    qlpiter  = 0
    xnorm_   = zero
    xl2norm  = zero
    axnorm   = zero
    anorm_   = zero
    acond_   = one
    pnorm    = zero
    relresl  = zero
    relaresl = zero
    x        = zero
    xl2      = zero
    x1last   = x(1)
 
    if (nout_ > 0) then
       write(nout_, headerstr) enter, n, beta1, precon_, itnlim_, rtol_, &
       shift_, maxxnorm_, acondlim_, trancond_
    end if
 
    !------------------------------------------------------------------
    ! 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_ ) then
       call msolve( n, b, y )
    end if
 
    beta1  = dot_product(b, y)
 
    if (beta1 < zero .and. dnrm2(n, y, 1) > eps) then     ! M must be indefinite.
       istop_ = 11
    end if
 
    if (beta1 == zero) then    ! b = 0 exactly.  Stop with x = 0.
       istop_ = 3
    end if
 
    beta1  = sqrt( beta1 )     ! Normalize y to get v1 later.
 
    if (debug) then
       write(*,ddebugstr1) ' y_', itn_, ' = ', (y(j), j=1,nprint)
       write(*,*) 'beta1 ', beta1
    end if
 
    !------------------------------------------------------------------
    ! 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   = (abs(s) + eps) * eps**0.33333_dp
       if (z > epsa) then
          istop_ = 10
       end if
    end if
 
    !------------------------------------------------------------------
    ! See if Aprod  is symmetric.
    !------------------------------------------------------------------
    if (checka_) then
       call aprod ( n, y, w  )  ! w  = A*y
       call aprod ( n, w, r2 )  ! r2 = A*w
       s      = dot_product( w, w )
       t      = dot_product( y, r2)
       z      = abs( s - t )
       epsa   = (abs(s) + eps) * eps**0.33333_dp
       if (z > epsa) then
          istop_ = 9
       end if
    end if
 
    !------------------------------------------------------------------
    ! Initialize other quantities.
    !------------------------------------------------------------------
    tau    = zero
    taul   = zero
    gmin   = zero
    beta   = zero
    betan  = beta1
    dbar   = zero
    dltan  = zero
    eta    = zero
    etal   = zero
    etal2  = zero
    vepln  = zero
    veplnl = zero
    veplnl2= zero
    eplnn  = zero
    gama   = zero
    gamal  = zero
    gamal2 = zero
    phi    = beta1
    cr1    = -one
    sr1    = zero
    cr2    = -one
    sr2    = zero
    cs     = -one
    sn     = zero
    ul3    = zero
    ul2    = zero
    ul     = zero
    u      = zero
    lines  = 1
    headlines = 30 * lines
    rnorm_ = betan
    relres = rnorm_ / (anorm_*xnorm_ + beta1)
    relares= zero
    r2     = b
    w      = zero    ! vector of zeros
    wl     = zero
    done   = .false.
 
    if (debug) then
       write(*,*)
       write(*,*) 'Checking variable values before main loop'
       write(*,*) 'istop ', istop_, ' done ', done, ' itn ', itn_, ' QLPiter' , qlpiter
       write(*,*) 'lastiter', lastiter, ' lines ', lines, ' headlines ', headlines
       write(*,*) 'beta ', beta, ' tau ', tau, ' taul ', taul, ' phi ', phi
       write(*,*) 'betan ', betan, ' gmin ', gmin, ' cs ', cs, ' sn ', sn
       write(*,*) 'cr1 ', cr1, ' sr1 ', sr1, ' cr2 ', cr2, ' sr2 ', sr2
       write(*,*) 'dltan ', dltan, ' eplnn ', eplnn, ' gama ', gama, ' gamal ', gamal
       write(*,*) 'gamal2 ', gamal2, ' eta ', eta, ' etal ', etal, 'etal2', etal2
       write(*,*) 'vepln ', vepln, ' veplnl', veplnl, ' veplnl2 ', veplnl2, ' ul3 ', ul3
       write(*,*) 'ul2 ', ul2, ' ul ', ul, ' u ', u, ' rnorm ', rnorm_
       write(*,*) 'xnorm ', xnorm_, ' xl2norm ', xl2norm, ' pnorm ', pnorm, ' Axnorm ', axnorm
       write(*,*) 'Anorm ', anorm_, ' Acond ', acond_, ' relres ', relres
       write(*,ddebugstr1) 'w_', itn_-1, ' = ', (wl(j), j=1,nprint)
       write(*,ddebugstr1) 'w_', itn_,   ' = ', (w(j),  j=1,nprint)
       write(*,ddebugstr1) 'x_', itn_,   ' = ', (x(j),  j=1,nprint)
    end if
 
 
    !------------------------------------------------------------------
    ! Main iteration loop.
    !------------------------------------------------------------------
    do while (istop_ <= flag0)
       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.
       !-----------------------------------------------------------------
 
       betal  = beta;               ! betal = betak
       beta   = betan;
       s      = one / beta          ! Normalize previous vector (in y).
       v      = s*y;                ! v = vk if P = I.
       call aprod ( n, v, y )
       if (shift_ /= zero) then
          y = y - shift_ * v
       end if
       if (itn_ >= 2) then
          y = y + (- beta/betal) * r1
       end if
       alfa = dot_product(v, y)                  ! alphak
       y    = y + (- alfa/beta) * r2
       r1   = r2
       r2   = y
 
       if ( .not. precon_ ) then
          betan = dnrm2(n, y, 1)                 ! betan = ||y||_2
       else
          call msolve( n, r2, y )
          betan = dot_product(r2, y)             ! betan = betak+1^2
          if (betan > zero) then
             betan = sqrt(betan)
          elseif ( dnrm2(n, y, 1) > eps ) then   ! M must be indefinite.
             istop_ = 11
             exit
          end if
       end if
 
       if (itn_ == 1) then
          vec2(1) = alfa
          vec2(2) = betan
          pnorm   = dnrm2(2, vec2, 1)
       else
          vec3(1) = beta
          vec3(2) = alfa
          vec3(3) = betan
          pnorm   = dnrm2(3, vec3, 1)
       end if
 
 
       if (debug) then
          write(*,*)
          write(*,*) 'Lanczos iteration ', itn_
          write(*,ddebugstr1) 'v_',  itn_, ' = ', (v(j),  j=1,nprint)
          write(*,ddebugstr1) 'r1_', itn_, ' = ', (r1(j), j=1,nprint)
          write(*,ddebugstr1) 'r2_', itn_, ' = ', (r2(j), j=1,nprint)
          write(*,ddebugstr1) 'y_',  itn_, ' = ', (y(j),  j=1,nprint)
          write(*,ddebugstr2) 'alpha_', itn_,   ' = ', alfa,  ', ', &
                              'beta_',  itn_,   ' = ', beta,  ', ', &
                              'beta_',  itn_+1, ' = ', betan, ', ', &
                              'pnorm_', itn_,   ' = ', pnorm
       end if
 
       ! Apply previous left reflection Qk-1 to get
       !   [deltak epslnk+1] = [cs  sn][dbark    0   ]
       !   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
 
       dbar   = dltan
       dlta   = cs * dbar  +  sn * alfa  ! dlta1  = 0         deltak
       epln   = eplnn;
       gbar   = sn * dbar  -  cs * alfa  ! gbar 1 = alfa1     gbar k
       eplnn  =               sn * betan ! eplnn2 = 0         epslnk+1
       dltan  =            -  cs * betan ! dbar 2 = beta2     dbar k+1
       dlta_qlp = dlta;
 
       if (debug) then
          write(*,*)
          write(*,*) 'Apply previous left reflection Q_{', itn_-1, ',', itn_, '}'
          write(*,ddebugstr2) 'c_',    itn_-1, ' = ', cs,    ', ', &
                              's_',    itn_-1, ' = ', sn
          write(*,ddebugstr2) 'dlta_', itn_,   ' = ', dlta,  ', ', &
                              'gbar_', itn_,   ' = ', gbar,  ', ', &
                              'epln_', itn_+1, ' = ', eplnn, ', ', &
                              'dbar_', itn_+1, ' = ', dltan
       end if
 
       ! Compute the current left reflection Qk
       gamal3 = gamal2
       gamal2 = gamal
       gamal  = gama
       call symortho(gbar, betan, cs, sn, gama)
       gama_tmp = gama;
       taul2  = taul
       taul   = tau
       tau    = cs * phi
       phi    = sn * phi                     ! phik
       axnorm = sqrt( axnorm**2 + tau**2 );
 
       if (debug) then
          write(*,*)
          write(*,*) 'Compute the current left reflection Q_{', itn_, ',', itn_+1, '}'
          write(*,ddebugstr2) 'c_',    itn_, ' = ', cs,   ', ', &
                              's_',    itn_, ' = ', sn
          write(*,ddebugstr2) 'tau_',  itn_, ' = ', tau,  ', ', &
                              'phi_',  itn_, ' = ', phi,  ', ', &
                              'gama_', itn_, ' = ', gama
       end if
 
       ! Apply the previous right reflection P{k-2,k}
 
       if (itn_ > 2) then
          veplnl2  = veplnl
          etal2    = etal
          etal     = eta
          dlta_tmp = sr2 * vepln - cr2 * dlta
          veplnl   = cr2 * vepln + sr2 * dlta
          dlta     = dlta_tmp
          eta      = sr2 * gama
          gama     = -cr2 * gama
          if (debug) then
             write(*,*)
             write(*,*) 'Apply the previous right reflection P_{', itn_-2, ',', itn_, '}'
             write(*,ddebugstr2) 'cr2_',    itn_,   ' = ', cr2,    ', ', &
                                 'sr2_',    itn_,   ' = ', sr2
             write(*,ddebugstr2) 'gamal2_', itn_,   ' = ', gamal2, ', ', &
                                 'gamal_',  itn_,   ' = ', gamal,  ', ', &
                                 'gama_',   itn_,   ' = ', gama
             write(*,ddebugstr2) 'dlta_',   itn_,   ' = ', dlta,   ', ', &
                                 'vepln_',  itn_-1, ' = ', veplnl
          end if
       end if
 
 
       ! Compute the current right reflection P{k-1,k}, P_12, P_23,...
       if (itn_ > 1) then
          call symortho(gamal, dlta, cr1, sr1, gamal_tmp)
          gamal     = gamal_tmp
          vepln     = sr1 * gama
          gama      = -cr1 * gama
 
          if (debug) then
             write(*,*)
             write(*,*) 'Apply the current right reflection  P_{', itn_-1, ',', itn_, '}'
             write(*,ddebugstr2) 'cr1_ ',  itn_,   ' = ', cr1,    ', ', &
                                 'sr1_' ,  itn_,   ' = ', sr1
             write(*,ddebugstr2) 'gama_',  itn_-2, ' = ', gamal2, ', ', &
                                 'gama_',  itn_-1, ' = ', gamal,  ', ', &
                                 'gama_',  itn_,   ' = ', gama
             write(*,ddebugstr2) 'dlta_',  itn_,   ' = ', dlta,   ', ', &
                                 'vepln_', itn_-1, ' = ', veplnl, ', ', &
                                 'eta_',   itn_,   ' = ', eta
          end if
       end if
 
       ! Update xnorm
 
       xnorml  = xnorm_
       ul4     = ul3
       ul3     = ul2
 
       if (itn_ > 2) then
          ul2 = ( taul2 - etal2 * ul4 - veplnl2 * ul3 ) / gamal2
          if (debug) then
             write(*,ddebugstr2) 'tau_',   itn_-2, ' = ', taul2,   ', ', &
                                 'eta_',   itn_-2, ' = ', etal2,   ', ', &
                                 'vepln_', itn_-2, ' = ', veplnl2, ', ', &
                                 'gama_',  itn_-2, ' = ', gamal2
          end if
       end if
       if (itn_ > 1) then
          ul  = ( taul  - etal  * ul3 - veplnl  * ul2) / gamal
          if (debug) then
             write(*,ddebugstr2) 'tau_',   itn_-1, ' = ', taul,   ', ', &
                                 'eta_',   itn_-1, ' = ', etal,   ', ', &
                                 'vepln_', itn_-1, ' = ', veplnl, ', ', &
                                 'gamal_', itn_-1, ' = ', gamal
          end if
       end if
 
       vec3(1) = xl2norm
       vec3(2) = ul2
       vec3(3) = ul
       xnorm_tmp = dnrm2(3, vec3, 1)     ! norm([xl2norm ul2 ul]);
       if (abs(gama) > eps) then         ! .and.  xnorm_tmp < maxxnorm_) then
          if (debug) then
             write(*,ddebugstr2) 'tau_',   itn_, ' = ', tau,   ', ', &
                                 'eta_',   itn_, ' = ', eta,   ', ', &
                                 'vepln_', itn_, ' = ', vepln, ', ', &
                                 'gama_',  itn_, ' = ', gama
          end if
          u       = (tau - eta*ul2 - vepln*ul) / gama
          likels  = relaresl < relresl
          vec2(1) = xnorm_tmp
          vec2(2) = u
          if (likels  .and.  dnrm2(2, vec2, 1) > maxxnorm_) then
             u      = zero
             istop_ = 12
          end if
       else
          u      = zero
          istop_ = 14
       end if
       vec2(1)   = xl2norm
       vec2(2)   = ul2
       xl2norm   = dnrm2(2, vec2, 1)
       vec3(1)   = xl2norm
       vec3(2)   = ul
       vec3(3)   = u
       xnorm_    = dnrm2(3, vec3, 1)
 
       if (acond_ < trancond_  .and.  istop_ == flag0  .and.  qlpiter == 0) then !! MINRES updates
          wl2  = wl
          wl   = w
          if (gama_tmp > eps) then
             s = one / gama_tmp
             w = (v - epln*wl2 - dlta_qlp*wl) * s
          end if
 
          if (xnorm_ < maxxnorm_) then
             x1last = x(1)
             x = x + tau*w
          else
             istop_   = 12
             lastiter = .true.
          end if
 
          if (debug) then
             write(*,*)
             write(*,*) 'MINRES updates'
             write(*,ddebugstr2) 'gama_tmp_', itn_, ' = ', gama_tmp, ', ', &
                                 'tau_',      itn_, ' = ', tau,      ', ', &
                                 'epln_',     itn_, ' = ', epln,     ', ', &
                                 'dlta_QLP_', itn_, ' = ', dlta_qlp
             write(*,ddebugstr1) 'v_', itn_ , ' = ', (v(j), j=1,nprint)
             write(*,ddebugstr1) 'w_', itn_ , ' = ', (w(j), j=1,nprint)
          end if
       else                                         !! MINRES-QLP updates
          qlpiter = qlpiter + 1;
          if (qlpiter == 1) then
             xl2 = zero     ! vector
             if (itn_ > 1) then                     ! construct w_{k-3}, w_{k-2}, w_{k-1}
                if (itn_ > 3) then
                   wl2 = gamal3*wl2 + veplnl2*wl + etal*w
                end if                              ! w_{k-3}
                if (itn_ > 2) then
                   wl = gamal_qlp*wl + vepln_qlp*w
                end if                              ! w_{k-2}
                w = gama_qlp*w
                xl2 =  x - ul_qlp*wl - u_qlp*w
             end if
          end if
 
          if (itn_ == 1) then
             wl2 =  wl
             wl  =  sr1*v
             w   = -cr1*v
          else if (itn_ == 2) then
             wl2 = wl
             wl  = cr1*w + sr1*v
             w   = sr1*w - cr1*v
          else
             wl2 = wl
             wl  = w
             w   = sr2*wl2 - cr2*v
             wl2 = cr2*wl2 + sr2*v
             v   = cr1*wl  + sr1*w
             w   = sr1*wl  - cr1*w
             wl  = v
          end if
          x1last = x(1)
          xl2 = xl2 + ul2*wl2
          x   = xl2 + ul *wl + u*w
 
          if (debug) then
             write(*,*)
             write(*,*) 'MINRESQLP updates'
          end if
       end if
 
       if (debug) then
          write(*,*)
          write(*,*) 'Update u, w, x and xnorm'
          write(*,ddebugstr2) 'u_', itn_-2, ' = ', ul2, ', ', &
                              'u_', itn_-1, ' = ', ul,  ', ', &
                              'u_', itn_,   ' = ', u
          write(*,ddebugstr1) 'w_',   itn_-2, ' = ',  (wl2(j), j=1,nprint)
          write(*,ddebugstr1) 'w_',   itn_-1, ' = ',  (wl(j),  j=1,nprint)
          write(*,ddebugstr1) 'w_',   itn_,   ' = ',  (w(j),   j=1,nprint)
          write(*,ddebugstr1) 'x_',   itn_,   ' = ',  (x(j),   j=1,nprint)
          write(*,ddebugstr2) 'xnorm_', itn_-2, ' = ', xl2norm, ', ', &
                              'xnorm_', itn_-1, ' = ', xnorml,  ', ', &
                              'xnorm_', itn_,   ' = ', xnorm_
       end if
 
       ! Compute the next right reflection P{k-1,k+1}
 
       if (debug) then
          write(*,*)
          write(*,*) 'Compute the next right reflection P{', itn_-1, itn_+1,'}'
          write(*,ddebugstr2) 'gama_', itn_-1, ' = ', gamal, ', ', &
                              'epln_', itn_+1, ' = ', eplnn
       end if
 
       gamal_tmp = gamal
       call symortho(gamal_tmp, eplnn, cr2, sr2, gamal)
 
       if (debug) then
          write(*,ddebugstr2) 'cr2_',  itn_+1, ' = ', cr2, ', ', &
                              'sr2_',  itn_+1, ' = ', sr2, ', ', &
                              'gama_', itn_-1, ' = ', gamal
       end if
 
       ! Store quantities for transfering from MINRES to MINRES-QLP
 
       gamal_qlp = gamal_tmp
       vepln_qlp = vepln
       gama_qlp  = gama
       ul2_qlp   = ul2
       ul_qlp    = ul
       u_qlp     = u
 
       if (debug) then
          write(*,*)
          write(*,*) 'Store quantities for transfering from MINRES to MINRES-QLP '
          write(*,ddebugstr2) 'gama_QLP_', itn_-1, ' = ', gamal_qlp, ', ', &
                              'vepln_QLP_',itn_,   ' = ', vepln_qlp, ', ', &
                              'gama_QLP_', itn_,   ' = ', gama_qlp
          write(*,ddebugstr2) 'u_QLP_', itn_-2,    ' = ', ul2_qlp,   ', ', &
                              'u_QLP_', itn_-1,    ' = ', ul_qlp,    ', ', &
                              'u_QLP_', itn_,      ' = ', u_qlp
       end if
 
       ! Estimate various norms
 
       abs_gama = abs(gama)
       anorml   = anorm_
       anorm_   = max(anorm_, pnorm, gamal, abs_gama)
       if (itn_ == 1) then
          gmin  = gama
          gminl = gmin
       else if (itn_ > 1) then
          gminl2  = gminl
          gminl   = gmin
          vec3(1) = gminl2
          vec3(2) = gamal
          vec3(3) = abs_gama
          gmin    = min(gminl2, gamal, abs_gama)
       end if
       acondl   = acond_
       acond_   = anorm_ / gmin
       rnorml   = rnorm_
       relresl  = relres
       if (istop_ /= 14) rnorm_ = phi
       relres   = rnorm_ / (anorm_ * xnorm_ + beta1)
       vec2(1)  = gbar
       vec2(2)  = dltan
       rootl    = dnrm2(2, vec2, 1)
       arnorml  = rnorml * rootl
       relaresl = rootl / anorm_
 
       if (debug) then
          write(*,*)
          write(*,*) 'Estimate various norms '
          write(*,ddebugstr2) 'gmin_',   itn_,   ' = ', gmin,     ', ', &
                              'pnorm_',  itn_,   ' = ', pnorm,    ', ', &
                              'rnorm_',  itn_,   ' = ', rnorm_,   ', ', &
                              'Arnorm_', itn_-1, ' = ', arnorml
          write(*,ddebugstr2) 'Axnorm_', itn_,   ' = ', axnorm,   ', ', &
                              'Anorm_',  itn_,   ' = ', anorm_,   ', ', &
                              'Acond_',  itn_,   ' = ', acond_
       end if
 
       ! See if any of the stopping criteria are satisfied.
 
       epsx = anorm_*xnorm_*eps
       if (istop_ == flag0  .or.  istop_ == 14) then
          t1   = one + relres
          t2   = one + relaresl
       end if
       if (t1 <= one                  ) then
          istop_ = 5                           ! Accurate Ax=b solution
       else if (t2 <= one             ) then
          istop_ = 7                           ! Accurate LS solution
       else if (relres   <= rtol_     ) then
          istop_ = 4                           ! Good enough Ax=b solution
       else if (relaresl <= rtol_     ) then
          istop_ = 6                           ! Good enough LS solution
       else if (epsx     >= beta1     ) then
          istop_ = 2                           ! x is an eigenvector
       else if (xnorm_    >= maxxnorm_) then
          istop_ = 12                          ! xnorm exceeded its limit
       else if (acond_    >= acondlim_ .or. acond_ >= ieps) then
          istop_ = 13                          ! Huge Acond
       else if (itn_      >= itnlim_  ) then
          istop_ = 8                           ! Too many itns
       else if (betan     <  eps      ) then
          istop_ = 1                           ! Last iteration of Lanczos, rarely happens
       end if
 
       if (disable_  .and.  itn_ < itnlim_) then
          istop_ = flag0
          done   = .false.
          if (axnorm < rtol_*anorm_*xnorm_) then
             istop_   = 15
             lastiter = .false.
          end if
       end if
 
       if (istop_ /= flag0) then
          done = .true.
          if (istop_ == 6  .or.  istop_ == 7  .or.  istop_ == 12  .or.  istop_ == 13) then
             lastiter = .true.
          end if
          if (lastiter) then
             itn_   = itn_ - 1
             acond_ = acondl
             rnorm_ = rnorml
             relres = relresl
          end if
 
          call aprod ( n, x, r1 )
          r1  = b - r1 + shift_*x       ! r1 to temporarily store residual vector
          call aprod ( n, r1, wl2 )     ! wl2 to temporarily store A*r1
          wl2 = wl2 - shift_*r1
          arnorm_ = dnrm2(n, wl2, 1)
          if (rnorm_ > zero  .and.  anorm_ > zero) then
             relares = arnorm_ / (anorm_*rnorm_)
          end if
       end if
 
       if (nout_ > 0  .and.  .not. lastiter  .and.  mod(itn_-1,lines) == 0) then
          if (itn_ == 101) then
             lines     = 10
             headlines = 30*lines
          else if (itn_ == 1001) then
             lines = 100
             headlines = 30*lines
          end if
 
          if (qlpiter == 1) then
              qlpstr = ' P'
          else
              qlpstr = ' '
          end if
       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_-1,10) == 0           ) prnt = .true.
          if (qlpiter        == 1           ) prnt = .true.
          if (acond_         >= 0.01_dp/eps ) prnt = .true.
          if (istop_         /= flag0       ) prnt = .true.
 
          if ( prnt ) then
             if (itn_ == 1) write(nout_, tableheaderstr)
             write(nout_, itnstr) itn_-1, x1last, xnorml, &
             rnorml, arnorml, relresl, relaresl, anorml, acondl, qlpstr
             if (mod(itn_,10) == 0) write(nout_, "(1x)")
          end if
       end if
 
       if (debug) then
          write(*,*) 'istop = ', istop_
       end if
       if (istop_ /= flag0) exit
    enddo
    !===================================================================
    ! End of iteration loop.
    !===================================================================
 
    ! Optional outputs
 
    if (present(istop)) then
       istop = istop_
    end if
 
    if (present(itn)) then
       itn = itn_
    end if
 
    if (present(rnorm)) then
       rnorm = rnorm_
    end if
 
    if (present(arnorm)) then
       arnorm = arnorm_
    end if
 
    if (present(xnorm)) then
       xnorm = xnorm_
    end if
 
    if (present(anorm)) then
       anorm = anorm_
    end if
 
    if (present(acond)) then
       acond = acond_
    end if
 
    if ( prnt ) then
       write(nout_, itnstr) itn_, x(1), xnorm_, &
             rnorm_, arnorm_, relres, relares, anorm_, acond_, qlpstr
    end if
 
    ! Display final status.
 
    if (nout_ > 0) then
       write(nout_, finalstr1) exitt, 'istop =', istop_, 'itn    =', itn_,    &
                               exitt, 'Anorm =', anorm_, 'Acond  =', acond_,  &
                               exitt, 'rnorm =', rnorm_, 'Arnorm =', arnorm_, &
                               exitt, 'xnorm =', xnorm_
       write(nout_, finalstr2) exitt, msg(istop_)
    end if
 
    return
end subroutine minresqlp
 
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  subroutine symortho(a, b, c, s, r)
 
    real(dp), intent(in)  :: a, b
    real(dp), intent(out) :: c, s, r
 
    intrinsic :: abs, sqrt
 
    ! Local variables
    logical, parameter :: debug = .false.
    real(dp)           :: t
    real(dp)           :: abs_a, abs_b
 
    abs_a = abs(a)
    abs_b = abs(b)
 
    if (abs_b <= realmin) then
       s = zero
       r = abs_a
       if (a == zero) then
          c = one
       else
          c = a / abs_a
       end if
 
    else if (abs_a <= realmin) then
       c = zero;
       r = abs_b
       s = b / abs_b
 
    else if (abs_b > abs_a) then
       t = a / b
       s = (b / abs_b) / sqrt(one + t**2)
       c = s * t
       r = b / s  ! computationally better than r = a / c since |c| <= |s|
 
    else
       t = b / a
       c = (a / abs_a) / sqrt(one + t**2)
       s = c * t;
       r = a / c  ! computationally better than r = b / s since |s| <= |c|
    end if
 
    if (debug) then
       write(*,*)  'c = ',  c, ', s = ',  s, ', r = ',  r
    end if
 
  end subroutine symortho
 
end module minresqlpmodule