symmlq.cpp

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

    alfa = ddot_(n, &v[1], &c__1, &y[1], &c__1);
    d__1 = -alfa / beta1;
    daxpy_(n, &d__1, &r1[1], &c__1, &y[1], &c__1);
/*     Make sure  r2  will be orthogonal to the first  v. */
    z = ddot_(n, &v[1], &c__1, &y[1], &c__1);
    s = ddot_(n, &v[1], &c__1, &v[1], &c__1);
    d__1 = -z / s;
    daxpy_(n, &d__1, &v[1], &c__1, &y[1], &c__1);
    dcopy_(n, &y[1], &c__1, &r2[1], &c__1);
    if (*precon) {
	msolve_(n, &r2[1], &y[1], a, &vwsqrt[1], &work[1]);
    }
    oldb = beta1;
    beta = ddot_(n, &r2[1], &c__1, &y[1], &c__1);
    if (beta < 0.) {
	*istop = 8;
	goto L900;
    }
/*     Cause termination (later) if beta is essentially zero. */
    beta = sqrt(beta);
    if (beta <= eps) {
	*istop = -1;
    }
/*     See if the local reorthogonalization achieved anything. */
    denom = sqrt(s) * dnrm2_(n, &r2[1], &c__1) + eps;
    s = z / denom;
    t = ddot_(n, &v[1], &c__1, &r2[1], &c__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 = 0.;
    bstep = 0.;
    snprod = 1.;
/* Computing 2nd power */
    d__1 = alfa;
    tnorm = d__1 * d__1;
    ynorm2 = 0.;
    gmax = abs(alfa);
    gmin = gmax;
    if (*goodb) {
	i__1 = *n;
	for (i = 1; i <= i__1; ++i) {
	    w[i] = 0.;
/* L200: */
	}
    } else {
	dcopy_(n, &v[1], &c__1, &w[1], &c__1);
    }
/*     ------------------------------------------------------------------ 
*/
/*     Main iteration loop. */
/*     ------------------------------------------------------------------ 
*/
/*     Estimate various norms and test for convergence. */
L300:
    *anorm = sqrt(tnorm);
    *ynorm = sqrt(ynorm2);
    epsa = *anorm * eps;
    epsx = *anorm * *ynorm * eps;
    epsr = *anorm * *ynorm * *rtol;
    diag = gbar;
    if (diag == 0.) {
	diag = epsa;
    }
/* Computing 2nd power */
    d__1 = rhs1;
/* Computing 2nd power */
    d__2 = rhs2;
    lqnorm = sqrt(d__1 * d__1 + d__2 * d__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 <= cgnorm) {
	*acond = gmax / gmin;
    } else {
/* Computing MIN */
	d__1 = gmin, d__2 = abs(diag);
	denom = min(d__1,d__2);
	*acond = gmax / denom;
    }
/*     See if any of the stopping criteria are satisfied. */
/*     In rare cases, istop is already -1 from above (Abar = const * I). 
*/
    if (*istop == 0) {
	if (*itn >= *itnlim) {
	    *istop = 5;
	}
/*         if (acond  .ge. 0.1/eps) istop = 4 */
	if (epsx >= beta1) {
	    *istop = 3;
	}
	if (cgnorm <= epsx) {
	    *istop = 2;
	}
	if (cgnorm <= epsr) {
	    *istop = 1;
	}
	if (*istop != 4 && *ynorm >= *normxlim && *normxlim != 0.) {
	    *istop = 9;
	}
    }
    if (*itn < *itnmin) {
	*istop = 0;
    }
/*     ================================================================== 
*/
/*     See if it is time to print something. */
    if (*nout <= 0) {
	goto L600;
    }
    if (*n <= 40) {
	goto L400;
    }
    if (*itn <= 10) {
	goto L400;
    }
    if (*itn >= *itnlim - 10) {
	goto L400;
    }
    if (*itn % 10 == 0) {
	goto L400;
    }
    if (cgnorm <= epsx * (float)10.) {
	goto L400;
    }
    if (cgnorm <= epsr * (float)10.) {
	goto L400;
    }
    if (*acond >= (float).01 / eps) {
	goto L400;
    }
    if (*istop != 0) {
	goto L400;
    }
    goto L600;
/*     Print a line for this iteration. */
L400:
    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. */
L600:
    if (*istop != 0) {
	goto L800;
    }
    s = 1. / beta;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	v[i] = s * y[i];
/* L620: */
    }
    aprod_(n, &v[1], &y[1], a, &vwsqrt[1], &work[1], orthlist);

	
    d__1 = -(*shift);
    daxpy_(n, &d__1, &v[1], &c__1, &y[1], &c__1);
    d__1 = -beta / oldb;
    daxpy_(n, &d__1, &r1[1], &c__1, &y[1], &c__1);
    alfa = ddot_(n, &v[1], &c__1, &y[1], &c__1);
/* Computing 2nd power */
    d__1 = alfa;
/* Computing 2nd power */
    d__2 = beta;
    tnorm = tnorm + d__1 * d__1 + d__2 * d__2 * 2.;
    d__1 = -alfa / beta;
    daxpy_(n, &d__1, &r2[1], &c__1, &y[1], &c__1);
    dcopy_(n, &r2[1], &c__1, &r1[1], &c__1);
    dcopy_(n, &y[1], &c__1, &r2[1], &c__1);
    if (*precon) {
	msolve_(n, &r2[1], &y[1], a, &vwsqrt[1], &work[1]);
    }
    oldb = beta;
    beta = ddot_(n, &r2[1], &c__1, &y[1], &c__1);
    if (beta < 0.) {
	*istop = 6;
	goto L800;
    }
    beta = sqrt(beta);
/*     Compute the next plane rotation for  Q. */
/* Computing 2nd power */
    d__1 = gbar;
/* Computing 2nd power */
    d__2 = oldb;
    gamma = sqrt(d__1 * d__1 + d__2 * d__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;
    i__1 = *n;
    for (i = 1; i <= i__1; ++i) {
	x[i] = w[i] * s + v[i] * t + x[i];
	w[i] = w[i] * sn - v[i] * cs;
/* L700: */
    }
/*     Accumulate the step along the direction  b, */
/*     and go round again. */
    bstep = snprod * cs * z + bstep;
    snprod *= sn;
    gmax = max(gmax,gamma);
    gmin = min(gmin,gamma);
/* Computing 2nd power */
    d__1 = z;
    ynorm2 = d__1 * d__1 + ynorm2;
    rhs1 = rhs2 - delta * z;
    rhs2 = -epsln * z;
    ++(*itn);
    goto L300;
/*     ------------------------------------------------------------------ 
*/
/*     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. */
L800:
    if (cgnorm <= lqnorm) {
	zbar = rhs1 / diag;
	bstep = snprod * zbar + bstep;
/* Computing 2nd power */
	d__1 = zbar;
	*ynorm = sqrt(ynorm2 + d__1 * d__1);
	*rnorm = cgnorm;
	daxpy_(n, &zbar, &w[1], &c__1, &x[1], &c__1);
    } else {
	*rnorm = lqnorm;
    }
    if (*goodb) {
/*        Add the step along  b. */
	bstep /= beta1;
	dcopy_(n, &b[1], &c__1, &y[1], &c__1);
	if (*precon) {
	    msolve_(n, &b[1], &y[1], a, &vwsqrt[1], &work[1]);
	}
	daxpy_(n, &bstep, &y[1], &c__1, &x[1], &c__1);
    }
/*     ================================================================== 
*/
/*     Display final status. */
/*     ================================================================== 
*/
L900:
/*     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 0;
/*     ------------------------------------------------------------------ 
*/
/* 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 */
} /* symmlq_ */