lsqr.f90

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

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!     File lsqr.f90  (double precision)
!
!     LSQR     d2norm
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine LSQR                              &
         (m,n,Aprod,damp,wantse,leniw,lenrw,iw,rw, &
          u,v,w,x,se,                              &
          atol,btol,conlim,itnlim,nout,            &
          istop,itn,Anorm,Acond,rnorm,Arnorm,xnorm)

      implicit           none
      external           Aprod
      logical            wantse
      integer            m, n, leniw, lenrw, itnlim, nout, istop, itn
      integer            iw(leniw)
      double precision   rw(lenrw), u(m), v(n), w(n), x(n), se(*)
      double precision   atol, btol, conlim, damp
      double precision   Anorm, Acond, rnorm, Arnorm, xnorm

!-----------------------------------------------------------------------
!
!     LSQR  finds a solution x to the following problems:
!
!     1. Unsymmetric equations --    solve  A*x = b
!
!     2. Linear least squares  --    solve  A*x = b
!                                    in the least-squares sense
!
!     3. Damped least squares  --    solve  (   A    )*x = ( b )
!                                           ( damp*I )     ( 0 )
!                                    in the least-squares sense
!
!     where A is a matrix with m rows and n columns, b is an
!     m-vector, and damp is a scalar.  (All quantities are real.)
!     The matrix A is intended to be large and sparse.  It is accessed
!     by means of subroutine calls of the form
!
!                call Aprod ( mode, m, n, x, y, leniw, lenrw, iw, rw )
!
!     which must perform the following functions:
!
!                If mode = 1, compute  y = y + A*x.
!                If mode = 2, compute  x = x + A(transpose)*y.
!
!     The vectors x and y are input parameters in both cases.
!     If  mode = 1,  y should be altered without changing x.
!     If  mode = 2,  x should be altered without changing y.
!     The parameters leniw, lenrw, iw, rw may be used for workspace
!     as described below.
!
!     The rhs vector b is input via u, and subsequently overwritten.
!
!     Note:  LSQR uses an iterative method to approximate the solution.
!     The number of iterations required to reach a certain accuracy
!     depends strongly on the scaling of the problem.  Poor scaling of
!     the rows or columns of A should therefore be avoided where
!     possible.
!
!     For example, in problem 1 the solution is unaltered by
!     row-scaling.  If a row of A is very small or large compared to
!     the other rows of A, the corresponding row of ( A  b ) should be
!     scaled up or down.
!
!     In problems 1 and 2, the solution x is easily recovered
!     following column-scaling.  Unless better information is known,
!     the nonzero columns of A should be scaled so that they all have
!     the same Euclidean norm (e.g., 1.0).
!
!     In problem 3, there is no freedom to re-scale if damp is
!     nonzero.  However, the value of damp should be assigned only
!     after attention has been paid to the scaling of A.
!
!     The parameter damp is intended to help regularize
!     ill-conditioned systems, by preventing the true solution from
!     being very large.  Another aid to regularization is provided by
!     the parameter Acond, which may be used to terminate iterations
!     before the computed solution becomes very large.
!
!     Note that x is not an input parameter.
!     If some initial estimate x0 is known and if damp = 0,
!     one could proceed as follows:
!
!       1. Compute a residual vector     r0 = b - A*x0.
!       2. Use LSQR to solve the system  A*dx = r0.
!       3. Add the correction dx to obtain a final solution x = x0 + dx.
!
!     This requires that x0 be available before and after the call
!     to LSQR.  To judge the benefits, suppose LSQR takes k1 iterations
!     to solve A*x = b and k2 iterations to solve A*dx = r0.
!     If x0 is "good", norm(r0) will be smaller than norm(b).
!     If the same stopping tolerances atol and btol are used for each
!     system, k1 and k2 will be similar, but the final solution x0 + dx
!     should be more accurate.  The only way to reduce the total work
!     is to use a larger stopping tolerance for the second system.
!     If some value btol is suitable for A*x = b, the larger value
!     btol*norm(b)/norm(r0)  should be suitable for A*dx = r0.
!
!     Preconditioning is another way to reduce the number of iterations.
!     If it is possible to solve a related system M*x = b efficiently,
!     where M approximates A in some helpful way
!     (e.g. M - A has low rank or its elements are small relative to
!     those of A), LSQR may converge more rapidly on the system
!           A*M(inverse)*z = b,
!     after which x can be recovered by solving M*x = z.
!
!     NOTE: If A is symmetric, LSQR should not be used!
!     Alternatives are the symmetric conjugate-gradient method (cg)
!     and/or SYMMLQ.
!     SYMMLQ is an implementation of symmetric cg that applies to
!     any symmetric A and will converge more rapidly than LSQR.
!     If A is positive definite, there are other implementations of
!     symmetric cg that require slightly less work per iteration
!     than SYMMLQ (but will take the same number of iterations).
!
!
!     Notation
!     --------
!     The following quantities are used in discussing the subroutine
!     parameters:
!
!     Abar   =  (   A    ),          bbar  =  ( b )
!               ( damp*I )                    ( 0 )
!
!     r      =  b  -  A*x,           rbar  =  bbar  -  Abar*x
!
!     rnorm  =  sqrt( norm(r)**2  +  damp**2 * norm(x)**2 )
!            =  norm( rbar )
!
!     relpr  =  the relative precision of floating-point arithmetic
!               on the machine being used.  On most machines,
!               relpr is about 1.0e-7 and 1.0d-16 in single and double
!               precision respectively.
!
!     LSQR  minimizes the function rnorm with respect to x.
!
!
!     Parameters
!     ----------
!     m       input      m, the number of rows in A.
!
!     n       input      n, the number of columns in A.
!
!     Aprod   external   See above.
!
!     damp    input      The damping parameter for problem 3 above.
!                        (damp should be 0.0 for problems 1 and 2.)
!                        If the system A*x = b is incompatible, values
!                        of damp in the range 0 to sqrt(relpr)*norm(A)
!                        will probably have a negligible effect.
!                        Larger values of damp will tend to decrease
!                        the norm of x and reduce the number of 
!                        iterations required by LSQR.
!
!                        The work per iteration and the storage needed
!                        by LSQR are the same for all values of damp.
!
!     wantse  input      A logical variable to say if the array se(*)
!                        of standard error estimates should be computed.
!                        If m > n  or  damp > 0,  the system is
!                        overdetermined and the standard errors may be
!                        useful.  (See the first LSQR reference.)
!                        Otherwise (m <= n  and  damp = 0) they do not
!                        mean much.  Some time and storage can be saved
!                        by setting wantse = .false. and using any
!                        convenient array for se(*), which won't be
!                        touched.
!
!     leniw   input      The length of the workspace array iw.
!     lenrw   input      The length of the workspace array rw.
!     iw      workspace  An integer array of length leniw.
!     rw      workspace  A real array of length lenrw.
!
!             Note:  LSQR  does not explicitly use the previous four
!             parameters, but passes them to subroutine Aprod for
!             possible use as workspace.  If Aprod does not need
!             iw or rw, the values leniw = 1 or lenrw = 1 should
!             be used, and the actual parameters corresponding to
!             iw or rw  may be any convenient array of suitable type.
!
!     u(m)    input      The rhs vector b.  Beware that u is
!                        over-written by LSQR.
!
!     v(n)    workspace
!
!     w(n)    workspace
!
!     x(n)    output     Returns the computed solution x.
!
!     se(*)   output     If wantse is true, the dimension of se must be
!             (maybe)    n or more.  se(*) then returns standard error
!                        estimates for the components of x.
!                        For each i, se(i) is set to the value
!                           rnorm * sqrt( sigma(i,i) / t ),
!                        where sigma(i,i) is an estimate of the i-th
!                        diagonal of the inverse of Abar(transpose)*Abar
!                        and  t = 1      if  m .le. n,
!                             t = m - n  if  m .gt. n  and  damp = 0,
!                             t = m      if  damp .ne. 0.
!
!                        If wantse is false, se(*) will not be touched.
!                        The actual parameter can be any suitable array
!                        of any length.
!
!     atol    input      An estimate of the relative error in the data
!                        defining the matrix A.  For example,
!                        if A is accurate to about 6 digits, set
!                        atol = 1.0e-6 .
!
!     btol    input      An estimate of the relative error in the data
!                        defining the rhs vector b.  For example,
!                        if b is accurate to about 6 digits, set
!                        btol = 1.0e-6 .
!
!     conlim  input      An upper limit on cond(Abar), the apparent
!                        condition number of the matrix Abar.
!                        Iterations will be terminated if a computed
!                        estimate of cond(Abar) exceeds conlim.
!                        This is intended to prevent certain small or
!                        zero singular values of A or Abar from
!                        coming into effect and causing unwanted growth
!                        in the computed solution.
!
!                        conlim and damp may be used separately or
!                        together to regularize ill-conditioned systems.
!
!                        Normally, conlim should be in the range
!                        1000 to 1/relpr.
!                        Suggested value:
!                        conlim = 1/(100*relpr)  for compatible systems,
!                        conlim = 1/(10*sqrt(relpr)) for least squares.
!
!             Note:  If the user is not concerned about the parameters
!             atol, btol and conlim, any or all of them may be set
!             to zero.  The effect will be the same as the values
!             relpr, relpr and 1/relpr respectively.
!
!     itnlim  input      An upper limit on the number of iterations.
!                        Suggested value:
!                        itnlim = n/2   for well-conditioned systems
!                                       with clustered singular values,
!                        itnlim = 4*n   otherwise.
!
!     nout    input      File number for printed output.  If positive,
!                        a summary will be printed on file nout.
!
!     istop   output     An integer giving the reason for termination:
!
!                0       x = 0  is the exact solution.
!                        No iterations were performed.
!
!                1       The equations A*x = b are probably
!                        compatible.  Norm(A*x - b) is sufficiently
!                        small, given the values of atol and btol.
!
!                2       damp is zero.  The system A*x = b is probably