snlaso.c

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

            r__1 = -sigma;
            saxpy_(n, &r__1, &eigvec[j * *lde], &c__1, vtemp, &c__1);
        }

/*  LOOP OVER SHIFTS */

/*  COMPUTE RAYLEIGH QUOTIENT, RESIDUAL NORM, AND CURRENT TOLERANCE */

L20:
        vnorm = snrm2_(n, &eigvec[j * *lde], &c__1);
        if (vnorm == 0.f) {
            slaran_(n, &eigvec[j * *lde]);
            goto L10;
        }

        rq = sigma + sdot_(n, &eigvec[j * *lde], &c__1, vtemp, &c__1) / vnorm / vnorm;
        r__1 = sigma - rq;
        saxpy_(n, &r__1, &eigvec[j * *lde], &c__1, vtemp, &c__1);
        resid = max(*atol,snrm2_(n, vtemp, &c__1) / vnorm);
        r__1 = 1.f / vnorm;
        sscal_(n, &r__1, &eigvec[j * *lde], &c__1);

/*  ACCEPT EIGENVALUE IF THE INTERVAL IS SMALL ENOUGH */

        if (bound[(j << 1) + 3] - bound[(j << 1) + 2] < *atol * 3.f) {
            goto L300;
        }

/*  COMPUTE MINIMAL ERROR BOUND */

        errb = resid;
        gap = min(bound[(j << 1) + 4] - rq,rq - bound[(j << 1) + 1]);
        if (gap > resid) {
            errb = max(*atol,resid * resid / gap);
        }

/*  TENTATIVE NEW SHIFT */

        sigma = (bound[(j << 1) + 2] + bound[(j << 1) + 3]) * .5f;

/*  CHECK FOR TERMINALTION */

        if (resid > *atol * 2.f) {
            goto L40;
        }
        if (rq - errb > bound[(j << 1) + 1] && rq + errb < bound[(j << 1) + 4]) {
            goto L310;
        }

/*  RQ IS TO THE LEFT OF THE INTERVAL */

L40:
        if (rq >= bound[(j << 1) + 2]) {
            goto L50;
        }
        if (rq - errb > bound[(j << 1) + 1]) {
            goto L100;
        }
        if (rq + errb < bound[(j << 1) + 2]) {
            slaran_(n, &eigvec[j * *lde]);
        }
        goto L200;

/*  RQ IS TO THE RIGHT OF THE INTERVAL */

L50:
        if (rq <= bound[(j << 1) + 3]) {
            goto L100;
        }
        if (rq + errb < bound[(j << 1) + 4]) {
            goto L100;
        }

/*  SAVE THE REJECTED VECTOR IF INDICATED */

        if (rq - errb <= bound[(j << 1) + 3]) {
            goto L200;
        }
        for (i = j; i < nval; ++i) {
            if (bound[(i << 1) + 3] > rq) {
                scopy_(n, &eigvec[j * *lde], &c__1, &eigvec[i * *lde], &c__1);
                break;
            }
        }
        slaran_(n, &eigvec[j * *lde]);
        goto L200;

/*  PERTURB RQ TOWARD THE MIDDLE */

L100:
        if (sigma < rq-errb) {
            sigma = rq-errb;
        }
        if (sigma > rq+errb) {
            sigma = rq+errb;
        }

/*  FACTOR AND SOLVE */

L200:
        for (i = j; i < nval; ++i) {
            if (sigma < bound[(i << 1) + 2]) {
                break;
            }
        }
        numvec = i - j;
        numvec = min(numvec,*nband+2);
        if (resid < *artol) {
            numvec = min(1,numvec);
        }
        scopy_(n, &eigvec[j * *lde], &c__1, vtemp, &c__1);
        i__1 = (*nband << 1) - 1;
        slabfc_(n, nband, a, &sigma, &numvec, lde, &eigvec[j * *lde], &numl, &i__1, atemp, d, atol);

/*  PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW */

        for (i = j+1; i < numvec+j; ++i) {
            r__1 = 1.f / vnorm;
            sscal_(n, &r__1, &eigvec[i * *lde], &c__1);
        }

/*  UPDATE INTERVALS */

        numl -= *nl - 1;
        if (numl >= 0) {
            bound[1] = min(bound[1],sigma);
        }
        for (i = j; i < nval; ++i) {
            if (sigma < bound[(i << 1) + 2]) {
                goto L20;
            }
            if (numl <= i)
                bound[(i << 1) + 2] = sigma;
            else
                bound[(i << 1) + 3] = sigma;
        }
        if (numl < nval + 1) {
            if (sigma > bound[(nval << 1) + 2])
                bound[(nval << 1) + 2] = sigma;
        }
        goto L20;

/*  ACCEPT AN EIGENPAIR */

L300:
        slaran_(n, &eigvec[j * *lde]);
        flag = TRUE_;
        goto L310;

L305:
        flag = FALSE_;
        rq = (bound[(j << 1) + 2] + bound[(j << 1) + 3]) * .5f;
        i__1 = (*nband << 1) - 1;
        slabfc_(n, nband, a, &rq, &numvec, lde, &eigvec[j * *lde], &numl, &i__1, atemp, d, atol);
        vnorm = snrm2_(n, &eigvec[j * *lde], &c__1);
        if (vnorm != 0.f) {
            r__1 = 1.f / vnorm;
            sscal_(n, &r__1, &eigvec[j * *lde], &c__1);
        }

/*  ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES */

L310:
        eigval[j] = rq;
        for (i = 0; i < j; ++i) {
            r__1 = -sdot_(n, &eigvec[i * *lde], &c__1, &eigvec[j * *lde], &c__1);
            saxpy_(n, &r__1, &eigvec[i * *lde], &c__1, &eigvec[j * *lde], &c__1);
        }
        vnorm = snrm2_(n, &eigvec[j * *lde], &c__1);
        if (vnorm == 0.f) {
            goto L305;
        }
        r__1 = 1.f / vnorm;
        sscal_(n, &r__1, &eigvec[j * *lde], &c__1);

/*   ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE */

        if (flag) {
            goto L305;
        }
        for (i = j+1; i < nval; ++i) {
            r__1 = -sdot_(n, &eigvec[j * *lde], &c__1, &eigvec[i * *lde], &c__1);
            saxpy_(n, &r__1, &eigvec[j * *lde], &c__1, &eigvec[i * *lde], &c__1);
        }
    }
} /* slabcm_ */


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

/* Subroutine */ void slabfc_(n, nband, a, sigma, number, lde, eigvec, numl, ldad, atemp, d, atol)
const integer *n, *nband;
real *a, *sigma;
const integer *number, *lde;
real *eigvec;
integer *numl, *ldad;
real *atemp, *d, *atol;
{
    /* System generated locals */
    integer i__1;
    real r__1;

    /* Local variables */
    static real zero=0.f;
    static integer i, j, k, l, m;
    static integer la, ld, nb1, lpm;

/*  THIS SUBROUTINE FACTORS (A-SIGMA*I) WHERE A IS A GIVEN BAND */
/*  MATRIX AND SIGMA IS AN INPUT PARAMETER.  IT ALSO SOLVES ZERO */
/*  OR MORE SYSTEMS OF LINEAR EQUATIONS.  IT RETURNS THE NUMBER */
/*  OF EIGENVALUES OF A LESS THAN SIGMA BY COUNTING THE STURM */
/*  SEQUENCE DURING THE FACTORIZATION.  TO OBTAIN THE STURM */
/*  SEQUENCE COUNT WHILE ALLOWING NON-SYMMETRIC PIVOTING FOR */
/*  STABILITY, THE CODE USES A GUPTA'S MULTIPLE PIVOTING */
/*  ALGORITHM. */

/*  INITIALIZE */

    nb1 = *nband - 1;
    *numl = 0;
    i__1 = *ldad * *nband;
    scopy_(&i__1, &zero, &c__0, d, &c__1);

/*   LOOP OVER COLUMNS OF A */

    for (k = 0; k < *n; ++k) {

/*   ADD A COLUMN OF A TO D */

        d[nb1 + nb1 * *ldad] = a[k * *nband] - *sigma;
        m = min(k,nb1);
        for (i = 0; i < m; ++i) {
            la = k - i - 1;
            ld = nb1 - i - 1;
            d[ld + nb1 * *ldad] = a[i + 1 + la * *nband];
        }

        m = min(*n-k-1,nb1);
        for (i = 0; i < m; ++i) {
            ld = *nband + i;
            d[ld + nb1 * *ldad] = a[i + 1 + k * *nband];
        }

/*   TERMINATE */

        lpm = 1;
        for (i = 0; i < nb1; ++i) {
            l = k - nb1 + i;
            if (d[i + nb1 * *ldad] == 0.f) {
                continue;
            }
            if (abs(d[i + i * *ldad]) >= abs(d[i + nb1 * *ldad])) {
                goto L50;
            }
            if ( (d[i + nb1 * *ldad] < 0.f && d[i + i * *ldad] < 0.f ) ||
                 (d[i + nb1 * *ldad] > 0.f && d[i + i * *ldad] >= 0.f) ) {
                lpm = -lpm;
            }
            i__1 = *ldad - i;
            sswap_(&i__1, &d[i + i * *ldad], &c__1, &d[i + nb1 * *ldad], &c__1);
            sswap_(number, &eigvec[l], lde, &eigvec[k], lde);
L50:
            i__1 = *ldad - i - 1;
            r__1 = -d[i + nb1 * *ldad] / d[i + i * *ldad];
            saxpy_(&i__1, &r__1, &d[i + 1 + i * *ldad], &c__1, &d[i + 1 + nb1 * *ldad], &c__1);
            r__1 = -d[i + nb1 * *ldad] / d[i + i * *ldad];
            saxpy_(number, &r__1, &eigvec[l], lde, &eigvec[k], lde);
        }

/*  UPDATE STURM SEQUENCE COUNT */

        if (d[nb1 + nb1 * *ldad] < 0.f) {
            lpm = -lpm;
        }
        if (lpm < 0) {
            ++(*numl);
        }
        if (k == *n-1) {
            goto L110;
        }

/*   COPY FIRST COLUMN OF D INTO ATEMP */
        if (k >= nb1) {
            l = k - nb1;
            scopy_(ldad, d, &c__1, &atemp[l * *ldad], &c__1);
        }

/*   SHIFT THE COLUMNS OF D OVER AND UP */

        for (i = 0; i < nb1; ++i) {
            i__1 = *ldad - i - 1;
            scopy_(&i__1, &d[i + 1 + (i + 1) * *ldad], &c__1, &d[i + i * *ldad], &c__1);
            d[*ldad - 1 + i * *ldad] = 0.f;
        }
    }

/*  TRANSFER D TO ATEMP */

L110:
    for (i = 0; i < *nband; ++i) {
        i__1 = *nband - i;
        l = *n - i__1;
        scopy_(&i__1, &d[i + i * *ldad], &c__1, &atemp[l * *ldad], &c__1);
    }

/*   BACK SUBSTITUTION */

    if (*number == 0) {
        return;
    }
    for (k = *n-1; k >= 0; --k) {
        if (abs(atemp[k * *ldad]) <= *atol) {
            atemp[k * *ldad] = r_sign(atol, &atemp[k * *ldad]);
        }

        for (i = 0; i < *number; ++i) {
            eigvec[k + i * *lde] /= atemp[k * *ldad];
            m = min(*ldad-1,k);
            for (j = 0; j < m; ++j) {
                l = k - j - 1;
                eigvec[l + i * *lde] -= atemp[j + 1 + l * *ldad] * eigvec[k + i * *lde];
            }
        }
    }
} /* slabfc_ */


/* Subroutine */ void slaeig_(n, nband, nl, nr, a, eigval, lde, eigvec, bound, atemp, d, vtemp, eps, tmin, tmax)
const integer *n, *nband, *nl, *nr;
real *a, *eigval;
const integer *lde;
real *eigvec, *bound, *atemp, *d, *vtemp, *eps, *tmin, *tmax;
{
    /* Local variables */
    static real atol;
    static integer nval, i;
    static real artol;

/*  THIS IS A SPECIALIZED VERSION OF THE SUBROUTINE BNDEIG TAILORED */
/*  SPECIFICALLY FOR USE BY THE LASO PACKAGE. */

/*  SET PARAMETERS */

    atol = *n * *eps * max(*tmax,-(*tmin));
    artol = atol / sqrtf(*eps);
    nval = *nr - *nl + 1;

/*   CHECK FOR SPECIAL CASE OF N = 1 */

    if (*n == 1) {
        eigval[0] = a[0];
        eigvec[0] = 1.f;
        return;
    }

/*   SET UP INITIAL EIGENVALUE BOUNDS */

    for (i = 1; i <= nval; ++i) {
        bound[(i << 1)] = *tmin;
        bound[(i << 1) + 1] = *tmax;
    }
    bound[1] = *tmax;
    bound[(nval << 1) + 2] = *tmin;
    if (*nl == 1) {
        bound[1] = *tmin;
    }
    if (*nr == *n) {
        bound[(nval << 1) + 2] = *tmax;
    }