From 19809d85fcbb926a8a947572db7c0a0203d38cc9 Mon Sep 17 00:00:00 2001 From: Roman Werpachowski Date: Thu, 16 Dec 2004 18:39:42 +0000 Subject: *** empty log message *** Changed files: symmlq.f -> 1.1 symmlq_f77.README -> 1.1 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 -- cgit v0.10.2