cqrsl.c

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

/*                                                                      */
/*     cqrsl uses the following functions and subprograms.              */
/*                                                                      */
/*     blas caxpy,ccopy,cdotc                                           */
/*     fortran aimag,min0,mod,real                                      */
/*                                                                      */
/************************************************************************/

/*     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.f && x[0].i == 0.f) {
                *info = 1;
            }
            else {
                c_div(b, y, x);
            }
        }
        if (cr) {
            rsd[0].r = 0.f, rsd[0].i = 0.f;
        }
        return;
    }

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

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

/*           compute qy. */

    if (cqy)
    for (j = ju-1; j >= 0; --j) {
        if (qraux[j].r == 0.f && qraux[j].i == 0.f) {
            continue; /* next j */
        }
        i__1 = j * *ldx + j; /* index [j,j] */
        temp.r = x[i__1].r, temp.i = x[i__1].i;
        ((complex*)x)[i__1].r = qraux[j].r, ((complex*)x)[i__1].i = qraux[j].i; /* temporarily */
        i__2 = *n - j;
        cdotc_(&t, &i__2, &x[i__1], &c__1, &qy[j], &c__1);
        c_div(&t, &t, &x[i__1]);
        t.r = -t.r, t.i = -t.i;
        caxpy_(&i__2, &t, &x[i__1], &c__1, &qy[j], &c__1);
        ((complex*)x)[i__1].r = temp.r, ((complex*)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.f && qraux[j].i == 0.f) {
            continue; /* next j */
        }
        i__1 = j * *ldx + j; /* index [j,j] */
        temp.r = x[i__1].r, temp.i = x[i__1].i;
        ((complex*)x)[i__1].r = qraux[j].r, ((complex*)x)[i__1].i = qraux[j].i; /* temporarily */
        i__2 = *n - j;
        cdotc_(&t, &i__2, &x[i__1], &c__1, &qty[j], &c__1);
        c_div(&t, &t, &x[i__1]);
        t.r = -t.r, t.i = -t.i;
        caxpy_(&i__2, &t, &x[i__1], &c__1, &qty[j], &c__1);
        ((complex*)x)[i__1].r = temp.r, ((complex*)x)[i__1].i = temp.i; /* restore original */
    }

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

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

/*           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.f && x[i__1].i == 0.f) {
            *info = j+1;
            break; /* last j */
        }
        c_div(&b[j], &b[j], &x[i__1]);
        if (j == 0) {
            break; /* last j */
        }
        t.r = -b[j].r, t.i = -b[j].i;
        caxpy_(&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.f && qraux[j].i == 0.f) {
            continue; /* next j */
        }
        i__1 = j * *ldx + j; /* index [j,j] */
        temp.r = x[i__1].r, temp.i = x[i__1].i;
        ((complex*)x)[i__1].r = qraux[j].r, ((complex*)x)[i__1].i = qraux[j].i; /* temporarily */
        i__2 = *n - j;
        if (cr) {
            cdotc_(&t, &i__2, &x[i__1], &c__1, &rsd[j], &c__1);
            c_div(&t, &t, &x[i__1]);
            t.r = -t.r, t.i = -t.i;
            caxpy_(&i__2, &t, &x[i__1], &c__1, &rsd[j], &c__1);
        }
        if (cxb) {
            cdotc_(&t, &i__2, &x[i__1], &c__1, &xb[j], &c__1);
            c_div(&t, &t, &x[i__1]);
            t.r = -t.r, t.i = -t.i;
            caxpy_(&i__2, &t, &x[i__1], &c__1, &xb[j], &c__1);
        }
        ((complex*)x)[i__1].r = temp.r, ((complex*)x)[i__1].i = temp.i; /* restore original */
    }
} /* cqrsl_ */