snlaso.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 1,841 行 · 第 1/5 页

/* THIS SECTION CALLS SNPPLA (THE POST PROCESSOR). */

    if (*nperm == 0) {
        ind[0] = nop;
        return;
    }
    i1 = *n * *nblock;
    i2 = i1 + *nperm * *nperm;
    i3 = i2 + *nperm * *nperm;
    i4 = i3 + max(*n * *nblock, 2 * *nperm * *nperm);
    i5 = i4 + *n * *nblock;
    i6 = i5 + (*nperm << 1) + 4;
    snppla_(op, iovect, n, nperm, &nop, nmval, val, nmvec, vec, nblock,
            &work[i1], &work[i2], &work[i3], &work[i4], &work[i5], &work[i6],
            &delta, &small, &raritz, &eps);

    ind[0] = nop;
    return;
} /* snlaso_ */


/* ------------------------------------------------------------------ */

/* Subroutine */ void snwla_(op, iovect, n, nband, nval, nfig, nperm, val,
        nmvec, vec, nblock, maxop, maxj, nop, p1, p0, res, tau, otau, t, alp,
        bet, s, p2, bound, atemp, vtemp, d, ind, small, raritz, delta, eps, ierr)
/* Subroutine */ void (*op) (const integer*,const integer*,const real*,real*);
/* Subroutine */ void (*iovect) (const integer*,const integer*,real*,const integer*,const integer*);
const integer *n, *nband, *nval, *nfig;
integer *nperm;
real *val;
const integer *nmvec;
real *vec;
const integer *nblock, *maxop, *maxj;
integer *nop;
real *p1, *p0, *res, *tau, *otau, *t, *alp, *bet, *s, *p2, *bound, *atemp, *vtemp, *d;
integer *ind;
logical *small, *raritz;
real *delta, *eps;
integer *ierr;
{
    /* System generated locals */
    integer i__1;
    real r__1;

    /* Local variables */
    static real tola, temp, tolg, tmin, tmax, tarr;
    static logical test;
    static real zero=0.f, utol;
    static integer i, j, k, l, m;
    static integer ngood, nleft;
    static real anorm;
    static integer mtemp;
    static integer i1;
    static real pnorm, epsrt, rnorm;
    static integer ng;
    static real betmin, alpmin, betmax, alpmax;
    static integer ntheta;
    static logical enough;
    static integer number, nstart;

/* SNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE              */
/* ORTHOGONALIZATION.                                                 */
/*                                                                    */
/* NBAND  NBLOCK + 1  THE BAND WIDTH OF T.                            */
/*                                                                    */
/* NVAL   THE NUMBER OF DESIRED EIGENVALUES.                          */
/*                                                                    */
/* NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS        */
/*         INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE     */
/*         ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL  */
/*         TO THE CURRENT KRYLOV SUBSPACE.                            */
/*                                                                    */
/* NOP   THE NUMBER OF CALLS TO OP.                                   */
/*                                                                    */
/* P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS.            */
/*                                                                    */
/* RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS.   */
/*                                                                    */
/* TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION.     */
/*                                                                    */
/* T     THE BAND MATRIX.                                             */
/*                                                                    */
/* ALP   THE CURRENT DIAGONAL BLOCK.                                  */
/*                                                                    */
/* BET   THE CURRENT OFF DIAGONAL BLOCK.                              */
/*                                                                    */
/* BOUND, ATEMP, VTEMP, D                                             */
/*       TEMPORARY STORAGE USED BY THE BAND EIGENVALUE SOLVER SLAEIG. */
/*                                                                    */
/* S     EIGENVECTORS OF T.                                           */
/*                                                                    */
/* SMALL   .TRUE.  IF THE SMALL EIGENVALUES ARE DESIRED.              */
/*                                                                    */
/* RARITZ  RETURNED AS  .TRUE.  IF A FINAL RAYLEIGH-RITZ PROCEDURE    */
/*         IS TO BE DONE.                                             */
/*                                                                    */
/* DELTA   RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE         */
/*         OF THE MATRIX.  USED IN ESTIMATING THE ACCURACY OF THE     */
/*         COMPUTED EIGENVALUES.                                      */
/*                                                                    */
/* INTERNAL VARIABLES USED                                            */
/*                                                                    */
/* J      THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT  */
/*        KRYLOV SUBSPACE.                                            */
/*                                                                    */
/* NGOOD  THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS               */
/*        LIE IN THE CURRENT KRYLOV SUBSPACE).                        */
/*                                                                    */
/* NLEFT  THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED,         */
/*        I.E.  NLEFT = NVAL - NPERM.                                 */
/*                                                                    */
/* NUMBER = NPERM + NGOOD.                                            */
/*                                                                    */
/* ANORM  AN ESTIMATE OF THE NORM OF THE MATRIX.                      */
/*                                                                    */
/* EPS    THE RELATIVE MACHINE PRECISION.                             */
/*                                                                    */
/* UTOL   THE USER TOLERANCE.                                         */
/*                                                                    */
/* TARR   AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN   */
/*        CALLS TO SLAEIG                                             */
/*                                                                    */
/* ZERO   AN ARRAY OF LENGTH ONE CONTAINING ZERO, USED TO INSURE TYPE */
/*        CONSISTENCY IN CALLS TO SCOPY                               */

    rnorm = 0.f;
    if (*nperm != 0) {
        rnorm = max(-val[0],val[*nperm-1]);
    }
    pnorm = rnorm;
    *delta = 1e31f;
    epsrt = sqrtf(*eps);
    nleft = *nval - *nperm;
    *nop = 0;
    number = *nperm;
    *raritz = FALSE_;
    utol = max((*n) * *eps, pow_ri(&c__10, nfig));
    j = *maxj;

/* ------------------------------------------------------------------ */

/* ANY ITERATION OF THE ALGORITHM BEGINS HERE. */

L30:
    for (i = 0; i < *nblock; ++i) {
        temp = snrm2_(n, &p1[i * *n], &c__1);
        if (temp == 0.f) {
            slaran_(n, &p1[i * *n]);
        }
    }
    for (i = 0; i < *nperm; ++i) {
        tau[i] = 1.f;
        otau[i] = 0.f;
    }
    i__1 = *n * *nblock;
    scopy_(&i__1, &zero, &c__0, p0, &c__1);
    i__1 = *nblock * *nblock;
    scopy_(&i__1, &zero, &c__0, bet, &c__1);
    i__1 = j * *nband;
    scopy_(&i__1, &zero, &c__0, t, &c__1);
    mtemp = *nval + 1;
    for (i = 0; i < mtemp; ++i) {
        scopy_(&j, &zero, &c__0, &s[i * *maxj], &c__1);
    }
    ngood = 0;
    tmin = 1e30f;
    tmax = -1e30f;
    test = TRUE_;
    enough = FALSE_;
    betmax = 0.f;
    j = 0;

/* ------------------------------------------------------------------ */

/* THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. */

L80:
    j += *nblock;

/* THIS IS THE SELECTIVE ORTHOGONALIZATION. */

    for (i = 0; i < number; ++i) {
        if (tau[i] < epsrt) {
            continue;
        }
        test = TRUE_;
        tau[i] = 0.f;
        if (otau[i] != 0.f) {
            otau[i] = 1.f;
        }
        for (k = 0; k < *nblock; ++k) {
            temp = -sdot_(n, &vec[i * *nmvec], &c__1, &p1[k * *n], &c__1);
            saxpy_(n, &temp, &vec[i * *nmvec], &c__1, &p1[k * *n], &c__1);

/* THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A */
/* NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS */
/* TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. */

            if (abs(temp * bet[k + k * *nblock]) > (*n) * epsrt * anorm && i >= *nperm) {
                goto L380;
            }
        }
    }

/* IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. */

    if (test)
        sortqr_(n, n, nblock, p1, alp);
    if (test && j != *nblock)
    for (i = 0; i < *nblock; ++i) {
        if (alp[i + i * *nblock] > 0.f) {
            continue;
        }
        m = j - (*nblock << 1) + i;
        l = *nblock;
        for (k = i; k < *nblock; ++k, --l, ++m) {
            bet[i + k * *nblock] = -bet[i + k * *nblock];
            t[l + m * *nband] = -t[l + m * *nband];
        }
    }
    test = FALSE_;

/* THIS IS THE LANCZOS STEP. */

    (*op)(n, nblock, p1, p2);
    ++(*nop);
    (*iovect)(n, nblock, p1, &j, &c__0);

/* THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) */

    for (i = 0; i < *nblock; ++i) {
        for (k = i; k < *nblock; ++k) {
            r__1 = -bet[i + k * *nblock];
            saxpy_(n, &r__1, &p0[k * *n], &c__1, &p2[i * *n], &c__1);
        }
    }

/* THIS COMPUTES ALP AND P2=P2-P1*ALP. */

    for (i = 0; i < *nblock; ++i) {
        for (k = 0; k <= i; ++k) {
            i1 = i - k;
            alp[i1 + k * *nblock] = sdot_(n, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
            r__1 = -alp[i1 + k * *nblock];
            saxpy_(n, &r__1, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
            if (k != i) {
                r__1 = -alp[i1 + k * *nblock];
                saxpy_(n, &r__1, &p1[k * *n], &c__1, &p2[i * *n], &c__1);
            }
        }
    }

/*  REORTHOGONALIZATION OF THE SECOND BLOCK */

    if (j == *nblock)
    for (i = 0; i < *nblock; ++i) {
        for (k = 0; k <= i; ++k) {
            temp = -sdot_(n, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
            saxpy_(n, &temp, &p1[i * *n], &c__1, &p2[k * *n], &c__1);
            if (k != i) {
                saxpy_(n, &temp, &p1[k * *n], &c__1, &p2[i * *n], &c__1);
            }
            i1 = i - k;
            alp[i1 + k * *nblock] += temp;
        }
    }

/* THIS ORTHONORMALIZES THE NEXT BLOCK */

    sortqr_(n, n, nblock, p2, bet);

/* THIS STORES ALP AND BET IN T. */

    for (i = 0; i < *nblock; ++i) {
        m = j - *nblock + i;
        for (k = i; k < *nblock; ++k) {
            l = k - i;
            t[l + m * *nband] = alp[l + i * *nblock];
        }
        for (k = 0; k <= i; ++k) {
            l = *nblock - i + k;
            t[l + m * *nband] = bet[k + i * *nblock];
        }
    }

/* THIS NEGATES T IF SMALL IS FALSE. */

    if (! *small)
    for (i = j - *nblock; i < j; ++i) {
        for (k = 0; k <= l; ++k) { /* FIXME *** This must be an error! (already in the fortran code) -- l is undefined *** */
            t[k + i * *nband] = -t[k + i * *nband];
        }
    }

/* THIS SHIFTS THE LANCZOS VECTORS */

    i__1 = *nblock * *n;
    scopy_(&i__1, p1, &c__1, p0, &c__1);
    scopy_(&i__1, p2, &c__1, p1, &c__1);
    i__1 = j - *nblock + 1;
    slager_(&j, nband, &i__1, t, &tmin, &tmax);
    anorm = max(max(rnorm,tmax),-tmin);

/* THIS COMPUTES THE EXTREME EIGENVALUES OF ALP. */

    if (number != 0) {
        scopy_(nblock, &zero, &c__0, p2, &c__1);
        slaeig_(nblock, nblock, &c__1, &c__1, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &tmin, &tmax);
        alpmin = tarr;
        scopy_(nblock, &zero, &c__0, p2, &c__1);
        slaeig_(nblock, nblock, nblock, nblock, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &tmin, &tmax);
        alpmax = tarr;
    }

/* THIS COMPUTES ALP = BET(TRANSPOSE)*BET. */

    for (i = 0; i < *nblock; ++i) {
        for (k = 0; k <= i; ++k) {
            l = i - k;
            i__1 = *nblock - i;
            alp[l + k * *nblock] = sdot_(&i__1, &bet[i + i * *nblock], nblock, &bet[k + i * *nblock], nblock);
        }
    }
    if (number == 0) {
        goto L330;
    }

/* THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET. */

    scopy_(nblock, &zero, &c__0, p2, &c__1);
    r__1 = anorm * anorm;
    slaeig_(nblock, nblock, &c__1, &c__1, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &c__00, &r__1);
    betmin = sqrtf(tarr);

/* THIS UPDATES TAU AND OTAU. */

    for (i = 0; i < number; ++i) {
        temp = (tau[i] * max(alpmax-val[i],val[i]-alpmin) + otau[i] * betmax + *eps * anorm) / betmin;
        if (i < *nperm) {
            temp += res[i] / betmin;
        }
        otau[i] = tau[i];
        tau[i] = temp;
    }

/* THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET. */

L330:
    scopy_(nblock, &zero, &c__0, p2, &c__1);
    r__1 = anorm * anorm;
    slaeig_(nblock, nblock, nblock, nblock, alp, &tarr, nblock, p2, bound, atemp, d, vtemp, eps, &c__00, &r__1);
    betmax = sqrtf(tarr);
    if (j <= *nblock << 1) {
        goto L80;
    }

/* ------------------------------------------------------------------ */

/* THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND */
/* LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK */
/* IS JUSTIFIED. */

    tolg = epsrt * anorm;
    tola = utol * rnorm;
    if (*maxj - j < *nblock || ( *nop >= *maxop && nleft != 0 ) ) {
        goto L390;
    }
    else