dtgexc.c

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

        if (a[*ilst + 1 + *ilst * a_dim1] != 0.) {
            nbl = 2;
        }
    }
    if (*ifst == *ilst) {
        return;
    }

    if (*ifst < *ilst) {

/*        Update ILST. */

        if (nbf == 2 && nbl == 1) {
            --(*ilst);
        }
        if (nbf == 1 && nbl == 2) {
            ++(*ilst);
        }

        here = *ifst;

L10:

/*        Swap with next one below. */

        if (nbf == 1 || nbf == 2) {

/*           Current block either 1-by-1 or 2-by-2. */

            nbnext = 1;
            if (here + nbf + 1 <= *n) {
                if (a[here + nbf + 1 + (here + nbf) * a_dim1] != 0.) {
                    nbnext = 2;
                }
            }
            dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb,
                    &q[q_offset], ldq, &z[z_offset], ldz, &here, &nbf, &nbnext, &work[1], lwork, info);
            if (*info != 0) {
                *ilst = here;
                return;
            }
            here += nbnext;

/*           Test if 2-by-2 block breaks into two 1-by-1 blocks. */

            if (nbf == 2) {
                if (a[here + 1 + here * a_dim1] == 0.) {
                    nbf = 3;
                }
            }
        } else {

/*           Current block consists of two 1-by-1 blocks, each of which */
/*           must be swapped individually. */

            nbnext = 1;
            if (here + 3 <= *n) {
                if (a[here + 3 + (here + 2) * a_dim1] != 0.) {
                    nbnext = 2;
                }
            }
            i__1 = here + 1;
            dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb,
                    &q[q_offset], ldq, &z[z_offset], ldz, &i__1, &c__1, &nbnext, &work[1], lwork, info);
            if (*info != 0) {
                *ilst = here;
                return;
            }
            if (nbnext == 1) {

/*              Swap two 1-by-1 blocks. */

                dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb,
                        &q[q_offset], ldq, &z[z_offset], ldz, &here, &c__1,
                        &c__1, &work[1], lwork, info);
                if (*info != 0) {
                    *ilst = here;
                    return;
                }
                ++here;

            } else {

/*              Recompute NBNEXT in case of 2-by-2 split. */

                if (a[here + 2 + (here + 1) * a_dim1] == 0.) {
                    nbnext = 1;
                }
                if (nbnext == 2) {

/*                 2-by-2 block did not split. */

                    dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                            &z[z_offset], ldz, &here, &c__1, &nbnext, &work[1], lwork, info);
                    if (*info != 0) {
                        *ilst = here;
                        return;
                    }
                    here += 2;
                } else {

/*                 2-by-2 block did split. */

                    dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                            &z[z_offset], ldz, &here, &c__1, &c__1, &work[1], lwork, info);
                    if (*info != 0) {
                        *ilst = here;
                        return;
                    }
                    ++here;
                    dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                            &z[z_offset], ldz, &here, &c__1, &c__1, &work[1], lwork, info);
                    if (*info != 0) {
                        *ilst = here;
                        return;
                    }
                    ++here;
                }
            }
        }
        if (here < *ilst) {
            goto L10;
        }
    } else {
        here = *ifst;

L20:

/*        Swap with next one below. */

        if (nbf == 1 || nbf == 2) {

/*           Current block either 1-by-1 or 2-by-2. */

            nbnext = 1;
            if (here >= 3) {
                if (a[here - 1 + (here - 2) * a_dim1] != 0.) {
                    nbnext = 2;
                }
            }
            i__1 = here - nbnext;
            dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                    &z[z_offset], ldz, &i__1, &nbnext, &nbf, &work[1], lwork, info);
            if (*info != 0) {
                *ilst = here;
                return;
            }
            here -= nbnext;

/*           Test if 2-by-2 block breaks into two 1-by-1 blocks. */

            if (nbf == 2) {
                if (a[here + 1 + here * a_dim1] == 0.) {
                    nbf = 3;
                }
            }

        } else {

/*           Current block consists of two 1-by-1 blocks, each of which */
/*           must be swapped individually. */

            nbnext = 1;
            if (here >= 3) {
                if (a[here - 1 + (here - 2) * a_dim1] != 0.) {
                    nbnext = 2;
                }
            }
            i__1 = here - nbnext;
            dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                    &z[z_offset], ldz, &i__1, &nbnext, &c__1, &work[1], lwork, info);
            if (*info != 0) {
                *ilst = here;
                return;
            }
            if (nbnext == 1) {

/*              Swap two 1-by-1 blocks. */

                dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                        &z[z_offset], ldz, &here, &nbnext, &c__1, &work[1], lwork, info);
                if (*info != 0) {
                    *ilst = here;
                    return;
                }
                --here;
            } else {

/*             Recompute NBNEXT in case of 2-by-2 split. */

                if (a[here + (here - 1) * a_dim1] == 0.) {
                    nbnext = 1;
                }
                if (nbnext == 2) {

/*                 2-by-2 block did not split. */

                    i__1 = here - 1;
                    dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                            &z[z_offset], ldz, &i__1, &c__2, &c__1, &work[1], lwork, info);
                    if (*info != 0) {
                        *ilst = here;
                        return;
                    }
                    here += -2;
                } else {

/*                 2-by-2 block did split. */

                    dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                            &z[z_offset], ldz, &here, &c__1, &c__1, &work[1], lwork, info);
                    if (*info != 0) {
                        *ilst = here;
                        return;
                    }
                    --here;
                    dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq,
                            &z[z_offset], ldz, &here, &c__1, &c__1, &work[1], lwork, info);
                    if (*info != 0) {
                        *ilst = here;
                        return;
                    }
                    --here;
                }
            }
        }
        if (here > *ilst) {
            goto L20;
        }
    }
    *ilst = here;
    work[1] = (doublereal) lwmin;
} /* dtgexc_ */