dhgeqz.c

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

        ilz = TRUE_;
        icompz = 2;
    } else if (lsame_(compz, "I")) {
        ilz = TRUE_;
        icompz = 3;
    } else {
        icompz = 0;
    }

/*     Check Argument Values */

    *info = 0;
    work[1] = (doublereal) max(1,*n);
    lquery = *lwork == -1;
    if (ischur == 0) {
        *info = -1;
    } else if (icompq == 0) {
        *info = -2;
    } else if (icompz == 0) {
        *info = -3;
    } else if (*n < 0) {
        *info = -4;
    } else if (*ilo < 1) {
        *info = -5;
    } else if (*ihi > *n || *ihi < *ilo - 1) {
        *info = -6;
    } else if (*lda < *n) {
        *info = -8;
    } else if (*ldb < *n) {
        *info = -10;
    } else if (*ldq < 1 || ( ilq && *ldq < *n) ) {
        *info = -15;
    } else if (*ldz < 1 || ( ilz && *ldz < *n) ) {
        *info = -17;
    } else if (*lwork < max(1,*n) && ! lquery) {
        *info = -19;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("DHGEQZ", &i__1);
        return;
    } else if (lquery) {
        return;
    }

/*     Quick return if possible */

    if (*n <= 0) {
        work[1] = 1.;
        return;
    }

/*     Initialize Q and Z */

    if (icompq == 3) {
        dlaset_("Full", n, n, &c_b12, &c_b13, &q[q_offset], ldq);
    }
    if (icompz == 3) {
        dlaset_("Full", n, n, &c_b12, &c_b13, &z[z_offset], ldz);
    }

/*     Machine Constants */

    in = *ihi + 1 - *ilo;
    safmin = dlamch_("S");
    safmax = 1. / safmin;
    ulp = dlamch_("E") * dlamch_("B");
    anorm = dlanhs_("F", &in, &a[*ilo + *ilo * a_dim1], lda, &work[1]);
    bnorm = dlanhs_("F", &in, &b[*ilo + *ilo * b_dim1], ldb, &work[1]);
    atol = max(safmin, ulp*anorm);
    btol = max(safmin, ulp*bnorm);
    ascale = 1. / max(safmin,anorm);
    bscale = 1. / max(safmin,bnorm);

/*     Set Eigenvalues IHI+1:N */

    for (j = *ihi + 1; j <= *n; ++j) {
        if (b[j + j * b_dim1] < 0.) {
            if (ilschr) {
                for (jr = 1; jr <= j; ++jr) {
                    a[jr + j * a_dim1] = -a[jr + j * a_dim1];
                    b[jr + j * b_dim1] = -b[jr + j * b_dim1];
                }
            } else {
                a[j + j * a_dim1] = -a[j + j * a_dim1];
                b[j + j * b_dim1] = -b[j + j * b_dim1];
            }
            if (ilz) {
                for (jr = 1; jr <= *n; ++jr) {
                    z[jr + j * z_dim1] = -z[jr + j * z_dim1];
                }
            }
        }
        alphar[j] = a[j + j * a_dim1];
        alphai[j] = 0.;
        beta[j] = b[j + j * b_dim1];
    }

/*     If IHI < ILO, skip QZ steps */

    if (*ihi < *ilo) {
        goto L380;
    }

/*     MAIN QZ ITERATION LOOP */

/*     Initialize dynamic indices */

/*     Eigenvalues ILAST+1:N have been found. */
/*        Column operations modify rows IFRSTM:whatever. */
/*        Row operations modify columns whatever:ILASTM. */

/*     If only eigenvalues are being computed, then */
/*        IFRSTM is the row of the last splitting row above row ILAST; */
/*        this is always at least ILO. */
/*     IITER counts iterations since the last eigenvalue was found, */
/*        to tell when to use an extraordinary shift. */
/*     MAXIT is the maximum number of QZ sweeps allowed. */

    ilast = *ihi;
    if (ilschr) {
        ifrstm = 1;
        ilastm = *n;
    } else {
        ifrstm = *ilo;
        ilastm = *ihi;
    }
    iiter = 0;
    eshift = 0.;
    maxit = (*ihi - *ilo + 1) * 30;

    for (jiter = 1; jiter <= maxit; ++jiter) {

/*        Split the matrix if possible. */

/*        Two tests: */
/*           1: A(j,j-1)=0  or  j=ILO */
/*           2: B(j,j)=0 */

        if (ilast == *ilo) {

/*           Special case: j=ILAST */

            goto L80;
        } else {
            if (abs(a[ilast + (ilast - 1) * a_dim1]) <= atol) {
                a[ilast + (ilast - 1) * a_dim1] = 0.;
                goto L80;
            }
        }

        if (abs(b[ilast + ilast * b_dim1]) <= btol) {
            b[ilast + ilast * b_dim1] = 0.;
            goto L70;
        }

/*        General case: j<ILAST */

        for (j = ilast - 1; j >= *ilo; --j) {

/*           Test 1: for A(j,j-1)=0 or j=ILO */

            if (j == *ilo) {
                ilazro = TRUE_;
            } else {
                if (abs(a[j + (j - 1) * a_dim1]) <= atol) {
                    a[j + (j - 1) * a_dim1] = 0.;
                    ilazro = TRUE_;
                } else {
                    ilazro = FALSE_;
                }
            }

/*           Test 2: for B(j,j)=0 */

            if (abs(b[j + j * b_dim1]) < btol) {
                b[j + j * b_dim1] = 0.;

/*              Test 1a: Check for 2 consecutive small subdiagonals in A */

                ilazr2 = FALSE_;
                if (! ilazro) {
                    temp = abs(a[j + (j - 1) * a_dim1]);
                    temp2 = abs(a[j + j * a_dim1]);
                    tempr = max(temp,temp2);
                    if (tempr < 1. && tempr != 0.) {
                        temp /= tempr;
                        temp2 /= tempr;
                    }
                    if (temp * (ascale * abs(a[j + 1 + j * a_dim1])) <= temp2 * (ascale * atol)) {
                        ilazr2 = TRUE_;
                    }
                }

/*              If both tests pass (1 & 2), i.e., the leading diagonal */
/*              element of B in the block is zero, split a 1x1 block off */
/*              at the top. (I.e., at the J-th row/column) The leading */
/*              diagonal element of the remainder can also be zero, so */
/*              this may have to be done repeatedly. */

                if (ilazro || ilazr2) {
                    for (jch = j; jch < ilast; ++jch) {
                        temp = a[jch + jch * a_dim1];
                        dlartg_(&temp, &a[jch + 1 + jch * a_dim1], &c, &s, &a[jch + jch * a_dim1]);
                        a[jch + 1 + jch * a_dim1] = 0.;
                        i__1 = ilastm - jch;
                        drot_(&i__1, &a[jch + (jch + 1) * a_dim1], lda, &a[jch + 1 + (jch + 1) * a_dim1], lda, &c, &s);
                        i__1 = ilastm - jch;
                        drot_(&i__1, &b[jch + (jch + 1) * b_dim1], ldb, &b[jch + 1 + (jch + 1) * b_dim1], ldb, &c, &s);
                        if (ilq) {
                            drot_(n, &q[jch * q_dim1 + 1], &c__1, &q[(jch + 1) * q_dim1 + 1], &c__1, &c, &s);
                        }
                        if (ilazr2) {
                            a[jch + (jch - 1) * a_dim1] *= c;
                        }
                        ilazr2 = FALSE_;
                        if (abs(b[jch + 1 + (jch + 1) * b_dim1]) >= btol) {
                            if (jch + 1 >= ilast) {
                                goto L80;
                            } else {
                                ifirst = jch + 1;
                                goto L110;
                            }
                        }
                        b[jch + 1 + (jch + 1) * b_dim1] = 0.;
                    }
                    goto L70;
                } else {

/*                 Only test 2 passed -- chase the zero to B(ILAST,ILAST) */
/*                 Then process as in the case B(ILAST,ILAST)=0 */

                    for (jch = j; jch <= ilast; ++jch) {
                        temp = b[jch + (jch + 1) * b_dim1];
                        dlartg_(&temp, &b[jch + 1 + (jch + 1) * b_dim1], &c, &s, &b[jch + (jch + 1) * b_dim1]);
                        b[jch + 1 + (jch + 1) * b_dim1] = 0.;
                        if (jch < ilastm - 1) {
                            i__1 = ilastm - jch - 1;
                            drot_(&i__1, &b[jch + (jch + 2) * b_dim1], ldb, &b[jch + 1 + (jch + 2) * b_dim1], ldb, &c, &s);
                        }
                        i__1 = ilastm - jch + 2;
                        drot_(&i__1, &a[jch + (jch - 1) * a_dim1], lda, &a[jch + 1 + (jch - 1) * a_dim1], lda, &c, &s);
                        if (ilq) {
                            drot_(n, &q[jch * q_dim1 + 1], &c__1, &q[(jch + 1) * q_dim1 + 1], &c__1, &c, &s);
                        }
                        temp = a[jch + 1 + jch * a_dim1];
                        dlartg_(&temp, &a[jch + 1 + (jch - 1) * a_dim1], &c, &s, &a[jch + 1 + jch * a_dim1]);
                        a[jch + 1 + (jch - 1) * a_dim1] = 0.;
                        i__1 = jch + 1 - ifrstm;
                        drot_(&i__1, &a[ifrstm + jch * a_dim1], &c__1, &a[ifrstm + (jch - 1) * a_dim1], &c__1, &c, &s);
                        i__1 = jch - ifrstm;
                        drot_(&i__1, &b[ifrstm + jch * b_dim1], &c__1, &b[ifrstm + (jch - 1) * b_dim1], &c__1, &c, &s);
                        if (ilz) {
                            drot_(n, &z[jch * z_dim1 + 1], &c__1, &z[(jch - 1) * z_dim1 + 1], &c__1, &c, &s);
                        }
                    }
                    goto L70;
                }
            } else if (ilazro) {

/*              Only test 1 passed -- work on J:ILAST */

                ifirst = j;
                goto L110;
            }
/*           Neither test passed -- try next J */
        }

/*        (Drop-through is "impossible") */

        *info = *n + 1;
        goto L420;

/*        B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a */
/*        1x1 block. */

L70:
        temp = a[ilast + ilast * a_dim1];
        dlartg_(&temp, &a[ilast + (ilast - 1) * a_dim1], &c, &s, &a[ilast + ilast * a_dim1]);
        a[ilast + (ilast - 1) * a_dim1] = 0.;
        i__1 = ilast - ifrstm;
        drot_(&i__1, &a[ifrstm + ilast * a_dim1], &c__1, &a[ifrstm + (ilast - 1) * a_dim1], &c__1, &c, &s);
        i__1 = ilast - ifrstm;
        drot_(&i__1, &b[ifrstm + ilast * b_dim1], &c__1, &b[ifrstm + (ilast - 1) * b_dim1], &c__1, &c, &s);
        if (ilz) {
            drot_(n, &z[ilast * z_dim1 + 1], &c__1, &z[(ilast - 1) * z_dim1 + 1], &c__1, &c, &s);
        }

/*        A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI, */
/*                              and BETA */

L80:
        if (b[ilast + ilast * b_dim1] < 0.) {
            if (ilschr) {
                for (j = ifrstm; j <= ilast; ++j) {
                    a[j + ilast * a_dim1] = -a[j + ilast * a_dim1];
                    b[j + ilast * b_dim1] = -b[j + ilast * b_dim1];
                }
            } else {
                a[ilast + ilast * a_dim1] = -a[ilast + ilast * a_dim1];
                b[ilast + ilast * b_dim1] = -b[ilast + ilast * b_dim1];
            }
            if (ilz) {
                for (j = 1; j <= *n; ++j) {
                    z[j + ilast * z_dim1] = -z[j + ilast * z_dim1];
                }
            }
        }
        alphar[ilast] = a[ilast + ilast * a_dim1];
        alphai[ilast] = 0.;
        beta[ilast] = b[ilast + ilast * b_dim1];

/*        Go to next block -- exit if finished. */

        --ilast;
        if (ilast < *ilo) {
            goto L380;