dnlaso.c

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

/*  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 / sqrt(*eps);
    nval = *nr - *nl + 1;

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

    if (*n == 1) {
        eigval[0] = a[0];
        eigvec[0] = 1.;
        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;
    }

    dlabcm_(n, nband, nl, nr, a, eigval, lde, eigvec, &atol, &artol, bound, atemp, d, vtemp);
} /* dlaeig_ */


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

/* Subroutine */ void dlager_(n, nband, nstart, a, tmin, tmax)
const integer *n, *nband, *nstart;
doublereal *a, *tmin, *tmax;
{
    /* Local variables */
    static doublereal temp;
    static integer i, k, l;

/*  THIS SUBROUTINE COMPUTES BOUNDS ON THE SPECTRUM OF A BY */
/*  EXAMINING THE GERSCHGORIN CIRCLES. ONLY THE NEWLY CREATED */
/*  CIRCLES ARE EXAMINED */

    for (k = *nstart - 1; k < *n; ++k) {
        temp = 0.;
        for (i = 1; i < *nband; ++i) {
            temp += abs(a[i + k * *nband]);
        }
        l = min(k,*nband-1);
        for (i = 1; i <= l; ++i) {
            temp += abs(a[i + (k-i) * *nband]);
        }
        *tmin = min(*tmin,a[k * *nband] - temp);
        *tmax = max(*tmax,a[k * *nband] + temp);
    }
} /* dlager_ */


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

/* Subroutine */ void dlaran_(n, x)
const integer *n;
doublereal *x;
{
    /* Initialized data */
    static integer iurand = 0;

    /* Local variables */
    static integer i;

/*  THIS SUBROUTINE SETS THE VECTOR X TO RANDOM NUMBERS */

/*  INITIALIZE SEED */

    for (i = 0; i < *n; ++i) {
        x[i] = urand_(&iurand) - .5;
    }
} /* dlaran_ */


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

/* Subroutine */ void dmvpc_(nblock, bet, maxj, j, s, number, resnrm, orthcf, rv)
const integer *nblock;
const doublereal *bet;
const integer *maxj, *j;
const doublereal *s;
const integer *number;
doublereal *resnrm, *orthcf, *rv;
{
    /* Local variables */
    static integer i, k, m;

/* THIS SUBROUTINE COMPUTES THE NORM AND THE SMALLEST ELEMENT */
/* (IN ABSOLUTE VALUE) OF THE VECTOR BET*SJI, WHERE SJI */
/* IS AN NBLOCK VECTOR OF THE LAST NBLOCK ELEMENTS OF THE ITH */
/* EIGENVECTOR OF T.  THESE QUANTITIES ARE THE RESIDUAL NORM */
/* AND THE ORTHOGONALITY COEFFICIENT RESPECTIVELY FOR THE */
/* CORRESPONDING RITZ PAIR.  THE ORTHOGONALITY COEFFICIENT IS */
/* NORMALIZED TO ACCOUNT FOR THE LOCAL REORTHOGONALIZATION. */

    m = *j - *nblock;
    for (i = 0; i < *number; ++i) {
        rv[0] = ddot_(nblock, &s[m + i * *maxj], &c__1, &bet[0], nblock);
        orthcf[i] = abs(rv[0]);
        for (k = 1; k < *nblock; ++k) {
            rv[k] = ddot_(nblock, &s[m + i * *maxj], &c__1, &bet[k], nblock);
            orthcf[i] = min(orthcf[i], abs(rv[k]));
        }
        resnrm[i] = dnrm2_(nblock, rv, &c__1);
    }
} /* dmvpc_ */


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

/* Subroutine */ void dnppla_(op, iovect, n, nperm, nop, nmval, val, nmvec,
        vec, nblock, h, hv, p, q, bound, d, delta, small, raritz, eps)
/* Subroutine */ void (*op) (const integer*,const integer*,const doublereal*,doublereal*);
/* Subroutine */ void (*iovect) (const integer*,const integer*,doublereal*,const integer*,const integer*);
const integer *n, *nperm, *nmval;
integer *nop;
doublereal *val;
const integer *nmvec;
doublereal *vec;
const integer *nblock;
doublereal *h, *hv, *p, *q, *bound, *d, *delta;
logical *small, *raritz;
doublereal *eps;
{
    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Local variables */
    static doublereal hmin, hmax, temp;
    static doublereal zero=0.;
    static integer i, j, k, l, m;
    static integer jj, kk;

/* THIS SUBROUTINE POST PROCESSES THE EIGENVECTORS.  BLOCK MATRIX */
/* VECTOR PRODUCTS ARE USED TO MINIMIZED THE NUMBER OF CALLS TO OP. */

/* IF RARITZ IS .TRUE.  A FINAL RAYLEIGH-RITZ PROCEDURE IS APPLIED */
/* TO THE EIGENVECTORS. */

    if (! (*raritz)) {
        goto L190;
    }

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

/* THIS CONSTRUCTS H=Q*AQ, WHERE THE COLUMNS OF Q ARE THE */
/* APPROXIMATE EIGENVECTORS.  TEMP = -1 IS USED WHEN SMALL IS */
/* FALSE TO AVOID HAVING TO RESORT THE EIGENVALUES AND EIGENVECTORS */
/* COMPUTED BY DLAEIG. */

    i__1 = *nperm * *nperm;
    dcopy_(&i__1, &zero, &c__0, h, &c__1);
    temp = -1.;
    if (*small) {
        temp = 1.;
    }
    m = *nperm % *nblock;
    if (m == 0) {
        goto L40;
    }
    for (i = 0; i < m; ++i) {
        dcopy_(n, &vec[i * *nmvec], &c__1, &p[i * *n], &c__1);
    }
    (*iovect)(n, &m, p, &m, &c__0);
    (*op)(n, &m, p, q);
    ++(*nop);
    for (i = 0; i < m; ++i) {
        for (j = i; j < *nperm; ++j) {
            jj = j - i;
            h[jj + i * *nperm] = temp * ddot_(n, &vec[j * *nmvec], &c__1, &q[i * *n], &c__1);
        }
    }
    if (*nperm < *nblock) {
        goto L90;
    }
L40:
    m += *nblock;
    for (i = m; *nblock < 0 ? i >= *nperm : i <= *nperm; i += *nblock) {
        for (j = 0; j < *nblock; ++j) {
            l = i - *nblock + j;
            dcopy_(n, &vec[l * *nmvec], &c__1, &p[j * *n], &c__1);
        }
        (*iovect)(n, nblock, p, &i, &c__0);
        (*op)(n, nblock, p, q);
        ++(*nop);
        for (j = 0; j < *nblock; ++j) {
            l = i - *nblock + j;
            for (k = l; k < *nperm; ++k) {
                kk = k - l;
                h[kk + l * *nperm] = temp * ddot_(n, &vec[k * *nmvec], &c__1, &q[j * *n], &c__1);
            }
        }
    }

/* THIS COMPUTES THE SPECTRAL DECOMPOSITION OF H. */

L90:
    hmin = h[0];
    hmax = h[0];
    dlager_(nperm, nperm, &c__1, h, &hmin, &hmax);
    dlaeig_(nperm, nperm, &c__1, nperm, h, val, nperm, hv, bound, p, d, q, eps, &hmin, &hmax);

/* THIS COMPUTES THE RITZ VECTORS--THE COLUMNS OF */
/* Y = QS WHERE S IS THE MATRIX OF EIGENVECTORS OF H. */

    for (i = 0; i < *nperm; ++i) {
        dcopy_(n, &zero, &c__0, &vec[i * *nmvec], &c__1);
    }
    m = *nperm % *nblock;
    if (m == 0) {
        goto L150;
    }
    (*iovect)(n, &m, p, &m, &c__1);
    for (i = 0; i < m; ++i) {
        for (j = 0; j < *nperm; ++j) {
            daxpy_(n, &hv[i + j * *nperm], &p[i * *n], &c__1, &vec[j * *nmvec], &c__1);
        }
    }
    if (*nperm < *nblock) {
        goto L190;
    }
L150:
    m += *nblock;
    for (i = m; *nblock < 0 ? i >= *nperm : i <= *nperm; i += *nblock) {
        (*iovect)(n, nblock, p, &i, &c__1);
        for (j = 0; j < *nblock; ++j) {
            l = i - *nblock + j;
            for (k = 0; k < *nperm; ++k) {
                daxpy_(n, &hv[l + k * *nperm], &p[j * *n], &c__1, &vec[k * *nmvec], &c__1);
            }
        }
    }

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

/* THIS SECTION COMPUTES THE RAYLEIGH QUOTIENTS (IN VAL(*,1)) */
/* AND RESIDUAL NORMS (IN VAL(*,2)) OF THE EIGENVECTORS. */

L190:
    if (! (*small)) {
        *delta = -(*delta);
    }
    m = *nperm % *nblock;
    if (m == 0) {
        goto L220;
    }
    for (i = 0; i < m; ++i) {
        dcopy_(n, &vec[i * *nmvec], &c__1, &p[i * *n], &c__1);
    }
    (*op)(n, &m, p, q);
    ++(*nop);
    for (i = 0; i < m; ++i) {
        val[i] = ddot_(n, &p[i * *n], &c__1, &q[i * *n], &c__1);
        d__1 = -val[i];
        daxpy_(n, &d__1, &p[i * *n], &c__1, &q[i * *n], &c__1);
        val[i + *nmval] = dnrm2_(n, &q[i * *n], &c__1);
    }
    if (*nperm < *nblock) {
        goto L260;
    }
L220:
    ++m;
    for (i = m; *nblock < 0 ? i >= *nperm : i <= *nperm; i += *nblock) {
        for (j = 0; j < *nblock; ++j) {
            l = i - 1 + j;
            dcopy_(n, &vec[l * *nmvec], &c__1, &p[j * *n], &c__1);
        }
        (*op)(n, nblock, p, q);
        ++(*nop);
        for (j = 0; j < *nblock; ++j) {
            l = i - 1 + j;
            val[l] = ddot_(n, &p[j * *n], &c__1, &q[j * *n], &c__1);
            d__1 = -val[l];
            daxpy_(n, &d__1, &p[j * *n], &c__1, &q[j * *n], &c__1);
            val[l + *nmval] = dnrm2_(n, &q[j * *n], &c__1);
        }
    }

/* THIS COMPUTES THE ACCURACY ESTIMATES.  FOR CONSISTENCY WITH DILASO */

L260:
    for (i = 0; i < *nperm; ++i) {
        temp = *delta - val[i];
        if (! (*small)) {
            temp = -temp;
        }
        val[i + *nmval * 3] = 0.;
        if (temp > 0.) {
            val[i + *nmval * 3] = val[i + *nmval] / temp;
        }
        val[i + *nmval * 2] = val[i + *nmval * 3] * val[i + *nmval];
    }

} /* dnppla_ */

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

/* Subroutine */ void dortqr_(nz, n, nblock, z, b)
const integer *nz, *n, *nblock;
doublereal *z, *b;
{
    /* System generated locals */
    doublereal d__1;

    /* Local variables */
    static doublereal temp;
    static integer i, k;
    static doublereal sigma;
    static integer length;
    static doublereal tau;

/* THIS SUBROUTINE COMPUTES THE QR FACTORIZATION OF THE N X NBLOCK */
/* MATRIX Z.  Q IS FORMED IN PLACE AND RETURNED IN Z.  R IS */
/* RETURNED IN B. */

/* THIS SECTION REDUCES Z TO TRIANGULAR FORM. */

    for (i = 0; i < *nblock; ++i) {

/* THIS FORMS THE ITH REFLECTION. */

        length = *n - i;
        d__1 = dnrm2_(&length, &z[i + i * *nz], &c__1);
        sigma = d_sign(&d__1, &z[i + i * *nz]);
        b[i + i * *nblock] = -sigma;
        z[i + i * *nz] += sigma;
        tau = sigma * z[i + i * *nz];

/* THIS APPLIES THE ROTATION TO THE REST OF THE COLUMNS. */

        for (k = i+1; k < *nblock; ++k) {
            if (tau != 0.) {
                temp = -ddot_(&length, &z[i + i * *nz], &c__1, &z[i + k * *nz], &c__1) / tau;
                daxpy_(&length, &temp, &z[i + i * *nz], &c__1, &z[i + k * *nz], &c__1);
            }
            b[i + k * *nblock] = z[i + k * *nz];
            z[i + k * *nz] = 0.;
        }
    }

/* THIS ACCUMULATES THE REFLECTIONS IN REVERSE ORDER. */

    for (i = *nblock-1; i >= 0; --i) {

/* THIS RECREATES THE ITH = NBLOCK-M+1)TH REFLECTIO