snlaso.c

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

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


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

/* Subroutine */ void slager_(n, nband, nstart, a, tmin, tmax)
const integer *n, *nband, *nstart;
real *a, *tmin, *tmax;
{
    /* Local variables */
    static real 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.f;
        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);
    }
} /* slager_ */


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

/* Subroutine */ void slaran_(n, x)
const integer *n;
real *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] = (real)urand_(&iurand) - .5f;
    }
} /* slaran_ */


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

/* Subroutine */ void smvpc_(nblock, bet, maxj, j, s, number, resnrm, orthcf, rv)
const integer *nblock;
const real *bet;
const integer *maxj, *j;
const real *s;
const integer *number;
real *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] = sdot_(nblock, &s[m + i * *maxj], &c__1, &bet[0], nblock);
        orthcf[i] = abs(rv[0]);
        for (k = 1; k < *nblock; ++k) {
            rv[k] = sdot_(nblock, &s[m + i * *maxj], &c__1, &bet[k], nblock);
            orthcf[i] = min(orthcf[i], abs(rv[k]));
        }
        resnrm[i] = snrm2_(nblock, rv, &c__1);
    }
} /* smvpc_ */


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

/* Subroutine */ void snppla_(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 real*,real*);
/* Subroutine */ void (*iovect) (const integer*,const integer*,real*,const integer*,const integer*);
const integer *n, *nperm, *nmval;
integer *nop;
real *val;
const integer *nmvec;
real *vec;
const integer *nblock;
real *h, *hv, *p, *q, *bound, *d, *delta;
logical *small, *raritz;
real *eps;
{
    /* System generated locals */
    integer i__1;
    real r__1;

    /* Local variables */
    static real hmin, hmax, temp;
    static real zero=0.f;
    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 SLAEIG. */

    i__1 = *nperm * *nperm;
    scopy_(&i__1, &zero, &c__0, h, &c__1);
    temp = -1.f;
    if (*small) {
        temp = 1.f;
    }
    m = *nperm % *nblock;
    if (m == 0) {
        goto L40;
    }
    for (i = 0; i < m; ++i) {
        scopy_(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 * sdot_(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;
            scopy_(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 * sdot_(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];
    slager_(nperm, nperm, &c__1, h, &hmin, &hmax);
    slaeig_(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) {
        scopy_(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) {
            saxpy_(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) {
                saxpy_(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) {
        scopy_(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] = sdot_(n, &p[i * *n], &c__1, &q[i * *n], &c__1);
        r__1 = -val[i];
        saxpy_(n, &r__1, &p[i * *n], &c__1, &q[i * *n], &c__1);
        val[i + *nmval] = snrm2_(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;
            scopy_(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] = sdot_(n, &p[j * *n], &c__1, &q[j * *n], &c__1);
            r__1 = -val[l];
            saxpy_(n, &r__1, &p[j * *n], &c__1, &q[j * *n], &c__1);
            val[l + *nmval] = snrm2_(n, &q[j * *n], &c__1);
        }
    }

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

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

} /* snppla_ */

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

/* Subroutine */ void sortqr_(nz, n, nblock, z, b)
const integer *nz, *n, *nblock;
real *z, *b;
{
    /* System generated locals */
    real r__1;

    /* Local variables */
    static real temp;
    static integer i, k;
    static real sigma;
    static integer length;
    static real 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;
        r__1 = snrm2_(&length, &z[i + i * *nz], &c__1);
        sigma = r_sign(&r__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.f) {
                temp = -sdot_(&length, &z[i + i * *nz], &c__1, &z[i + k * *nz], &c__1) / tau;
                saxpy_(&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.f;
        }
    }

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

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

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

        sigma = -b[i + i * *nblock];
        tau = z[i + i * *nz] * sigma;
        if (tau == 0.f) {
            goto L60;
        }
        length = *n - i;

/* THIS APPLIES IT TO THE LATER COLUMNS. */

        for (k = i+1; k < *nblock; ++k) {
            temp = -sdot_(&length, &z[i + i * *nz], &c__1, &z[i + k * *nz], &c__1) / tau;
            saxpy_(&length, &temp, &z[i + i * *nz], &c__1, &z[i + k * *nz], &c__1);
        }
        r__1 = -1.f / sigma;
        sscal_(&length, &r__1, &z[i + i * *nz], &c__1);
L60:
        z[i + i * *nz] += 1.f;
    }
} /* sortqr_ */


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

/* Subroutine */ void svsort_(num, val, res, iflag, v, nmvec, n, vec)
const integer *num;
real *val, *res;
const integer *iflag;
real *v;
const integer *nmvec, *n;
real *vec;
{
    /* Local variables */
    static real temp;
    static integer kk, k, m;

/*  THIS SUBROUTINE SORTS THE EIGE