dnlaso.c

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

    nval = *nr - *nl + 1;
    flag = FALSE_;
    for (i = 0; i < nval; ++i) {
        if (ddot_(n, &eigvec[i * *lde], &c__1, &eigvec[i * *lde], &c__1) == 0.) {
            dlaran_(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:
        dlabax_(n, nband, a, &eigvec[j * *lde], vtemp);
        vnorm = dnrm2_(n, vtemp, &c__1);
        if (vnorm != 0.) {
            d__1 = 1. / vnorm;
            dscal_(n, &d__1, vtemp, &c__1);
            dscal_(n, &d__1, &eigvec[j * *lde], &c__1);
            d__1 = -sigma;
            daxpy_(n, &d__1, &eigvec[j * *lde], &c__1, vtemp, &c__1);
        }

/*  LOOP OVER SHIFTS */

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

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

        rq = sigma + ddot_(n, &eigvec[j * *lde], &c__1, vtemp, &c__1) / vnorm / vnorm;
        d__1 = sigma - rq;
        daxpy_(n, &d__1, &eigvec[j * *lde], &c__1, vtemp, &c__1);
        resid = max(*atol,dnrm2_(n, vtemp, &c__1) / vnorm);
        d__1 = 1. / vnorm;
        dscal_(n, &d__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.) {
            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]) * .5;

/*  CHECK FOR TERMINALTION */

        if (resid > *atol * 2.) {
            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]) {
            dlaran_(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) {
                dcopy_(n, &eigvec[j * *lde], &c__1, &eigvec[i * *lde], &c__1);
                break;
            }
        }
        dlaran_(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);
        }
        dcopy_(n, &eigvec[j * *lde], &c__1, vtemp, &c__1);
        i__1 = (*nband << 1) - 1;
        dlabfc_(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) {
            d__1 = 1. / vnorm;
            dscal_(n, &d__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:
        dlaran_(n, &eigvec[j * *lde]);
        flag = TRUE_;
        goto L310;

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

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

L310:
        eigval[j] = rq;
        for (i = 0; i < j; ++i) {
            d__1 = -ddot_(n, &eigvec[i * *lde], &c__1, &eigvec[j * *lde], &c__1);
            daxpy_(n, &d__1, &eigvec[i * *lde], &c__1, &eigvec[j * *lde], &c__1);
        }
        vnorm = dnrm2_(n, &eigvec[j * *lde], &c__1);
        if (vnorm == 0.) {
            goto L305;
        }
        d__1 = 1. / vnorm;
        dscal_(n, &d__1, &eigvec[j * *lde], &c__1);

/*   ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE */

        if (flag) {
            goto L305;
        }
        for (i = j+1; i < nval; ++i) {
            d__1 = -ddot_(n, &eigvec[j * *lde], &c__1, &eigvec[i * *lde], &c__1);
            daxpy_(n, &d__1, &eigvec[j * *lde], &c__1, &eigvec[i * *lde], &c__1);
        }
    }
} /* dlabcm_ */


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

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

    /* Local variables */
    static doublereal zero=0.;
    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;
    dcopy_(&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.) {
                continue;
            }
            if (abs(d[i + i * *ldad]) >= abs(d[i + nb1 * *ldad])) {
                goto L50;
            }
            if ( (d[i + nb1 * *ldad] < 0. && d[i + i * *ldad] < 0. ) ||
                 (d[i + nb1 * *ldad] > 0. && d[i + i * *ldad] >= 0.) ) {
                lpm = -lpm;
            }
            i__1 = *ldad - i;
            dswap_(&i__1, &d[i + i * *ldad], &c__1, &d[i + nb1 * *ldad], &c__1);
            dswap_(number, &eigvec[l], lde, &eigvec[k], lde);
L50:
            i__1 = *ldad - i - 1;
            d__1 = -d[i + nb1 * *ldad] / d[i + i * *ldad];
            daxpy_(&i__1, &d__1, &d[i + 1 + i * *ldad], &c__1, &d[i + 1 + nb1 * *ldad], &c__1);
            d__1 = -d[i + nb1 * *ldad] / d[i + i * *ldad];
            daxpy_(number, &d__1, &eigvec[l], lde, &eigvec[k], lde);
        }

/*  UPDATE STURM SEQUENCE COUNT */

        if (d[nb1 + nb1 * *ldad] < 0.) {
            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;
            dcopy_(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;
            dcopy_(&i__1, &d[i + 1 + (i + 1) * *ldad], &c__1, &d[i + i * *ldad], &c__1);
            d[*ldad - 1 + i * *ldad] = 0.;
        }
    }

/*  TRANSFER D TO ATEMP */

L110:
    for (i = 0; i < *nband; ++i) {
        i__1 = *nband - i;
        l = *n - i__1;
        dcopy_(&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] = d_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];
            }
        }
    }
} /* dlabfc_ */


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