lsqr.f90

比较经典的求解线性方程的方法 原理是C.C. Paige and M.A. Sauders等你提出的

      temp   = d2norm(alpha, beta)
      temp   = d2norm(temp , damp)
      Anorm  = d2norm(Anorm, temp)

      if (beta > zero) then
         call dscal (m, (one/beta), u, 1)
         call dscal (n, (- beta), v, 1)
         call Aprod (2,m,n,v,u,leniw,lenrw,iw,rw)
         alpha  = dnrm2 (n, v, 1)
         if (alpha > zero) then
            call dscal (n, (one/alpha), v, 1)
         end if
      end if

!     ------------------------------------------------------------------
!     Use a plane rotation to eliminate the damping parameter.
!     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
!     ------------------------------------------------------------------
      rhbar1 = rhobar
      if (damped) then
         rhbar1 = d2norm( rhobar, damp )
         cs1    = rhobar/ rhbar1
         sn1    = damp  / rhbar1
         psi    = sn1 * phibar
         phibar = cs1 * phibar
      end if                   

!     ------------------------------------------------------------------
!     Use a plane rotation to eliminate the subdiagonal element (beta)
!     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
!     ------------------------------------------------------------------
      rho    =   d2norm( rhbar1, beta )
      cs     =   rhbar1/rho
      sn     =   beta  /rho
      theta  =   sn*alpha
      rhobar = - cs*alpha
      phi    =   cs*phibar
      phibar =   sn*phibar
      tau    =   sn*phi

!     ------------------------------------------------------------------
!     Update  x, w  and (perhaps) the standard error estimates.
!     ------------------------------------------------------------------
      t1     =   phi/rho
      t2     = - theta/rho
      t3     =   one/rho
      dknorm =   zero

      if (wantse) then
         do i=1,n
            t      = w(i)
            x(i)   = t1*t +  x(i)
            w(i)   = t2*t +  v(i)
            t      =(t3*t)**2
            se(i)  = t + se(i)
            dknorm = t + dknorm
         end do
      else
         do i=1,n
            t      = w(i)
            x(i)   = t1*t + x(i)
            w(i)   = t2*t + v(i)
            dknorm =(t3*t)**2 + dknorm
         end do
      end if

!     ------------------------------------------------------------------
!     Monitor the norm of d_k, the update to x.
!     dknorm = norm( d_k )
!     dnorm  = norm( D_k ),        where   D_k = (d_1, d_2, ..., d_k )
!     dxk    = norm( phi_k d_k ),  where new x = x_k + phi_k d_k.
!     ------------------------------------------------------------------
      dknorm = sqrt(dknorm)
      dnorm  = d2norm(dnorm,dknorm)
      dxk    = abs(phi*dknorm)
      if (dxmax < dxk) then
          dxmax   =  dxk
          maxdx   =  itn
      end if

!     ------------------------------------------------------------------
!     Use a plane rotation on the right to eliminate the
!     super-diagonal element (theta) of the upper-bidiagonal matrix.
!     Then use the result to estimate  norm(x).
!     ------------------------------------------------------------------
      delta  =   sn2*rho
      gambar = - cs2*rho
      rhs    =   phi - delta*z
      zbar   =   rhs/gambar
      xnorm  =   d2norm(xnorm1,zbar)
      gamma  =   d2norm(gambar,theta)
      cs2    =   gambar/gamma
      sn2    =   theta /gamma
      z      =   rhs   /gamma
      xnorm1 =   d2norm(xnorm1,z)

!     ------------------------------------------------------------------
!     Test for convergence.
!     First, estimate the norm and condition of the matrix  Abar,
!     and the norms of  rbar  and  Abar(transpose)*rbar.
!     ------------------------------------------------------------------
      Acond  = Anorm*dnorm
      res2   = d2norm(res2,psi)
      rnorm  = d2norm(res2,phibar)
      Arnorm = alpha*abs( tau )

!     Now use these norms to estimate certain other quantities,
!     some of which will be small near a solution.

      alfopt = sqrt( rnorm/(dnorm*xnorm) )
      test1  = rnorm/bnorm
      test2  = zero
      if (rnorm > zero) test2 = Arnorm/(Anorm*rnorm)
      test3  = one/Acond
      t1     = test1/(one + Anorm*xnorm/bnorm)
      rtol   = btol  + atol*Anorm*xnorm/bnorm

!     The following tests guard against extremely small values of
!     atol, btol  or  ctol.  (The user may have set any or all of
!     the parameters  atol, btol, conlim  to zero.)
!     The effect is equivalent to the normal tests using
!     atol = relpr,  btol = relpr,  conlim = 1/relpr.

      t3 = one + test3
      t2 = one + test2
      t1 = one + t1
      if (itn >= itnlim) istop = 5
      if (t3  <= one   ) istop = 4
      if (t2  <= one   ) istop = 2
      if (t1  <= one   ) istop = 1

!     Allow for tolerances set by the user.

      if (test3 <= ctol) istop = 4
      if (test2 <= atol) istop = 2
      if (test1 <= rtol) istop = 1

!     ------------------------------------------------------------------
!     See if it is time to print something.
!     ------------------------------------------------------------------
      prnt = .false.
      if (nout > 0) then
         if (n     <=        40) prnt = .true.
         if (itn   <=        10) prnt = .true.
         if (itn   >= itnlim-10) prnt = .true.
         if (mod(itn,10)  ==  0) prnt = .true.
         if (test3 <=  2.0*ctol) prnt = .true.
         if (test2 <= 10.0*atol) prnt = .true.
         if (test1 <= 10.0*rtol) prnt = .true.
         if (istop /=         0) prnt = .true.

         if (prnt) then   ! Print a line for this iteration.
            write(nout,1500) itn,x(1),rnorm,test1,test2,Anorm,Acond,alfopt
         end if
      end if

!     ------------------------------------------------------------------
!     Stop if appropriate.
!     The convergence criteria are required to be met on  nconv
!     consecutive iterations, where  nconv  is set below.
!     Suggested value:  nconv = 1, 2  or  3.
!     ------------------------------------------------------------------
      if (istop == 0) then
         nstop = 0
      else
         nconv = 1
         nstop = nstop + 1
         if (nstop < nconv  .and.  itn < itnlim) istop = 0
      end if
      if (istop == 0) go to 100

!     ==================================================================
!     End of iteration loop.
!     ==================================================================

!     Finish off the standard error estimates.

      if (wantse) then
         t = one
         if (m > n)  t = m-n
         if (damped) t = m
         t = rnorm/sqrt(t)
      
         do i=1,n
            se(i) = t*sqrt(se(i))
         end do
      end if

!     Decide if istop = 2 or 3.
!     Print the stopping condition.

  800 if (damped .and. istop==2) istop=3
      if (nout > 0) then
         write(nout, 2000)              &
            exitt,istop,itn,            &
            exitt,Anorm,Acond,          &
            exitt,bnorm, xnorm,         &
            exitt,rnorm,Arnorm
         write(nout, 2100)              &
            exitt,dxmax, maxdx,         &
            exitt,dxmax/(xnorm+1.0d-20)
         write(nout, 3000)              &
            exitt, msg(istop)
      end if

  900 return

!     ------------------------------------------------------------------
 1000 format(// 1p, a, '     Least-squares solution of  Ax = b'   &
      / ' The matrix  A  has', i7, ' rows   and', i7, ' columns'  &
      / ' damp   =', e22.14, 3x,        'wantse =', l10           &
      / ' atol   =', e10.2, 15x,        'conlim =', e10.2         &
      / ' btol   =', e10.2, 15x,        'itnlim =', i10)
 1200 format(// '   Itn       x(1)           Function',           &
      '     Compatible   LS        Norm A    Cond A')
 1300 format(// '   Itn       x(1)           Function',           &
      '     Compatible   LS     Norm Abar Cond Abar alfa_opt')
 1500 format(1p, i6, 2e17.9, 4e10.2, e9.1)
 2000 format(/ 1p, a, 5x, 'istop  =', i2,   15x, 'itn    =', i8   &
      /     a, 5x, 'Anorm  =', e12.5, 5x, 'Acond  =', e12.5       &
      /     a, 5x, 'bnorm  =', e12.5, 5x, 'xnorm  =', e12.5       &
      /     a, 5x, 'rnorm  =', e12.5, 5x, 'Arnorm =', e12.5)
 2100 format(1p, a, 5x, 'max dx =', e8.1 , ' occurred at itn ', i8 &
      /     a, 5x, '       =', e8.1 , '*xnorm')
 3000 format(a, 5x, a)

      end subroutine LSQR

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
      function          d2norm( a, b )

      implicit          none
      double precision  d2norm, a, b

!-----------------------------------------------------------------------
!     d2norm  returns  sqrt( a**2 + b**2 )  with precautions
!     to avoid overflow.
!
!     21 Mar 1990: First version.
!     17 Sep 2007: Fortran 90 version.
!-----------------------------------------------------------------------

      intrinsic         abs, sqrt
      double precision  scale
      double precision, parameter :: zero = 0.0d+0

      scale = abs(a) + abs(b)
      if (scale == zero) then
         d2norm = zero
      else
         d2norm = scale*sqrt((a/scale)**2 + (b/scale)**2)
      end if

      end function d2norm