snlaso.c

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

        goto L400;

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

/* THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO */
/* SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN */
/* ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY */
/* IS AT 400. */

L380:
    j -= *nblock;
    *ierr = -8;
L390:
    if (nleft == 0) {
        return;
    }
    test = TRUE_;
L400:
    ntheta = min(j/2, nleft+1);
    slaeig_(&j, nband, &c__1, &ntheta, t, &val[number], maxj, s, bound, atemp, d, vtemp, eps, &tmin, &tmax);
    smvpc_(nblock, bet, maxj, &j, s, &ntheta, atemp, vtemp, d);

/* THIS CHECKS FOR TERMINATION OF A CHECK RUN */

    if (nleft == 0 && j >= *nblock * 6) {
        if (val[number] - atemp[0] > val[*nperm-1] - tola) {
            goto L790;
        }
    }

/* THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T */
/* TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED. */

    if (ntheta <= nleft) {
        goto L470;
    }
    if (*nperm != 0 && val[number+nleft] < val[*nperm-1]) {
        --(*nperm);
        ngood = 0;
        number = *nperm;
        ++nleft;
        goto L400;
    }

/* THIS UPDATES DELTA. */

    *delta = min(*delta,val[number+nleft]);
    enough = TRUE_;
    if (nleft == 0) {
        goto L80;
    }
    ntheta = nleft;
    vtemp[ntheta] = 1.f;

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

/* THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL. */

/* THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES. */

    if (! (test || enough)) {
        goto L470;
    }
    *delta = min(*delta,anorm);
    pnorm = max(rnorm,max(-val[number],*delta));
    tola = utol * pnorm;
    nstart = 0;
    for (i = 0; i < ntheta; ++i) {
        if (min(atemp[i]*atemp[i]/(*delta-val[number+i]), atemp[i]) <= tola) {
            ind[i] = -1;
            continue;
        }
        enough = FALSE_;
        if (! test) {
            goto L470;
        }
        ind[i] = 1;
        ++nstart;
    }

/*  COPY VALUES OF IND INTO VTEMP */

    for (i = 0; i < ntheta; ++i) {
        vtemp[i] = (real) ind[i];
    }
    goto L500;

/* THIS CHECKS FOR NEW GOOD VECTORS. */

L470:
    ng = 0;
    for (i = 0; i < ntheta; ++i) {
        if (vtemp[i] > tolg) {
            vtemp[i] = 1.f;
        }
        else {
            ++ng;
            vtemp[i] = -1.f;
        }
    }

    if (ng <= ngood) {
        goto L80;
    }
    nstart = ntheta - ng;

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

/* THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS. */
/* IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED. */

L500:
    test = test && ! enough;
    ngood = ntheta - nstart;
    ++nstart;
    ++ntheta;

/* THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND */
/* EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR */
/* FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE */
/* EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS */
/* PAUSE. */

    scopy_(&ntheta, &val[number], &c__1, &val[*nperm], &c__1);
    if (nstart == 0) {
        goto L580;
    }
    if (nstart != ntheta) {
        svsort_(&ntheta, vtemp, atemp, &c__1, &val[*nperm], maxj, &j, s);
    }

/* THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING */
/* VECTORS. */

    if (! test) {
        nstart = 0;
    }
    if (! test) {
        goto L580;
    }

/*  FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW */

    temp = atemp[0];
    for (i = 0; i < nstart; ++i) {
        temp = min(temp,atemp[i]);
    }
    l = ngood + min(nstart,*nblock);
    for (i = ngood; i < l; ++i) {
        r__1 = temp / atemp[i];
        sscal_(&j, &r__1, &s[i * *maxj], &c__1);
    }
    m = (nstart - 1) / *nblock;
    l = ngood + *nblock;
    for (i = 0; i < m; ++i) {
        for (k = 0; k < *nblock; ++k, ++l) {
            if (l >= ntheta) {
                goto L570;
            }
            i1 = ngood + k;
            r__1 = temp / atemp[l];
            saxpy_(&j, &r__1, &s[l * *maxj], &c__1, &s[i1 * *maxj], &c__1);
        }
    }
L570:
    nstart = min(nstart,*nblock);

/* THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS. */

L580:
    if (test || enough)
    for (i = 0; i < ngood; ++i) {
        res[*nperm+i] = atemp[i];
    }

/* THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE */
/* LANCZOS VECTORS. */

    number = *nperm + ngood;
    if (test || enough) {
        i__1 = *n * *nblock;
        scopy_(&i__1, &zero, &c__0, p1, &c__1);
    }
    if (ngood != 0)
    for (i = *nperm; i < number; ++i) {
        scopy_(n, &zero, &c__0, &vec[i * *nmvec], &c__1);
    }
    for (i = *nblock; *nblock < 0 ? i >= j : i <= j; i += *nblock) {
        (*iovect)(n, nblock, p2, &i, &c__1);
        for (k = 0; k < *nblock; ++k) {
            m = i - *nblock + k;
            for (l = 0; l < nstart; ++l) {
                i1 = ngood + l;
                saxpy_(n, &s[m + i1 * *maxj], &p2[k * *n], &c__1, &p1[l * *n], &c__1);
            }
            for (l = 0; l < ngood; ++l) {
                i1 = l + *nperm;
                saxpy_(n, &s[m + l * *maxj], &p2[k * *n], &c__1, &vec[i1 * *nmvec], &c__1);
            }
        }
    }
    if (test || enough) {
        goto L690;
    }

/* THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE */
/* TAU RECURRENCE. */

    for (i = *nperm; i < number; ++i) {
        temp = 1.f / snrm2_(n, &vec[i * *nmvec], &c__1);
        sscal_(n, &temp, &vec[i * *nmvec], &c__1);
        tau[i] = 1.f;
        otau[i] = 1.f;
    }

/*  SHIFT S VECTORS TO ALIGN FOR LATER CALL TO SLAEIG */

    scopy_(&ntheta, &val[*nperm], &c__1, vtemp, &c__1);
    svsort_(&ntheta, vtemp, atemp, &c__0, &tarr, maxj, &j, s);
    goto L80;

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

/* THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE */
/* PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING */
/* THE PERMANENT VECTORS. */

L690:
    if (ngood == 0 && *nop >= *maxop) {
        *ierr = -2; /* THIS REPORTS THAT MAXOP WAS EXCEEDED. */
        goto L790;
    }
    if (ngood == 0) {
        goto L30;
    }

/* THIS ORTHONORMALIZES THE VECTORS */

    i__1 = *nperm + ngood;
    sortqr_(nmvec, n, &i__1, vec, s);

/* THIS SORTS THE VALUES AND VECTORS. */

    if (*nperm != 0) {
        i__1 = *nperm + ngood;
        svsort_(&i__1, val, res, &c__0, &temp, nmvec, n, vec);
    }
    *nperm += ngood;
    nleft -= ngood;
    rnorm = max(-val[0],val[*nperm-1]);

/* THIS DECIDES WHERE TO GO NEXT. */

    if (*nop >= *maxop && nleft != 0) {
        *ierr = -2; /* THIS REPORTS THAT MAXOP WAS EXCEEDED. */
        goto L790;
    }
    if (nleft != 0) {
        goto L30;
    }
    if (val[*nval-1] - val[0] < tola) {
        goto L790;
    }

/* THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED */
/* TO LOOK FOR UNDISCLOSED MULTIPLICITIES. */

    m = *nperm - *nblock;
    for (i = 0; i <= m; ++i) {
        if (val[i + *nblock - 1] - val[i] < tola) {
            goto L30;
        }
    }

/* THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ */
/* PROCEDURE IS NEEDED. */

L790:
    m = *nperm - *nblock;
    for (i = 0; i < m; ++i) {
        if (val[i + *nblock] - val[i] < tola) {
            *raritz = TRUE_;
            break;
        }
    }
} /* snwla_ */


/* *********************************************************************** */

/* Subroutine */ void slabax_(n, nband, a, x, y)
const integer *n, *nband;
real *a, *x, *y;
{
    /* Local variables */
    static real zero = 0.f;
    static integer i, k, m;

/*  THIS SUBROUTINE SETS Y = A*X */
/*  WHERE X AND Y ARE VECTORS OF LENGTH N */
/*  AND A IS AN  N X NBAND  SYMMETRIC BAND MATRIX */

    scopy_(n, &zero, &c__0, y, &c__1);
    for (k = 0; k < *n; ++k) {
        y[k] += a[k * *nband] * x[k];
        m = min(*n-k,*nband);
        for (i = 1; i < m; ++i) {
            y[k+i] += a[i + k * *nband] * x[k];
            y[k]   += a[i + k * *nband] * x[k+i];
        }
    }
} /* slabax_ */


/* *********************************************************************** */

/* Subroutine */ void slabcm_(n, nband, nl, nr, a, eigval, lde, eigvec, atol, artol, bound, atemp, d, vtemp)
const integer *n, *nband, *nl, *nr;
real *a, *eigval;
const integer *lde;
real *eigvec, *atol, *artol, *bound, *atemp, *d, *vtemp;
{
    /* System generated locals */
    integer i__1;
    real r__1;

    /* Local variables */
    static logical flag;
    static real errb;
    static integer nval, numl;
    static integer i, j;
    static real sigma, resid;
    static real vnorm;
    static real rq;
    static integer numvec;
    static real gap;

/*  THIS SUBROUTINE ORGANIZES THE CALCULATION OF THE EIGENVALUES */
/*  FOR THE BNDEIG PACKAGE.  EIGENVALUES ARE COMPUTED BY */
/*  A MODIFIED RAYLEIGH QUOTIENT ITERATION.  THE EIGENVALUE COUNT */
/*  OBTAINED BY EACH FACTORIZATION IS USED TO OCCASIONALLY OVERRIDE */
/*  THE COMPUTED RAYLEIGH QUOTIENT WITH A DIFFERENT SHIFT TO */
/*  INSURE CONVERGENCE TO THE DESIRED EIGENVALUES. */

/*  REPLACE ZERO VECTORS BY RANDOM */

    nval = *nr - *nl + 1;
    flag = FALSE_;
    for (i = 0; i < nval; ++i) {
        if (sdot_(n, &eigvec[i * *lde], &c__1, &eigvec[i * *lde], &c__1) == 0.f) {
            slaran_(n, &eigvec[i * *lde]);
        }
    }

/*  LOOP OVER EIGENVALUES */

    sigma = bound[(nval << 1) + 1];
    for (j = 0; j < nval; ++j) {
        numl = j+1;

/*  PREPARE TO COMPUTE FIRST RAYLEIGH QUOTIENT */

L10:
        slabax_(n, nband, a, &eigvec[j * *lde], vtemp);
        vnorm = snrm2_(n, vtemp, &c__1);
        if (vnorm != 0.f) {
            r__1 = 1.f / vnorm;
            sscal_(n, &r__1, vtemp, &c__1);
            sscal_(n, &r__1, &eigvec[j * *lde], &c__1);