lsqr.f90

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

!                        not compatible.  A least-squares solution has
!                        been obtained that is sufficiently accurate,
!                        given the value of atol.
!
!                3       damp is nonzero.  A damped least-squares
!                        solution has been obtained that is sufficiently
!                        accurate, given the value of atol.
!
!                4       An estimate of cond(Abar) has exceeded
!                        conlim.  The system A*x = b appears to be
!                        ill-conditioned.  Otherwise, there could be an
!                        error in subroutine Aprod.
!
!                5       The iteration limit itnlim was reached.
!
!     itn     output     The number of iterations performed.
!
!     Anorm   output     An estimate of the Frobenius norm of  Abar.
!                        This is the square-root of the sum of squares
!                        of the elements of Abar.
!                        If damp is small and if the columns of A
!                        have all been scaled to have length 1.0,
!                        Anorm should increase to roughly sqrt(n).
!                        A radically different value for Anorm may
!                        indicate an error in subroutine Aprod (there
!                        may be an inconsistency between modes 1 and 2).
!
!     Acond   output     An estimate of cond(Abar), the condition
!                        number of Abar.  A very high value of Acond
!                        may again indicate an error in Aprod.
!
!     rnorm   output     An estimate of the final value of norm(rbar),
!                        the function being minimized (see notation
!                        above).  This will be small if A*x = b has
!                        a solution.
!
!     Arnorm  output     An estimate of the final value of
!                        norm( Abar(transpose)*rbar ), the norm of
!                        the residual for the usual normal equations.
!                        This should be small in all cases.  (Arnorm
!                        will often be smaller than the true value
!                        computed from the output vector x.)
!
!     xnorm   output     An estimate of the norm of the final
!                        solution vector x.
!
!     Subroutines and functions used              
!     ------------------------------
!     USER               Aprod
!     LSQR               d2norm
!     BLAS               dcopy, dnrm2, dscal (see Lawson et al. below)
!
!
!     Precision
!     ---------
!     The number of iterations required by LSQR will usually decrease
!     if the computation is performed in higher precision.
!     At least 15-digit arithmetic should normally be used.
!     To convert LSQR and d2norm between single and double precision,
!     change
!                        double precision
!                        dcopy, dnrm2, dscal
!     to the appropriate FORTRAN and BLAS equivalents.
!     Also change 'd+' or 'e+' in the parameter statement.
!
!
!     References
!     ----------
!     C.C. Paige and M.A. Saunders,  LSQR: An algorithm for sparse
!          linear equations and sparse least squares,
!          ACM Transactions on Mathematical Software 8, 1 (March 1982),
!          pp. 43-71.
!
!     C.C. Paige and M.A. Saunders,  Algorithm 583, LSQR: Sparse
!          linear equations and least-squares problems,
!          ACM Transactions on Mathematical Software 8, 2 (June 1982),
!          pp. 195-209.
!
!     C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh,
!          Basic linear algebra subprograms for Fortran usage,
!          ACM Transactions on Mathematical Software 5, 3 (Sept 1979),
!          pp. 308-323 and 324-325.
!     ------------------------------------------------------------------
!
!     LSQR development:
!     22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583.
!     15 Sep 1985: Final F66 version.  LSQR sent to "misc" in netlib.
!     13 Oct 1987: Bug (Robert Davies, DSIR).  Have to delete
!                     if ( (one + dabs(t)) .le. one ) GO TO 200
!                  from loop 200.  The test was an attempt to reduce
!                  underflows, but caused w(i) not to be updated.
!     17 Mar 1989: First F77 version.
!     04 May 1989: Bug (David Gay, AT&T).  When the second beta is zero,
!                  rnorm = 0 and
!                  test2 = Arnorm / (Anorm * rnorm) overflows.
!                  Fixed by testing for rnorm = 0.
!     05 May 1989: Sent to "misc" in netlib.
!     14 Mar 1990: Bug (John Tomlin via IBM OSL testing).
!                  Setting rhbar2 = rhobar**2 + dampsq can give zero
!                  if rhobar underflows and damp = 0.
!                  Fixed by testing for damp = 0 specially.
!     15 Mar 1990: Converted to lower case.
!     21 Mar 1990: d2norm introduced to avoid overflow in numerous
!                  items like  c = sqrt( a**2 + b**2 ).
!     04 Sep 1991: wantse added as an argument to LSQR, to make
!                  standard errors optional.  This saves storage and
!                  time when se(*) is not wanted.
!     13 Feb 1992: istop now returns a value in [1,5], not [1,7].
!                  1, 2 or 3 means that x solves one of the problems
!                  Ax = b,  min norm(Ax - b)  or  damped least squares.
!                  4 means the limit on cond(A) was reached.
!                  5 means the limit on iterations was reached.
!     07 Dec 1994: Keep track of dxmax = max_k norm( phi_k * d_k ).
!                  So far, this is just printed at the end.
!                  A large value (relative to norm(x)) indicates
!                  significant cancellation in forming
!                  x  =  D*f  =  sum( phi_k * d_k ).
!                  A large column of D need NOT be serious if the
!                  corresponding phi_k is small.
!     27 Dec 1994: Include estimate of alfa_opt in iteration log.
!                  alfa_opt is the optimal scale factor for the
!                  residual in the "augmented system", as described by
!                  A. Bjorck (1992),
!                  Pivoting and stability in the augmented system method,
!                  in D. F. Griffiths and G. A. Watson (eds.),
!                  "Numerical Analysis 1991",
!                  Proceedings of the 14th Dundee Conference,
!                  Pitman Research Notes in Mathematics 260,
!                  Longman Scientific and Technical, Harlow, Essex, 1992.
!
! 07 Sep 2007: Line by line translation for Fortran 90 compilers
!              by Eric Badel <[e-mail]>.
! 16 Sep 2007: Polished a bit by Michael Saunders.
!
! Maintained by
!     Michael A. Saunders
!     Systems Optimization Laboratory (SOL)
!     Stanford University
!     Stanford, CA 94305-4026, USA
!     [e-mail]  [phone]
!-----------------------------------------------------------------------

      external           d2norm, dnrm2, dcopy, dscal
      double precision   d2norm, dnrm2

!     Local variables
      logical            damped, prnt
      integer            i, maxdx, nconv, nstop
      double precision   alfopt,                                   &
                         alpha, beta, bnorm, cs, cs1, cs2, ctol,   &
                         delta, dknorm, dnorm, dxk, dxmax, gamma,  &
                         gambar, phi, phibar, psi, res2, rho,      &
                         rhobar, rhbar1, rhs, rtol, sn, sn1, sn2,  &
                         t, tau, temp, test1, test2, test3, theta, &
                         t1, t2, t3, xnorm1, z, zbar

      double precision, parameter :: zero = 0.0d+0, one = 1.0d+0
      character(len=*), parameter :: enter = ' Enter LSQR.  '
      character(len=*), parameter :: exitt = ' Exit  LSQR.  '
      character(len=*), parameter :: msg(0:5) =                   &
      (/ 'The exact solution is  x = 0                         ', &
         'A solution to Ax = b was found, given atol, btol     ', &
         'A least-squares solution was found, given atol       ', &
         'A damped least-squares solution was found, given atol', &
         'Cond(Abar) seems to be too large, given conlim       ', &
         'The iteration limit was reached                      ' /)
!     ------------------------------------------------------------------

!     Initialize.

      if (nout > 0) then
         write(nout, 1000) enter, m, n, damp, wantse, &
                           atol, conlim, btol, itnlim
      end if

      damped =   damp > zero
      itn    =   0
      istop  =   0
      nstop  =   0
      maxdx  =   0
      ctol   =   zero
      if (conlim > zero) ctol = one / conlim
      Anorm  =   zero
      Acond  =   zero
      dnorm  =   zero
      dxmax  =   zero
      res2   =   zero
      psi    =   zero
      xnorm  =   zero
      xnorm1 =   zero
      cs2    = - one
      sn2    =   zero
      z      =   zero

!     ------------------------------------------------------------------
!     Set up the first vectors u and v for the bidiagonalization.
!     These satisfy  beta*u = b,  alpha*v = A(transpose)*u.
!     ------------------------------------------------------------------
      v(1:n) = zero
      x(1:n) = zero

      if (wantse) then
         se(1:n) = zero
      end if

      alpha  = zero
      beta   = dnrm2 (m, u, 1)

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

      if (alpha > zero) then
         call dscal (n, (one/alpha), v, 1)
         call dcopy (n, v, 1, w, 1)
      end if

      Arnorm = alpha*beta
      if (Arnorm == zero) go to 800

      rhobar = alpha
      phibar = beta
      bnorm  = beta
      rnorm  = beta

      if (nout > 0) then
         if (damped) then
             write(nout,1300)
         else
             write(nout,1200)
         end if
         test1 = one
         test2 = alpha/beta
         write(nout, 1500) itn,x(1),rnorm,test1,test2
      end if


!     ==================================================================
!     Main iteration loop.
!     ==================================================================
  100 itn = itn + 1
	  
!     ------------------------------------------------------------------
!     Perform the next step of the bidiagonalization to obtain the
!     next  beta, u, alpha, v.  These satisfy the relations
!                beta*u  =  A*v  -  alpha*u,
!               alpha*v  =  A(transpose)*u  -  beta*v.
!     ------------------------------------------------------------------
      call dscal (m,(- alpha), u, 1)
      call Aprod (1, m, n, v, u, leniw, lenrw, iw, rw)
      beta   = dnrm2 (m,u,1)

!     Accumulate  Anorm = || Bk ||
!                       =  sqrt( sum of  alpha**2 + beta**2 + damp**2 ).