symmlq.f

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

*     to their single or double equivalents.
*     ------------------------------------------------------------------
*
*
*     This routine is an implementation of the algorithm described in
*     the following references:
*
*     C.C. Paige and M.A. Saunders,  Solution of Sparse Indefinite
*          Systems of Linear Equations,
*          SIAM J. Numer. Anal. 12, 4, September 1975, pp. 617-629.
*
*     J.G. Lewis,  Algorithms for Sparse Matrix Eigenvalue Problems,
*          Report STAN-CS-77-595, Computer Science Department,
*          Stanford University, Stanford, California, March 1977.
*
*     Applications of SYMMLQ and the theory of preconditioning
*     are described in the following references:
*
*     D.B. Szyld and O.B. Widlund,  Applications of Conjugate Gradient
*          Type Methods to Eigenvalue Calculations,
*          in R. Vichnevetsky and R.S. Steplman (editors),
*          Advances in Computer Methods for Partial Differential
*          Equations -- III, IMACS, 1979, 167-173.
*
*     D.B. Szyld,  A Two-level Iterative Method for Large Sparse
*          Generalized Eigenvalue Calculations,
*          Ph. D. dissertation, Department of Mathematics,
*          New York University, New York, October 1983.
*
*     P.E. Gill, W. Murray, D.B. Ponceleon and M.A. Saunders,
*          Preconditioners for indefinite systems arising in
*          optimization, SIMAX 13, 1, 292--311, January 1992.
*          (SIAM J. on Matrix Analysis and Applications)
*     ------------------------------------------------------------------
*
*
*     SYMMLQ development:
*            1972: First version.
*            1975: John Lewis recommended modifications to help with
*                  inverse iteration:
*                  1. Reorthogonalize v1 and v2.
*                  2. Regard the solution as x = x1  +  bstep * b,
*                     with x1 and bstep accumulated separately
*                     and bstep * b added at the end.
*                     (In inverse iteration, b might be close to the
*                     required x already, so x1 may be a lot smaller
*                     than the multiple of b.)
*            1978: Daniel Szyld and Olof Widlund implemented the first
*                  form of preconditioning.
*                  This required both a solve and a multiply with M.
*            1979: Implemented present method for preconditioning.
*                  This requires only a solve with M.
*            1984: Sven Hammarling noted corrections to tnorm and x1lq.
*                  SYMMLQ added to NAG Fortran Library.
*     15 Sep 1985: Final F66 version.  SYMMLQ sent to "misc" in netlib.
*     16 Feb 1989: First F77 version.
*
*     22 Feb 1989: Hans Mittelmann observed beta2 = 0 (hence failure)
*                  if Abar = const*I.  istop = -1 added for this case.
*
*     01 Mar 1989: Hans Mittelmann observed premature termination on
*                  ( 1  1  1 )     (   )                   ( 1  1    )
*                  ( 1  1    ) x = ( 1 ),  for which  T3 = ( 1  1  1 ).
*                  ( 1     1 )     (   )                   (    1  1 )
*                  T2 is exactly singular, so estimating cond(A) from
*                  the diagonals of Lbar is unsafe.  We now use
*                  L       or  Lbar         depending on whether
*                  lqnorm  or  cgnorm       is least.
*
*     03 Mar 1989: eps computed internally instead of coming in as a
*                  parameter.
*     07 Jun 1989: ncheck added as a parameter to say if A and M
*                  should be checked for symmetry.
*                  Later changed to checka (see below).
*     20 Nov 1990: goodb added as a parameter to make Lewis's changes
*                  an option.  Usually b is NOT much like x.  Setting
*                  goodb = .false. saves a call to msolve at the end.
*     20 Nov 1990: Residual not computed exactly at end, to save time
*                  when only one or two iterations are required
*                  (e.g. if the preconditioner is very good).
*                  Beware, if precon is true, rnorm estimates the
*                  residual of the preconditioned system, not Ax = b.
*     04 Sep 1991: Parameter list changed and reordered.
*                  integer ncheck is now logical checka.
*     22 Jul 1992: Example from Lothar Reichel and Daniela Calvetti
*                  showed that beta2 = 0 (istop = -1) means that
*                  b is an eigenvector when M = I.
*                  More complicated if there is a preconditioner;
*                  not clear yet how to describe it.
*    14 Dec 1992:  Modified by Robert Leland, Sandia National Laboratories
*		   to integrate with a C application code. The matrix
*                  data is now passed by reference through symmlq to
*                  aprod and msolve. These are now just Fortran wrappers
*                  for C codes consistent with the matrix data passed
*                  via the pointers "A", "vwsqrt", "work" and "orthlist"
*    10 Feb 1993:  Modified by Robert Leland to return calculate machine 
*		   precision and terminate if the norm of the iterate gets
*		   above the limit normxlim. Relevant for inverse iteration.
*		   Also incorporated itnmin to enforce minimum number itns.
*    17 Aug 1993:  Observed that the Fortran i/o in this routine won't
*		   work because there is no main fortran program to open
*                  the standard i/o files. So for this (and other reasons)
*                  converted the Fortran to C, necessitating inclusion of
*                  the file f2c.h. To avoid a problem with maintaining
*		   Symmlq, I commented out the i/o within it and instead
*		   report its performance based only on the return value
*	   	   of various parameters. That means we can modifiy the
*		   Fortran source, run f2c and recompile without losing or
*		   re-writing any functionality. 
*
*     Michael A. Saunders                    [e-mail]
*     Department of Operations Research    [e-mail]
*     Stanford University
*     Stanford, CA 94305-4022                             [phone]
*     ------------------------------------------------------------------
*
*
*     Subroutines and functions
*
*     USER       aprod, msolve
*     BLAS       daxpy, dcopy, ddot , dnrm2
*
*
*     Intrinsics and local variables

      intrinsic          abs, max, min, mod, sqrt
      double precision   ddot, dnrm2
      double precision   alfa, b1, beta, beta1, bstep, cs,
     $                   cgnorm, dbar, delta, denom, diag,
     $                   eps, epsa, epsln, epsr, epsx,
     $                   gamma, gbar, gmax, gmin, 
     $                   lqnorm, oldb, qrnorm, rhs1, rhs2,
     $                   s, sn, snprod, t, tnorm,
     $                   x1cg, x1lq, ynorm2, zbar, z
      integer            i

      double precision   zero         ,  one         ,  two
      parameter        ( zero = 0.0d+0,  one = 1.0d+0,  two = 2.0d+0 )

      character*16       enter, exit
      character*52       msg(-1:9)

      data               enter /' Enter SYMMLQ.  '/,
     $                   exit  /' Exit  SYMMLQ.  '/

      data               msg
     $ / 'beta2 = 0.  If M = I, b and x are eigenvectors of A',
     $   'beta1 = 0.  The exact solution is  x = 0',
     $   'Requested accuracy achieved, as determined by rtol',
     $   'Reasonable accuracy achieved, given eps',
     $   'x has converged to an eigenvector',
     $   'acond has exceeded 0.1/eps',
     $   'The iteration limit was reached',
     $   'aprod  does not define a symmetric matrix',
     $   'msolve does not define a symmetric matrix',
     $   'msolve does not define a pos-def preconditioner',
     $   'Norm of iterate > max for well conditioned system' /
*     ------------------------------------------------------------------


*     Compute eps, the machine precision.  The call to daxpy is
*     intended to fool compilers that use extra-length registers.

      eps    = one / 16.0d+0

   10 eps    = eps / two
      x(1)   = eps
      y(1)   = one
      call daxpy ( 1, one, x, 1, y, 1 )
      if (y(1) .gt. one) go to 10

      eps    = eps * two

*     Return the calculated machine precision - rwl
      macheps = eps

*     Print heading and initialize.

c   This i/o won't work - see note in preamble..
c     if (nout .gt. 0) then
c        write(nout, 1000) enter, n, checka, goodb, precon,
c     $                     itnlim, rtol, shift
c     end if

      istop  = 0
      itn    = 0
      anorm  = zero
      acond  = zero
      rnorm  = zero
      ynorm  = zero

      do 50 i = 1, n
         x(i) = zero
   50 continue

*     Set up y for the first Lanczos vector v1.
*     y is really beta1 * P * v1  where  P = C**(-1).
*     y and beta1 will be zero if b = 0.

      call dcopy ( n, b, 1, y , 1 )
      call dcopy ( n, b, 1, r1, 1 )
      if ( precon ) call msolve( n, r1, y, A, vwsqrt, work )
c     if ( goodb  ) then
c        b1  = y(1)
c     else
c        b1  = zero
c     end if
      beta1  = ddot  ( n, r1, 1, y, 1 )

*     See if msolve is symmetric.

      if (checka  .and.  precon) then
         call msolve( n, y, r2, A, vwsqrt, work )
         s      = ddot  ( n, y, 1, y, 1 )
         t      = ddot  ( n,r1, 1,r2, 1 )
         z      = abs( s - t )
         epsa   = (s + eps) * eps**0.33333
         if (z .gt. epsa) then
            istop = 7
            go to 900
         end if
      end if

*     Test for an indefinite preconditioner.

      if (beta1 .lt. zero) then
         istop = 8
         go to 900
      end if

*     If b = 0 exactly, stop with x = 0.

      if (beta1 .eq. zero) then
         go to 900
      end if

*     Here and later, v is really P * (the Lanczos v).

      beta1  = sqrt( beta1 )
      s      = one / beta1
      do 100 i = 1, n
         v(i)  = s * y(i)
  100 continue

*     See if aprod  is symmetric.

      call aprod( n, v, y, A, vwsqrt, work, orthlist )
      if (checka) then
         call aprod ( n, y, r2, A, vwsqrt, work, orthlist )
         s      = ddot  ( n, y, 1, y, 1 )
         t      = ddot  ( n, v, 1,r2, 1 )
         z      = abs( s - t )
         epsa   = (s + eps) * eps**0.33333
         if (z .gt. epsa) then
            istop = 6
            go to 900
         end if
      end if

*     Set up y for the second Lanczos vector.
*     Again, y is beta * P * v2  where  P = C**(-1).