zqrsl.c

InsightToolkit-1.4.0(有大量的优化算法程序)

/*     g.w. stewart, university of maryland, argonne national lab.      */
/*                                                                      */
/*     zqrsl uses the following functions and subprograms.              */
/*                                                                      */
/*     blas zaxpy,zcopy,zdotc                                           */
/*     fortran min0,mod                                                 */
/*                                                                      */
/************************************************************************/

/*     set info flag. */
    *info = 0;

/*     determine what is to be computed. */

    cqy = *job / 10000 != 0;
    cqty = *job % 10000 != 0;
    cb = *job % 1000 / 100 != 0;
    cr = *job % 100 / 10 != 0;
    cxb = *job % 10 != 0;
    ju = min(*k, *n - 1);

/*     special action when n=1. */

    if (ju == 0) {
        if (cqy) {
            qy[0].r = y[0].r, qy[0].i = y[0].i;
        }
        if (cqty) {
            qty[0].r = y[0].r, qty[0].i = y[0].i;
        }
        if (cxb) {
            xb[0].r = y[0].r, xb[0].i = y[0].i;
        }
        if (cb) {
            if (x[0].r == 0. && x[0].i == 0.) {
                *info = 1;
            }
            else {
                z_div(b, y, x);
            }
        }
        if (cr) {
            rsd[0].r = 0., rsd[0].i = 0.;
        }
        return;
    }

/*        set up to compute qy or qty. */

    if (cqy) {
        zcopy_(n, y, &c__1, qy, &c__1);
    }
    if (cqty) {
        zcopy_(n, y, &c__1, qty, &c__1);
    }

/*           compute qy. */

    if (cqy)
    for (j = ju-1; j >= 0; --j) {
        if (qraux[j].r == 0. && qraux[j].i == 0.) {
            continue; /* next j */
        }
        i__1 = j * *ldx + j; /* index [j,j] */
        temp.r = x[i__1].r, temp.i = x[i__1].i;
        ((doublecomplex*)x)[i__1].r = qraux[j].r, ((doublecomplex*)x)[i__1].i = qraux[j].i; /* temporarily */
        i__2 = *n - j;
        zdotc_(&t, &i__2, &x[i__1], &c__1, &qy[j], &c__1);
        z_div(&t, &t, &x[i__1]);
        t.r = -t.r, t.i = -t.i;
        zaxpy_(&i__2, &t, &x[i__1], &c__1, &qy[j], &c__1);
        ((doublecomplex*)x)[i__1].r = temp.r, ((doublecomplex*)x)[i__1].i = temp.i; /* restore original */
    }

/*           compute ctrans(q)*y. */

    if (cqty)
    for (j = 0; j < ju; ++j) {
        if (qraux[j].r == 0. && qraux[j].i == 0.) {
            continue; /* next j */
        }
        i__1 = j * *ldx + j; /* index [j,j] */
        temp.r = x[i__1].r, temp.i = x[i__1].i;
        ((doublecomplex*)x)[i__1].r = qraux[j].r, ((doublecomplex*)x)[i__1].i = qraux[j].i; /* temporarily */
        i__2 = *n - j;
        zdotc_(&t, &i__2, &x[i__1], &c__1, &qty[j], &c__1);
        z_div(&t, &t, &x[i__1]);
        t.r = -t.r, t.i = -t.i;
        zaxpy_(&i__2, &t, &x[i__1], &c__1, &qty[j], &c__1);
        ((doublecomplex*)x)[i__1].r = temp.r, ((doublecomplex*)x)[i__1].i = temp.i; /* restore original */
    }

/*        set up to compute b, rsd, or xb. */

    if (cb) {
        zcopy_(k, qty, &c__1, b, &c__1);
    }
    if (cxb) {
        zcopy_(k, qty, &c__1, xb, &c__1);
    }
    if (cr && *k < *n) {
        i__2 = *n - *k;
        zcopy_(&i__2, &qty[*k], &c__1, &rsd[*k], &c__1);
    }
    if (cxb && *k < *n)
    for (i = *k; i < *n; ++i) {
        xb[i].r = 0., xb[i].i = 0.;
    }
    if (cr)
    for (i = 0; i < *k; ++i) {
        rsd[i].r = 0., rsd[i].i = 0.;
    }

/*           compute b. */

    if (cb)
    for (j = *k-1; j >= 0; --j) {
        i__1 = j * *ldx + j; /* index [j,j] */
        if (x[i__1].r == 0. && x[i__1].i == 0.) {
            *info = j+1;
            break; /* last j */
        }
        z_div(&b[j], &b[j], &x[i__1]);
        if (j == 0) {
            break; /* last j */
        }
        t.r = -b[j].r, t.i = -b[j].i;
        zaxpy_(&j, &t, &x[j* *ldx], &c__1, b, &c__1);
    }

/*           compute rsd or xb as required. */

    if (cr || cxb)
    for (j = ju-1; j >= 0; --j) {
        if (qraux[j].r == 0. && qraux[j].i == 0.) {
            continue; /* next j */
        }
        i__1 = j * *ldx + j; /* index [j,j] */
        temp.r = x[i__1].r, temp.i = x[i__1].i;
        ((doublecomplex*)x)[i__1].r = qraux[j].r, ((doublecomplex*)x)[i__1].i = qraux[j].i; /* temporarily */
        i__2 = *n - j;
        if (cr) {
            zdotc_(&t, &i__2, &x[i__1], &c__1, &rsd[j], &c__1);
            z_div(&t, &t, &x[i__1]);
            t.r = -t.r, t.i = -t.i;
            zaxpy_(&i__2, &t, &x[i__1], &c__1, &rsd[j], &c__1);
        }
        if (cxb) {
            zdotc_(&t, &i__2, &x[i__1], &c__1, &xb[j], &c__1);
            z_div(&t, &t, &x[i__1]);
            t.r = -t.r, t.i = -t.i;
            zaxpy_(&i__2, &t, &x[i__1], &c__1, &xb[j], &c__1);
        }
        ((doublecomplex*)x)[i__1].r = temp.r, ((doublecomplex*)x)[i__1].i = temp.i; /* restore original */
    }
} /* zqrsl_ */