zlahqr.c

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

#include "f2c.h"
#include "netlib.h"

/* Modified by Peter Vanroose, June 2001: manual optimisation and clean-up */

/* Table of constant values */
static integer c__2 = 2;
static integer c__1 = 1;

/* Subroutine */ void zlahqr_(wantt, wantz, n, ilo, ihi, h, ldh, w, iloz, ihiz, z, ldz, info)
const logical *wantt, *wantz;
const integer *n, *ilo, *ihi;
doublecomplex *h;
const integer *ldh;
doublecomplex *w;
integer *iloz, *ihiz;
doublecomplex *z;
const integer *ldz;
integer *info;
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1, d__2;
    doublecomplex z__1, z__2;

    /* Local variables */
    static doublereal unfl, ovfl;
    static doublecomplex temp;
    static integer i, j, k, l, m;
    static doublereal s;
    static doublecomplex t, u, v[2], x, y;
    static doublereal rtemp;
    static integer i1, i2;
    static doublereal rwork[1];
    static doublecomplex t1;
    static doublereal t2;
    static doublecomplex v2;
    static doublereal h10;
    static doublecomplex h11;
    static doublereal h21;
    static doublecomplex h22;
    static integer nh;
    static integer nz;
    static doublereal smlnum;
    static doublecomplex h11s;
    static integer itn, its;
    static doublereal ulp;
    static doublecomplex sum;
    static doublereal tst1;


/*  -- LAPACK auxiliary routine (version 2.0) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/*     Courant Institute, Argonne National Lab, and Rice University */
/*     September 30, 1994 */

/*  ===================================================================== */
/*                                                                        */
/*  Purpose                                                               */
/*  =======                                                               */
/*                                                                        */
/*  ZLAHQR is an auxiliary routine called by ZHSEQR to update the         */
/*  eigenvalues and Schur decomposition already computed by ZHSEQR, by    */
/*  dealing with the Hessenberg submatrix in rows and columns ILO to IHI. */
/*                                                                        */
/*  Arguments                                                             */
/*  =========                                                             */
/*                                                                        */
/*  WANTT   (input) LOGICAL                                               */
/*          = .TRUE. : the full Schur form T is required;                 */
/*          = .FALSE.: only eigenvalues are required.                     */
/*                                                                        */
/*  WANTZ   (input) LOGICAL                                               */
/*          = .TRUE. : the matrix of Schur vectors Z is required;         */
/*          = .FALSE.: Schur vectors are not required.                    */
/*                                                                        */
/*  N       (input) INTEGER                                               */
/*          The order of the matrix H.  N >= 0.                           */
/*                                                                        */
/*  ILO     (input) INTEGER                                               */
/*  IHI     (input) INTEGER                                               */
/*          It is assumed that H is already upper triangular in rows and  */
/*          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).  */
/*          ZLAHQR works primarily with the Hessenberg submatrix in rows  */
/*          and columns ILO to IHI, but applies transformations to all of */
/*          H if WANTT is .TRUE..                                         */
/*          1 <= ILO <= max(1,IHI); IHI <= N.                             */
/*                                                                        */
/*  H       (input/output) COMPLEX*16 array, dimension (LDH,N)            */
/*          On entry, the upper Hessenberg matrix H.                      */
/*          On exit, if WANTT is .TRUE., H is upper triangular in rows    */
/*          and columns ILO:IHI, with any 2-by-2 diagonal blocks in       */
/*          standard form. If WANTT is .FALSE., the contents of H are     */
/*          unspecified on exit.                                          */
/*                                                                        */
/*  LDH     (input) INTEGER                                               */
/*          The leading dimension of the array H. LDH >= max(1,N).        */
/*                                                                        */
/*  W       (output) COMPLEX*16 array, dimension (N)                      */
/*          The computed eigenvalues ILO to IHI are stored in the         */
/*          corresponding elements of W. If WANTT is .TRUE., the          */
/*          eigenvalues are stored in the same order as on the diagonal   */
/*          of the Schur form returned in H, with W(i) = H(i,i).          */
/*                                                                        */
/*  ILOZ    (input) INTEGER                                               */
/*  IHIZ    (input) INTEGER                                               */
/*          Specify the rows of Z to which transformations must be        */
/*          applied if WANTZ is .TRUE..                                   */
/*          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.                           */
/*                                                                        */
/*  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)            */
/*          If WANTZ is .TRUE., on entry Z must contain the current       */
/*          matrix Z of transformations accumulated by ZHSEQR, and on     */
/*          exit Z has been updated; transformations are applied only to  */
/*          the submatrix Z(ILOZ:IHIZ,ILO:IHI).                           */
/*          If WANTZ is .FALSE., Z is not referenced.                     */
/*                                                                        */
/*  LDZ     (input) INTEGER                                               */
/*          The leading dimension of the array Z. LDZ >= max(1,N).        */
/*                                                                        */
/*  INFO    (output) INTEGER                                              */
/*          = 0: successful exit                                          */
/*          > 0: if INFO = i, ZLAHQR failed to compute all the            */
/*               eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1)      */
/*               iterations; elements i+1:ihi of W contain those          */
/*               eigenvalues which have been successfully computed.       */
/*                                                                        */
/*  ===================================================================== */

    *info = 0;

/*     Quick return if possible */

    if (*n == 0) {
        return;
    }
    if (*ilo == *ihi) {
        i__1 = *ilo-1 + (*ilo-1) * *ldh;
        w[*ilo-1].r = h[i__1].r, w[*ilo-1].i = h[i__1].i;
        return;
    }

    nh = *ihi - *ilo + 1;
    nz = *ihiz - *iloz + 1;

/*     Set machine-dependent constants for the stopping criterion. */
/*     If norm(H) <= sqrt(OVFL), overflow should not occur. */

    unfl = dlamch_("Safe minimum");
    ovfl = 1. / unfl;
    dlabad_(&unfl, &ovfl);
    ulp = dlamch_("Precision");
    smlnum = unfl * (nh / ulp);

/*     I1 and I2 are the indices of the first row and last column of H */
/*     to which transformations must be applied. If eigenvalues only are */
/*     being computed, I1 and I2 are set inside the main loop. */

    if (*wantt) {
        i1 = 0;
        i2 = *n-1;
    }

/*     ITN is the total number of QR iterations allowed. */

    itn = nh * 30;

/*     The main loop begins here. I is the loop index and decreases from */
/*     IHI to ILO in steps of 1. Each iteration of the loop works */
/*     with the active submatrix in rows and columns L to I. */
/*     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or */
/*     H(L,L-1) is negligible so that the matrix splits. */

    i = *ihi-1;
L10:
    if (i < *ilo-1) {
        return; /* exit from zlahqr_ */
    }

/*     Perform QR iterations on rows and columns ILO to I until a */
/*     submatrix of order 1 splits off at the bottom because a */
/*     subdiagonal element has become negligible. */

    l = *ilo-1;
    for (its = 0; /*fsm*/ 1 /*its <= itn*/; ++its) {

/*        Look for a single small subdiagonal element. */

        for (k = i; k > l; --k) {
            i__1 = k - 1 + (k - 1) * *ldh;
            i__2 = k + k * *ldh;
            tst1 = abs(h[i__1].r) + abs(h[i__1].i) + abs(h[i__2].r) + abs(h[i__2].i);
            if (tst1 == 0.) {
                i__1 = i - l + 1;
                tst1 = zlanhs_("1", &i__1, &h[l + l * *ldh], ldh, rwork);
            }
            if (abs(h[k + (k - 1) * *ldh].r) <= max(ulp * tst1,smlnum)) {
                break;
            }
        }
        l = k;
        if (l > *ilo-1) {

/*           H(L,L-1) is negligible */

            i__1 = l + (l - 1) * *ldh;
            h[i__1].r = h[i__1].i = 0.;
        }

/*        Exit from loop if a submatrix of order 1 has split off. */

        if (l >= i) {
/*         H(I,I-1) is negligible: one eigenvalue has converged. */

            i__1 = i + i * *ldh;
            w[i].r = h[i__1].r, w[i].i = h[i__1].i;

/*         Decrement number of remaining iterations, and return to start of */
/*         the main loop with new value of I. */