*** empty log message ***

Changed files:
    symmlq.f -> 1.1
    symmlq_f77.README -> 1.1

2 files changed, 750 insertions, 0 deletions

diff --git a/symmlq.f b/symmlq.f
new file mode 100644
index 0000000..e4ae134
--- /dev/null
+++ b/symmlq.f
@@ -0,0 +1,732 @@
+*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+*     File symmlq.f
+*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+      subroutine SYMMLQ( n, b, r1, r2, v, w, x, y,
+     $                   Aprod, Msolve, checkA, goodb, precon, shift,
+     $                   nout , itnlim, rtol,
+     $                   istop, itn, Anorm, Acond, rnorm, ynorm )
+
+      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)
+*     ------------------------------------------------------------------
+*
+*     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
+*
+*                call Aprod ( n, x, y )
+*
+*     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
+*
+*                              call Aprod ( n, x, y )
+*
+*                        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
+*
+*                              call Msolve( n, x, y )
+*
+*                        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.
+*
+*     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).
+*
+*
+*
+*     To change precision
+*     -------------------
+*
+*     Alter the words
+*            double precision,
+*            daxpy, dcopy, ddot, dnrm2
+*     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.
+*     20 Oct 1999: Bug.  alfa1 = 0 caused Anorm = 0, divide by zero.
+*                  Need to estimate Anorm from column of Tk.
+*
+*     Michael A. Saunders                   na.msaunders@na-net.ornl.gov
+*     Department of EESOR                  mike@sol-michael.stanford.edu
+*     Stanford University
+*     Stanford, CA 94305-4023                             (650) 723-1875
+*     ------------------------------------------------------------------
+*
+*
+*     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, gpert,
+     $                   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:8)
+
+      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' /
+*     ------------------------------------------------------------------
+
+
+*     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
+
+*     Print heading and initialize.
+
+      if (nout .gt. 0) then
+         write(nout, 1000) enter, n, checkA, goodb, precon,
+     $                     itnlim, rtol, shift
+      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 )
+      if ( goodb  ) then
+         b1  = y(1)
+      else
+         b1  = zero
+      end if
+      beta1  = ddot  ( n, r1, 1, y, 1 )
+
+*     See if Msolve is symmetric.
+
+      if (checkA  .and.  precon) then
+         call Msolve( n, y, r2 )
+         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 )
+      if (checkA) then
+         call Aprod ( n, y, r2 )
+         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).
+*     y and beta will be zero or very small if b is an eigenvector.
+
+      call daxpy ( n, (- shift), v, 1, y, 1 )
+      alfa   = ddot  ( n, v, 1, y, 1 )
+      call daxpy ( n, (- alfa / beta1), r1, 1, y, 1 )
+
+*     Make sure  r2  will be orthogonal to the first  v.
+
+      z      = ddot  ( n, v, 1, y, 1 )
+      s      = ddot  ( n, v, 1, v, 1 )
+      call daxpy ( n, (- z / s), v, 1, y, 1 )
+
+      call dcopy ( n, y, 1, r2, 1 )
+      if ( precon ) call Msolve( n, r2, y )
+      oldb   = beta1
+      beta   = ddot  ( n, r2, 1, y, 1 )
+
+      if (beta .lt. zero) then
+         istop = 8
+         go to 900
+      end if
+
+*     Cause termination (later) if beta is essentially zero.
+
+      beta   = sqrt( beta )
+      if (beta .le. eps) then
+         istop = -1
+      end if
+
+*     See if the local reorthogonalization achieved anything.
+
+      denom  = sqrt( s ) * dnrm2( n, r2, 1 )  +  eps
+      s      = z / denom
+      t      = ddot  ( n, v, 1, r2, 1 ) / denom
+      if (nout .gt. 0  .and.  goodb) then
+         write(nout, 1100) beta1, alfa, s, t
+      end if
+
+*     Initialize other quantities.
+
+      cgnorm = beta1
+      gbar   = alfa
+      dbar   = beta
+      rhs1   = beta1
+      rhs2   = zero
+      bstep  = zero
+      snprod = one
+      tnorm  = alfa**2 + beta**2
+      ynorm2 = zero
+      gmax   = abs( alfa ) + eps
+      gmin   = gmax
+
+      if ( goodb ) then
+         do 200 i = 1, n
+            w(i)  = zero
+  200    continue
+      else
+         call dcopy ( n, v, 1, w, 1 )
+      end if
+
+*     ------------------------------------------------------------------
+*     Main iteration loop.
+*     ------------------------------------------------------------------
+
+*     Estimate various norms and test for convergence.
+
+  300 Anorm  = sqrt( tnorm  )
+      ynorm  = sqrt( ynorm2 )
+      epsa   = Anorm * eps
+      epsx   = Anorm * ynorm * eps
+      epsr   = Anorm * ynorm * rtol
+      diag   = gbar
+      if (diag .eq. zero) diag = epsa
+
+      lqnorm = sqrt( rhs1**2 + rhs2**2 )
+      qrnorm = snprod * beta1
+      cgnorm = qrnorm * beta / abs( diag )
+
+*     Estimate  cond(A).
+*     In this version we look at the diagonals of  L  in the
+*     factorization of the tridiagonal matrix,  T = L*Q.
+*     Sometimes, T(k) can be misleadingly ill-conditioned when
+*     T(k+1) is not, so we must be careful not to overestimate Acond.
+
+      if (lqnorm .le. cgnorm) then
+         Acond  = gmax / gmin
+      else
+         denom  = min( gmin, abs( diag ) )
+         Acond  = gmax / denom
+      end if
+
+*     See if any of the stopping criteria are satisfied.
+*     In rare cases, istop is already -1 from above (Abar = const * I).
+
+      if (istop .eq. 0) then
+         if (itn    .ge. itnlim ) istop = 5
+         if (Acond  .ge. 0.1/eps) istop = 4
+         if (epsx   .ge. beta1  ) istop = 3
+         if (cgnorm .le. epsx   ) istop = 2
+         if (cgnorm .le. epsr   ) istop = 1
+      end if
+*     ==================================================================
+
+*     See if it is time to print something.
+
+      if (nout .le.  0)          go to 600
+      if (n    .le. 40)          go to 400
+      if (itn  .le. 10)          go to 400
+      if (itn  .ge. itnlim - 10) go to 400
+      if (mod(itn,10)  .eq.   0) go to 400
+      if (cgnorm .le. 10.0*epsx) go to 400
+      if (cgnorm .le. 10.0*epsr) go to 400
+      if (Acond  .ge. 0.01/eps ) go to 400
+      if (istop  .ne. 0)         go to 400
+      go to 600
+
+*     Print a line for this iteration.
+
+  400 zbar   = rhs1 / diag
+      z      = (snprod * zbar  +  bstep) / beta1
+      x1lq   = x(1)  +  b1 * bstep / beta1
+      x1cg   = x(1)  +  w(1) * zbar  +  b1 * z
+
+      if (    itn     .eq. 0) write(nout, 1200)
+      write(nout, 1300) itn, x1cg, cgnorm, bstep/beta1, Anorm, Acond
+      if (mod(itn,10) .eq. 0) write(nout, 1500)
+*     ==================================================================
+
+
+*     Obtain the current Lanczos vector  v = (1 / beta)*y
+*     and set up  y  for the next iteration.
+
+  600 if (istop .ne. 0) go to 800
+      s      = one / beta
+
+      do 620 i = 1, n
+         v(i)  = s * y(i)
+  620 continue
+
+      call Aprod ( n, v, y )
+      call daxpy ( n, (- shift), v, 1, y, 1 )
+      call daxpy ( n, (- beta / oldb), r1, 1, y, 1 )
+      alfa   = ddot( n, v, 1, y, 1 )
+      call daxpy ( n, (- alfa / beta), r2, 1, y, 1 )
+      call dcopy ( n, r2, 1, r1, 1 )
+      call dcopy ( n, y, 1, r2, 1 )
+      if ( precon ) call Msolve( n, r2, y )
+      oldb   = beta
+      beta   = ddot  ( n, r2, 1, y, 1 )
+      if (beta .lt. zero) then
+         istop = 6
+         go to 800
+      end if
+      beta   = sqrt( beta )
+      tnorm  = tnorm  +  alfa**2  +  oldb**2  +  beta**2
+
+*     Compute the next plane rotation for  Q.
+
+      gamma  = sqrt( gbar**2 + oldb**2 )
+      cs     = gbar / gamma
+      sn     = oldb / gamma
+      delta  = cs * dbar  +  sn * alfa
+      gbar   = sn * dbar  -  cs * alfa
+      epsln  = sn * beta
+      dbar   =            -  cs * beta
+
+*     Update  x.
+
+      z      = rhs1 / gamma
+      s      = z * cs
+      t      = z * sn
+
+      do 700 i = 1, n
+         x(i)  = (w(i) * s   +   v(i) * t)  +  x(i)
+         w(i)  =  w(i) * sn  -   v(i) * cs
+  700 continue
+
+*     Accumulate the step along the direction  b,
+*     and go round again.
+
+      bstep  = snprod * cs * z  +  bstep
+      snprod = snprod * sn
+      gmax   = max( gmax, gamma )
+      gmin   = min( gmin, gamma )
+      ynorm2 = z**2  +  ynorm2
+      rhs1   = rhs2  -  delta * z
+      rhs2   =       -  epsln * z
+      itn    = itn   +  1
+      go to 300
+
+*     ------------------------------------------------------------------
+*     End of main iteration loop.
+*     ------------------------------------------------------------------
+
+*     Move to the CG point if it seems better.
+*     In this version of SYMMLQ, the convergence tests involve
+*     only cgnorm, so we're unlikely to stop at an LQ point,
+*     EXCEPT if the iteration limit interferes.
+
+  800 if (cgnorm .le. lqnorm) then
+         zbar   = rhs1 / diag
+         bstep  = snprod * zbar  +  bstep
+         ynorm  = sqrt( ynorm2  +  zbar**2 )
+         rnorm  = cgnorm
+         call daxpy ( n, zbar, w, 1, x, 1 )
+      else
+         rnorm  = lqnorm
+      end if
+
+      if ( goodb ) then
+
+*        Add the step along  b.
+
+         bstep  = bstep / beta1
+         call dcopy ( n, b, 1, y, 1 )
+         if ( precon ) call Msolve( n, b, y )
+         call daxpy ( n, bstep, y, 1, x, 1 )
+      end if
+
+*     ==================================================================
+*     Display final status.
+*     ==================================================================
+  900 if (nout  .gt. 0) then
+         write(nout, 2000) exit, istop, itn,
+     $                     exit, Anorm, Acond,
+     $                     exit, rnorm, ynorm
+         write(nout, 3000) exit, msg(istop)
+      end if
+
+      return
+
+*     ------------------------------------------------------------------
+ 1000 format(// 1p,    a, 5x, 'Solution of symmetric   Ax = b'
+     $       / ' n      =', i7, 5x, 'checkA =', l4, 12x,
+     $          'goodb  =', l4, 7x, 'precon =', l4
+     $       / ' itnlim =', i7, 5x, 'rtol   =', e11.2, 5x,
+     $          'shift  =', e23.14)
+ 1100 format(/ 1p, ' beta1  =', e10.2, 3x, 'alpha1 =', e10.2
+     $       / ' (v1,v2) before and after ', e14.2
+     $       / ' local reorthogonalization', e14.2)
+ 1200 format(// 5x, 'itn', 7x, 'x1(cg)', 10x,
+     $         'norm(r)', 5x, 'bstep', 7x, 'norm(A)', 3X, 'cond(A)')
+ 1300 format(1p, i8, e19.10, e11.2, e14.5, 2e10.2)
+ 1500 format(1x)
+ 2000 format(/ 1p, a, 6x, 'istop =', i3,   15x, 'itn   =', i8
+     $       /     a, 6x, 'Anorm =', e12.4, 6x, 'Acond =', e12.4
+     $       /     a, 6x, 'rnorm =', e12.4, 6x, 'ynorm =', e12.4)
+ 3000 format(      a, 6x, a )
+*     ------------------------------------------------------------------
+*     end of SYMMLQ
+      end
diff --git a/symmlq_f77.README b/symmlq_f77.README
new file mode 100644
index 0000000..a02b2b3
--- /dev/null
+++ b/symmlq_f77.README
@@ -0,0 +1,18 @@
+symmlq_f77.README
+
+The software for SYMMLQ (f77 version) is provided by SOL, Stanford University
+under the terms of the Common Public License (CPL):
+http://oss.software.ibm.com/developerworks/opensource/license-cpl.html
+
+
+11 Feb 2000: First set of files available for download from SOL.
+
+Please send comments to Michael Saunders, SOL, Stanford University
+                        saunders@stanford.edu  650-723-1875
+-----------------------------------------------------------------------------
+
+The f77 version of SYMMLQ involves the following files:
+
+   symmlq.f
+   symmlqblas.f   (not needed if you have BLAS-1)
+   symmlqtest.f