dtgsen.c

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

    pair = FALSE_;
    for (k = 1; k <= *n; ++k) {
        if (pair) {
            pair = FALSE_;
        } else {

            swap = select[k];
            if (k < *n) {
                if (a[k + 1 + k * a_dim1] != 0.) {
                    pair = TRUE_;
                    swap = swap || select[k + 1];
                }
            }

            if (swap) {
                ++ks;

/*              Swap the K-th block to position KS. */
/*              Perform the reordering of diagonal blocks in (A, B) */
/*              by orthogonal transformation matrices and update */
/*              Q and Z accordingly (if requested): */

                kk = k;
                if (k != ks) {
                    dtgexc_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset],
                            ldb, &q[q_offset], ldq, &z[z_offset], ldz, &kk,
                            &ks, &work[1], lwork, &ierr);
                }

                if (ierr > 0) {

/*                 Swap is rejected: exit. */

                    *info = 1;
                    if (wantp) {
                        *pl = 0.;
                        *pr = 0.;
                    }
                    if (wantd) {
                        dif[1] = 0.;
                        dif[2] = 0.;
                    }
                    goto L60;
                }

                if (pair) {
                    ++ks;
                }
            }
        }
    }
    if (wantp) {

/*        Solve generalized Sylvester equation for R and L */
/*        and compute PL and PR. */

        n1 = *m;
        n2 = *n - *m;
        i = n1 + 1;
        ijb = 0;
        dlacpy_("Full", &n1, &n2, &a[i * a_dim1 + 1], lda, &work[1], &n1);
        dlacpy_("Full", &n1, &n2, &b[i * b_dim1 + 1], ldb, &work[n1 * n2 + 1], &n1);
        i__1 = *lwork - (n1 << 1) * n2;
        dtgsyl_("N", &ijb, &n1, &n2, &a[a_offset], lda, &a[i + i * a_dim1], lda,
                &work[1], &n1, &b[b_offset], ldb, &b[i + i * b_dim1], ldb,
                &work[n1 * n2 + 1], &n1, &dscale, &dif[1],
                &work[(n1 * n2 << 1) + 1], &i__1, &iwork[1], &ierr);

/*        Estimate the reciprocal of norms of "projections" onto left */
/*        and right eigenspaces. */

        rdscal = 0.;
        dsum = 1.;
        i__1 = n1 * n2;
        dlassq_(&i__1, &work[1], &c__1, &rdscal, &dsum);
        *pl = rdscal * sqrt(dsum);
        if (*pl == 0.) {
            *pl = 1.;
        } else {
            *pl = dscale / (sqrt(dscale * dscale / *pl + *pl) * sqrt(*pl));
        }
        rdscal = 0.;
        dsum = 1.;
        i__1 = n1 * n2;
        dlassq_(&i__1, &work[n1 * n2 + 1], &c__1, &rdscal, &dsum);
        *pr = rdscal * sqrt(dsum);
        if (*pr == 0.) {
            *pr = 1.;
        } else {
            *pr = dscale / (sqrt(dscale * dscale / *pr + *pr) * sqrt(*pr));
        }
    }

    if (wantd) {

/*        Compute estimates of Difu and Difl. */

        if (wantd1) {
            n1 = *m;
            n2 = *n - *m;
            i = n1 + 1;
            ijb = 3;

/*           Frobenius norm-based Difu-estimate. */

            i__1 = *lwork - (n1 << 1) * n2;
            dtgsyl_("N", &ijb, &n1, &n2, &a[a_offset], lda, &a[i + i * a_dim1], lda,
                    &work[1], &n1, &b[b_offset], ldb, &b[i + i * b_dim1], ldb,
                    &work[n1 * n2 + 1], &n1, &dscale, &dif[1], &work[(n1 << 1) * n2 + 1],
                    &i__1, &iwork[1], &ierr);

/*           Frobenius norm-based Difl-estimate. */

            i__1 = *lwork - (n1 << 1) * n2;
            dtgsyl_("N", &ijb, &n2, &n1, &a[i + i * a_dim1], lda, &a[a_offset], lda,
                    &work[1], &n2, &b[i + i * b_dim1], ldb, &b[b_offset], ldb,
                    &work[n1 * n2 + 1], &n2, &dscale, &dif[2], &work[(n1 << 1) * n2 + 1],
                    &i__1, &iwork[1], &ierr);
        } else {


/*           Compute 1-norm-based estimates of Difu and Difl using */
/*           reversed communication with DLACON. In each step a */
/*           generalized Sylvester equation or a transposed variant */
/*           is solved. */

            kase = 0;
            n1 = *m;
            n2 = *n - *m;
            i = n1 + 1;
            ijb = 0;
            mn2 = (n1 << 1) * n2;

/*           1-norm-based estimate of Difu. */

L40:
            dlacon_(&mn2, &work[mn2 + 1], &work[1], &iwork[1], &dif[1], &kase);
            if (kase != 0) {
                if (kase == 1) {

/*                 Solve generalized Sylvester equation. */

                    i__1 = *lwork - (n1 << 1) * n2;
                    dtgsyl_("N", &ijb, &n1, &n2, &a[a_offset], lda, &a[i + i * a_dim1], lda,
                            &work[1], &n1, &b[b_offset], ldb, &b[i + i * b_dim1], ldb,
                            &work[n1 * n2 + 1], &n1, &dscale, &dif[1],
                            &work[(n1 << 1) * n2 + 1], &i__1, &iwork[1], &ierr);
                } else {

/*                 Solve the transposed variant. */

                    i__1 = *lwork - (n1 << 1) * n2;
                    dtgsyl_("T", &ijb, &n1, &n2, &a[a_offset], lda, &a[i + i * a_dim1], lda,
                            &work[1], &n1, &b[b_offset], ldb, &b[i + i * b_dim1], ldb,
                            &work[n1 * n2 + 1], &n1, &dscale, &dif[1], &work[(n1 << 1) * n2 + 1],
                            &i__1, &iwork[1], &ierr);
                }
                goto L40;
            }
            dif[1] = dscale / dif[1];

/*           1-norm-based estimate of Difl. */

L50:
            dlacon_(&mn2, &work[mn2 + 1], &work[1], &iwork[1], &dif[2], &kase);
            if (kase != 0) {
                if (kase == 1) {

/*                 Solve generalized Sylvester equation. */

                    i__1 = *lwork - (n1 << 1) * n2;
                    dtgsyl_("N", &ijb, &n2, &n1, &a[i + i * a_dim1], lda, &a[a_offset], lda,
                            &work[1], &n2, &b[i + i * b_dim1], ldb, &b[b_offset], ldb,
                            &work[n1 * n2 + 1], &n2, &dscale, &dif[2],
                            &work[(n1 << 1) * n2 + 1], &i__1, &iwork[1], &ierr);
                } else {

/*                 Solve the transposed variant. */

                    i__1 = *lwork - (n1 << 1) * n2;
                    dtgsyl_("T", &ijb, &n2, &n1, &a[i + i * a_dim1], lda, &a[a_offset], lda,
                            &work[1], &n2, &b[i + i * b_dim1], ldb, &b[b_offset], ldb,
                            &work[n1 * n2 + 1], &n2, &dscale, &dif[2],
                            &work[(n1 << 1) * n2 + 1], &i__1, &iwork[1], &ierr);
                }
                goto L50;
            }
            dif[2] = dscale / dif[2];
        }
    }

L60:

/*     Compute generalized eigenvalues of reordered pair (A, B) and */
/*     normalize the generalized Schur form. */

    pair = FALSE_;
    for (k = 1; k <= *n; ++k) {
        if (pair) {
            pair = FALSE_;
        } else {

            if (k < *n) {
                if (a[k + 1 + k * a_dim1] != 0.) {
                    pair = TRUE_;
                }
            }

            if (pair) {

/*             Compute the eigenvalue(s) at position K. */

                work[1] = a[k + k * a_dim1];
                work[2] = a[k + 1 + k * a_dim1];
                work[3] = a[k + (k + 1) * a_dim1];
                work[4] = a[k + 1 + (k + 1) * a_dim1];
                work[5] = b[k + k * b_dim1];
                work[6] = b[k + 1 + k * b_dim1];
                work[7] = b[k + (k + 1) * b_dim1];
                work[8] = b[k + 1 + (k + 1) * b_dim1];
                d__1 = smlnum * eps;
                dlag2_(&work[1], &c__2, &work[5], &c__2, &d__1,
                       &beta[k], &beta[k + 1], &alphar[k], &alphar[k + 1], &alphai[k]);
                alphai[k + 1] = -alphai[k];

            } else {

                if (d_sign(&c_b28, &b[k + k * b_dim1]) < 0.) {

/*                 If B(K,K) is negative, make it positive */

                    for (i = 1; i <= *n; ++i) {
                        a[k + i * a_dim1] = -a[k + i * a_dim1];
                        b[k + i * b_dim1] = -b[k + i * b_dim1];
                        q[i + k * q_dim1] = -q[i + k * q_dim1];
                    }
                }

                alphar[k] = a[k + k * a_dim1];
                alphai[k] = 0.;
                beta[k] = b[k + k * b_dim1];
            }
        }
    }

    work[1] = (doublereal) lwmin;
    iwork[1] = liwmin;

} /* dtgsen_ */