symmlq.cpp

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

/* 		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 */
/*     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 */
    /* Parameter adjustments */
    --vwsqrt;
    --work;
    --y;
    --x;
    --w;
    --v;
    --r2;
    --r1;
    --b;

    /* Function Body */
/*     ------------------------------------------------------------------ 
*/
/*     Compute eps, the machine precision.  The call to daxpy is */
/*     intended to fool compilers that use extra-length registers. */
    eps = .0625;
L10:
    eps /= 2.;
    x[1] = eps;
    y[1] = 1.;
    daxpy_(&c__1, &c_b4, &x[1], &c__1, &y[1], &c__1);
    if (y[1] > 1.) {
	goto L10;
    }
    eps *= 2.;
/*     Return the calculated machine precision - rwl */
    *macheps = eps;
/*     Print heading and initialize. */
/*   This i/o won't work - see note in preamble.. */
/*     if (nout .gt. 0) then */
/*        write(nout, 1000) enter, n, checka, goodb, precon, */
/*     $                     itnlim, rtol, shift */
/*     end if */
    *istop = 0;
    *itn = 0;
    *anorm = 0.;
    *acond = 0.;
    *rnorm = 0.;
    *ynorm = 0.;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	x[i] = 0.;
/* L50: */
    }
/*     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. */
    dcopy_(n, &b[1], &c__1, &y[1], &c__1);
    dcopy_(n, &b[1], &c__1, &r1[1], &c__1);
    if (*precon) {
	msolve_(n, &r1[1], &y[1], a, &vwsqrt[1], &work[1]);
    }
/*     if ( goodb  ) then */
/*        b1  = y(1) */
/*     else */
/*        b1  = zero */
/*     end if */
    beta1 = ddot_(n, &r1[1], &c__1, &y[1], &c__1);
/*     See if msolve is symmetric. */
    if (*checka && *precon) {
	msolve_(n, &y[1], &r2[1], a, &vwsqrt[1], &work[1]);
	s = ddot_(n, &y[1], &c__1, &y[1], &c__1);
	t = ddot_(n, &r1[1], &c__1, &r2[1], &c__1);
	z = (d__1 = s - t, abs(d__1));
	epsa = (s + eps) * pow_dd(&eps, &c_b18);
	if (z > epsa) {
	    *istop = 7;
	    goto L900;
	}
    }
/*     Test for an indefinite preconditioner. */
    if (beta1 < 0.) {
	*istop = 8;
	goto L900;
    }
/*     If b = 0 exactly, stop with x = 0. */
    if (beta1 == 0.) {
	goto L900;
    }
/*     Here and later, v is really P * (the Lanczos v). */
    beta1 = sqrt(beta1);
    s = 1. / beta1;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	v[i] = s * y[i];
/* L100: */
    }
/*     See if aprod  is symmetric. */
    aprod_(n, &v[1], &y[1], a, &vwsqrt[1], &work[1], orthlist);
    if (*checka) {
	aprod_(n, &y[1], &r2[1], a, &vwsqrt[1], &work[1], orthlist);
	s = ddot_(n, &y[1], &c__1, &y[1], &c__1);
	t = ddot_(n, &v[1], &c__1, &r2[1], &c__1);
	z = (d__1 = s - t, abs(d__1));
	epsa = (s + eps) * pow_dd(&eps, &c_b18);
	if (z > epsa) {
	    *istop = 6;
	    goto L900;
	}
    }
/*     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. */
    d__1 = -(*shift);
    daxpy_(n, &d__1, &v[1], &c__1, &y[1], &c__1);