symmlq.f

一个画有向图的程序。里面含有力导引画图算法等多个经典算法。

      subroutine symmlq( n, b, r1, r2, v, w, x, y, work,
     $                   checka, goodb, precon, shift,
     $                   nout , itnlim, rtol, istop, itn,
     $			 anorm, acond, rnorm, ynorm, A, vwsqrt,
     $                   orthlist, macheps,normxlim,itnmin)

      external           aprod, msolve
      integer            n, nout, itnlim, istop, itn
      logical            checka, goodb, precon
      double precision   shift, rtol, anorm, acond, rnorm, ynorm,
     $                   b(n), r1(n), r2(n), v(n), w(n), x(n), y(n)
      double precision   vwsqrt(n), work(n)
      double precision	 A, orthlist
      integer 		 itnmin 
      double precision   macheps, normxlim
*     ------------------------------------------------------------------
*
*     SYMMLQ  is designed to solve the system of linear equations
*
*                Ax = b
*
*     where A is an n by n symmetric matrix and b is a given vector.
*     The matrix A is not required to be positive definite.
*     (If A is known to be definite, the method of conjugate gradients
*     might be preferred, since it will require about the same number of
*     iterations as SYMMLQ but slightly less work per iteration.)
*
*
*     The matrix A is intended to be large and sparse.  It is accessed
*     by means of a subroutine call of the form
*
*     old        call aprod ( n, x, y )
*     new:       call aprod ( n, x, y, A, vwsqrt, work, orthlist ) -rwl
*
*     which must return the product y = Ax for any given vector x.
*
*
*     More generally, SYMMLQ is designed to solve the system
*
*                (A - shift*I) x = b
*
*     where  shift  is a specified scalar value.  If  shift  and  b
*     are suitably chosen, the computed vector x may approximate an
*     (unnormalized) eigenvector of A, as in the methods of
*     inverse iteration and/or Rayleigh-quotient iteration.
*     Again, the matrix (A - shift*I) need not be positive definite.
*     The work per iteration is very slightly less if  shift = 0.
*
*
*     A further option is that of preconditioning, which may reduce
*     the number of iterations required.  If M = C C' is a positive
*     definite matrix that is known to approximate  (A - shift*I)
*     in some sense, and if systems of the form  My = x  can be
*     solved efficiently, the parameters precon and msolve may be
*     used (see below).  When  precon = .true., SYMMLQ will
*     implicitly solve the system of equations
*
*             P (A - shift*I) P' xbar  =  P b,
*
*     i.e.                  Abar xbar  =  bbar
*     where                         P  =  C**(-1),
*                                Abar  =  P (A - shift*I) P',
*                                bbar  =  P b,
*
*     and return the solution       x  =  P' xbar.
*     The associated residual is rbar  =  bbar - Abar xbar
*                                      =  P (b - (A - shift*I)x)
*                                      =  P r.
*
*     In the discussion below, eps refers to the machine precision.
*     eps is computed by SYMMLQ.  A typical value is eps = 2.22e-16
*     for IBM mainframes and IEEE double-precision arithmetic.
*
*     Parameters
*     ----------
*
*     n       input      The dimension of the matrix A.
*
*     b(n)    input      The rhs vector b.
*
*     r1(n)   workspace
*     r2(n)   workspace
*     v(n)    workspace
*     w(n)    workspace
*
*     x(n)    output     Returns the computed solution  x.
*
*     y(n)    workspace
*
*     aprod   external   A subroutine defining the matrix A.
*                        For a given vector x, the statement
*
*                        old: call aprod ( n, x, y, )
*                        new: call aprod ( n, x, y, A, vwsqrt, work, orthlist ) -rwl
*
*                        must return the product y = Ax
*                        without altering the vector x.
*
*     msolve  external   An optional subroutine defining a
*                        preconditioning matrix M, which should
*                        approximate (A - shift*I) in some sense.
*                        M must be positive definite.
*                        For a given vector x, the statement
*
*                        old:      call msolve( n, x, y )
*                        new:      call msolve( n, x, y )  -rwl
*
*                        must solve the linear system My = x
*                        without altering the vector x.
*
*                        In general, M should be chosen so that Abar has
*                        clustered eigenvalues.  For example,
*                        if A is positive definite, Abar would ideally
*                        be close to a multiple of I.
*                        If A or A - shift*I is indefinite, Abar might
*                        be close to a multiple of diag( I  -I ).
*
*                        NOTE.  The program calling SYMMLQ must declare
*                        aprod and msolve to be external.
*
*     checka  input      If checka = .true., an extra call of aprod will
*                        be used to check if A is symmetric.  Also,
*                        if precon = .true., an extra call of msolve
*                        will be used to check if M is symmetric.
*
*     goodb   input      Usually, goodb should be .false.
*                        If x is expected to contain a large multiple of
*                        b (as in Rayleigh-quotient iteration),
*                        better precision may result if goodb = .true.
*                        See Lewis (1977) below.
*                        When goodb = .true., an extra call to msolve
*                        is required.
*
*     precon  input      If precon = .true., preconditioning will
*                        be invoked.  Otherwise, subroutine msolve
*                        will not be referenced; in this case the
*                        actual parameter corresponding to msolve may
*                        be the same as that corresponding to aprod.
*
*     shift   input      Should be zero if the system Ax = b is to be
*                        solved.  Otherwise, it could be an
*                        approximation to an eigenvalue of A, such as
*                        the Rayleigh quotient b'Ab / (b'b)
*                        corresponding to the vector b.
*                        If b is sufficiently like an eigenvector
*                        corresponding to an eigenvalue near shift,
*                        then the computed x may have very large
*                        components.  When normalized, x may be
*                        closer to an eigenvector than b.
*
*     nout    input      A file number.
*                        If nout .gt. 0, a summary of the iterations
*                        will be printed on unit nout.
*
*     itnlim  input      An upper limit on the number of iterations.
*
*     rtol    input      A user-specified tolerance.  SYMMLQ terminates
*                        if it appears that norm(rbar) is smaller than
*                              rtol * norm(Abar) * norm(xbar),
*                        where rbar is the transformed residual vector,
*                              rbar = bbar - Abar xbar.
*
*                        If shift = 0 and precon = .false., SYMMLQ
*                        terminates if norm(b - A*x) is smaller than
*                              rtol * norm(A) * norm(x).
*
*     istop   output     An integer giving the reason for termination...
*
*              -1        beta2 = 0 in the Lanczos iteration; i.e. the
*                        second Lanczos vector is zero.  This means the
*                        rhs is very special.
*                        If there is no preconditioner, b is an
*                        eigenvector of A.
*                        Otherwise (if precon is true), let My = b.
*                        If shift is zero, y is a solution of the
*                        generalized eigenvalue problem Ay = lambda My,
*                        with lambda = alpha1 from the Lanczos vectors.
*
*                        In general, (A - shift*I)x = b
*                        has the solution         x = (1/alpha1) y
*                        where My = b.
*                        
*               0        b = 0, so the exact solution is x = 0.
*                        No iterations were performed.
*
*               1        Norm(rbar) appears to be less than
*                        the value  rtol * norm(Abar) * norm(xbar).
*                        The solution in  x  should be acceptable.
*
*               2        Norm(rbar) appears to be less than
*                        the value  eps * norm(Abar) * norm(xbar).
*                        This means that the residual is as small as
*                        seems reasonable on this machine.
*
*               3        Norm(Abar) * norm(xbar) exceeds norm(b)/eps,
*                        which should indicate that x has essentially
*                        converged to an eigenvector of A
*                        corresponding to the eigenvalue shift.
*
*               4        acond (see below) has exceeded 0.1/eps, so
*                        the matrix Abar must be very ill-conditioned.
*                        x may not contain an acceptable solution.
*
*               5        The iteration limit was reached before any of
*                        the previous criteria were satisfied.
*
*               6        The matrix defined by aprod does not appear
*                        to be symmetric.
*                        For certain vectors y = Av and r = Ay, the
*                        products y'y and r'v differ significantly.
*
*               7        The matrix defined by msolve does not appear
*                        to be symmetric.
*                        For vectors satisfying My = v and Mr = y, the
*                        products y'y and r'v differ significantly.
*
*               8        An inner product of the form  x' M**(-1) x
*                        was not positive, so the preconditioning matrix
*                        M does not appear to be positive definite.
*
*                        If istop .ge. 5, the final x may not be an
*                        acceptable solution.
*
*		9	 The norm of the iterate is > than normxlim. 
*			 Termination enforced on presumption that inverse
*			 iteration is being performed -rwl.
*
*     itn     output     The number of iterations performed.
*
*     anorm   output     An estimate of the norm of the matrix operator
*                        Abar = P (A - shift*I) P',   where P = C**(-1).
*
*     acond   output     An estimate of the condition of Abar above.
*                        This will usually be a substantial
*                        under-estimate of the true condition.
*
*     rnorm   output     An estimate of the norm of the final
*                        transformed residual vector,
*                           P (b  -  (A - shift*I) x).
*
*     ynorm   output     An estimate of the norm of xbar.
*                        This is sqrt( x'Mx ).  If precon is false,
*                        ynorm is an estimate of norm(x).
*
*     A	      input      A pointer variable to the matrix data. Passed 
*	                 in to use in revised call to aprod and msolve
*			 added 14 Dec 92 by rwl.
*
*    macheps  output     Used to return the calculated machine precision.
*                        Added 10 Feb 93 by rwl.
*
*    normxlim input	Used as possible termination criterion. 10 Feb 93 rwl.
*
*    itnmin input	Used to enforce a minimum number of itns. 10 Feb 93 rwl.
*
*     To change precision
*     -------------------
*
*     Alter the words
*            double precision,
*            daxpy, dcopy, ddot, dnrm2