dgges.c

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

/*       NB refers to the optimal block size for the immediately */
/*       following subroutine, as returned by ILAENV.) */

    minwrk = 1;
    if (*info == 0 && (*lwork >= 1 || lquery)) {
        minwrk = (*n + 1) * 7 + 16;
        maxwrk = (*n + 1) * 7 + *n * ilaenv_(&c__1, "DGEQRF", " ", n, &c__1, n, &c__0) + 16;
        if (ilvsl) {
            i__1 = (*n + 1) * 7 + *n * ilaenv_(&c__1, "DORGQR", " ", n, &c__1, n, &c_n1);
            maxwrk = max(maxwrk,i__1);
        }
        work[0] = (doublereal) maxwrk;
    }

    if (*lwork < minwrk && ! lquery) {
        *info = -19;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("DGGES ", &i__1);
        return;
    } else if (lquery) {
        return;
    }

/*     Quick return if possible */

    if (*n == 0) {
        *sdim = 0;
        return;
    }

/*     Get machine constants */

    eps = dlamch_("P");
    safmin = dlamch_("S");
    safmax = 1. / safmin;
    dlabad_(&safmin, &safmax);
    smlnum = sqrt(safmin) / eps;
    bignum = 1. / smlnum;

/*     Scale A if max element outside range [SMLNUM,BIGNUM] */

    anrm = dlange_("M", n, n, a, lda, work);
    ilascl = FALSE_;
    if (anrm > 0. && anrm < smlnum) {
        anrmto = smlnum;
        ilascl = TRUE_;
    } else if (anrm > bignum) {
        anrmto = bignum;
        ilascl = TRUE_;
    }
    if (ilascl) {
        dlascl_("G", &c__0, &c__0, &anrm, &anrmto, n, n, a, lda, &ierr);
    }

/*     Scale B if max element outside range [SMLNUM,BIGNUM] */

    bnrm = dlange_("M", n, n, b, ldb, work);
    ilbscl = FALSE_;
    if (bnrm > 0. && bnrm < smlnum) {
        bnrmto = smlnum;
        ilbscl = TRUE_;
    } else if (bnrm > bignum) {
        bnrmto = bignum;
        ilbscl = TRUE_;
    }
    if (ilbscl) {
        dlascl_("G", &c__0, &c__0, &bnrm, &bnrmto, n, n, b, ldb, &ierr);
    }

/*     Permute the matrix to make it more nearly triangular */
/*     (Workspace: need 6*N + 2*N space for storing balancing factors) */

    ileft = 0;
    iright = *n;
    iwrk = iright + *n;
    dggbal_("P", n, a, lda, b, ldb, &ilo, &ihi, &work[ileft], &work[iright], &work[iwrk], &ierr);
    --ilo;

/*     Reduce B to triangular form (QR decomposition of B) */
/*     (Workspace: need N, prefer N*NB) */

    irows = ihi - ilo;
    icols = *n - ilo;
    itau = iwrk;
    iwrk = itau + irows;
    i__1 = *lwork - iwrk;
    dgeqrf_(&irows, &icols, &b[ilo + ilo * *ldb], ldb, &work[itau], &work[iwrk], &i__1, &ierr);

/*     Apply the orthogonal transformation to matrix A */
/*     (Workspace: need N, prefer N*NB) */

    i__1 = *lwork - iwrk;
    dormqr_("L", "T", &irows, &icols, &irows, &b[ilo + ilo * *ldb], ldb, &work[itau],
            &a[ilo + ilo * *lda], lda, &work[iwrk], &i__1, &ierr);

/*     Initialize VSL */
/*     (Workspace: need N, prefer N*NB) */

    if (ilvsl) {
        dlaset_("Full", n, n, &c_b33, &c_b34, vsl, ldvsl);
        i__1 = irows - 1;
        dlacpy_("L", &i__1, &i__1, &b[ilo + 1 + ilo * *ldb], ldb, &vsl[ilo + 1 + ilo * *ldvsl], ldvsl);
        i__1 = *lwork - iwrk;
        dorgqr_(&irows, &irows, &irows, &vsl[ilo + ilo * *ldvsl], ldvsl, &work[itau], &work[iwrk], &i__1, &ierr);
    }

/*     Initialize VSR */

    if (ilvsr) {
        dlaset_("Full", n, n, &c_b33, &c_b34, vsr, ldvsr);
    }

/*     Reduce to generalized Hessenberg form */
/*     (Workspace: none needed) */

    ++ilo;
    dgghrd_(jobvsl, jobvsr, n, &ilo, &ihi, a, lda, b, ldb, vsl, ldvsl, vsr, ldvsr, &ierr);

/*     Perform QZ algorithm, computing Schur vectors if desired */
/*     (Workspace: need N) */

    iwrk = itau;
    i__1 = *lwork - iwrk;
    dhgeqz_("S", jobvsl, jobvsr, n, &ilo, &ihi, a, lda, b, ldb, alphar, alphai, beta,
            vsl, ldvsl, vsr, ldvsr, &work[iwrk], &i__1, &ierr);
    if (ierr != 0) {
        if (ierr > 0 && ierr <= *n) {
            *info = ierr;
        } else if (ierr > *n && ierr <= *n << 1) {
            *info = ierr - *n;
        } else {
            *info = *n + 1;
        }
        work[0] = (doublereal) maxwrk;
        return;
    }

/*     Sort eigenvalues ALPHA/BETA if desired */
/*     (Workspace: need 4*N+16 ) */

    *sdim = 0;
    if (wantst) {

/*        Undo scaling on eigenvalues before DELZTGing */

        if (ilascl) {
            dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, alphar, n, &ierr);
            dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, alphai, n, &ierr);
        }
        if (ilbscl) {
            dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, beta, n, &ierr);
        }

/*        Select eigenvalues */

        for (i = 0; i < *n; ++i) {
            bwork[i] = (*delctg)(&alphar[i], &alphai[i], &beta[i]);
        }

        i__1 = *lwork - iwrk;
        dtgsen_(&c__0, &ilvsl, &ilvsr, bwork, n, a, lda, b, ldb, alphar, alphai, beta,
                vsl, ldvsl, vsr, ldvsr, sdim, &pvsl, &pvsr, dif, &work[iwrk], &i__1, idum, &c__1, &ierr);
        if (ierr == 1) {
            *info = *n + 3;
        }
    }

/*     Apply back-permutation to VSL and VSR */
/*     (Workspace: none needed) */

    if (ilvsl) {
        dggbak_("P", "L", n, &ilo, &ihi, &work[ileft], &work[iright], n, vsl, ldvsl, &ierr);
    }

    if (ilvsr) {
        dggbak_("P", "R", n, &ilo, &ihi, &work[ileft], &work[iright], n, vsr, ldvsr, &ierr);
    }

/*     Check if unscaling would cause over/underflow, if so, rescale */
/*     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of */
/*     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I) */

    if (ilascl) {
        for (i = 0; i < *n; ++i) {
            if (alphai[i] != 0.) {
                if (alphar[i] / safmax > anrmto / anrm || safmin / alphar[i] > anrm / anrmto) {
                    work[0] = abs(a[i + i * *lda] / alphar[i]);
                    beta[i] *= work[0];
                    alphar[i] *= work[0];
                    alphai[i] *= work[0];
                } else if (alphai[i] / safmax > anrmto / anrm || safmin / alphai[i] > anrm / anrmto) {
                    work[0] = abs(a[i + (i+1) * *lda] / alphai[i]);
                    beta[i] *= work[0];
                    alphar[i] *= work[0];
                    alphai[i] *= work[0];
                }
            }
        }
    }

    if (ilbscl) {
        for (i = 0; i < *n; ++i) {
            if (alphai[i] != 0.) {
                if (beta[i] / safmax > bnrmto / bnrm || safmin / beta[i] > bnrm / bnrmto) {
                    work[0] = abs(b[i + i * *ldb] / beta[i]);
                    beta[i] *= work[0];
                    alphar[i] *= work[0];
                    alphai[i] *= work[0];
                }
            }
        }
    }

/*     Undo scaling */

    if (ilascl) {
        dlascl_("H", &c__0, &c__0, &anrmto, &anrm, n, n, a, lda, &ierr);
        dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, alphar, n, &ierr);
        dlascl_("G", &c__0, &c__0, &anrmto, &anrm, n, &c__1, alphai, n, &ierr);
    }

    if (ilbscl) {
        dlascl_("U", &c__0, &c__0, &bnrmto, &bnrm, n, n, b, ldb, &ierr);
        dlascl_("G", &c__0, &c__0, &bnrmto, &bnrm, n, &c__1, beta, n, &ierr);
    }

    if (wantst) {

/*        Check if reordering is correct */

        lastsl = TRUE_;
        lst2sl = TRUE_;
        *sdim = 0;
        ip = 0;
        for (i = 0; i < *n; ++i) {
            cursl = (*delctg)(&alphar[i], &alphai[i], &beta[i]);
            if (alphai[i] == 0.) {
                if (cursl) {
                    ++(*sdim);
                }
                ip = 0;
                if (cursl && ! lastsl) {
                    *info = *n + 2;
                }
            } else {
                if (ip == 1) {

/*                 Last eigenvalue of conjugate pair */

                    cursl = cursl || lastsl;
                    lastsl = cursl;
                    if (cursl) {
                        *sdim += 2;
                    }
                    ip = -1;
                    if (cursl && ! lst2sl) {
                        *info = *n + 2;
                    }
                } else {

/*                 First eigenvalue of conjugate pair */

                    ip = 1;
                }
            }
            lst2sl = lastsl;
            lastsl = cursl;
        }
    }
    work[0] = (doublereal) maxwrk;

} /* dgges_ */