zlahqr.c

来自「算断裂的」· C语言 代码 · 共 576 行 · 第 1/2 页

#include "f2c.h"

/* Subroutine */ int zlahqr_(logical *wantt, logical *wantz, integer *n, 
	integer *ilo, integer *ihi, doublecomplex *h, integer *ldh, 
	doublecomplex *w, integer *iloz, integer *ihiz, doublecomplex *z, 
	integer *ldz, integer *info)
{
/*  -- 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.   

    ===================================================================== 
  


    
   Parameter adjustments   
       Function Body */
    /* Table of constant values */
    static integer c__2 = 2;
    static integer c__1 = 1;
    
    /* System generated locals */
    integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
    doublereal d__1, d__2, d__3, d__4, d__5, d__6;
    doublecomplex z__1, z__2, z__3, z__4;
    /* Builtin functions */
    double d_imag(doublecomplex *);
    void z_sqrt(doublecomplex *, doublecomplex *), d_cnjg(doublecomplex *, 
	    doublecomplex *);
    /* 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;
    extern /* Subroutine */ int zscal_(integer *, doublecomplex *, 
	    doublecomplex *, integer *);
    static doublereal rtemp;
    static integer i1, i2;
    static doublereal rwork[1];
    static doublecomplex t1;
    static doublereal t2;
    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *);
    static doublecomplex v2;
    extern doublereal dlapy2_(doublereal *, doublereal *);
    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
    static doublereal h10;
    static doublecomplex h11;
    static doublereal h21;
    static doublecomplex h22;
    static integer nh;
    extern doublereal dlamch_(char *);
    static integer nz;
    extern /* Subroutine */ int zlarfg_(integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *);
    extern /* Double Complex */ VOID zladiv_(doublecomplex *, doublecomplex *,
	     doublecomplex *);
    extern doublereal zlanhs_(char *, integer *, doublecomplex *, integer *, 
	    doublereal *);
    static doublereal smlnum;
    static doublecomplex h11s;
    static integer itn, its;
    static doublereal ulp;
    static doublecomplex sum;
    static doublereal tst1;



#define V(I) v[(I)]
#define RWORK(I) rwork[(I)]
#define W(I) w[(I)-1]

#define H(I,J) h[(I)-1 + ((J)-1)* ( *ldh)]
#define Z(I,J) z[(I)-1 + ((J)-1)* ( *ldz)]

    *info = 0;

/*     Quick return if possible */

    if (*n == 0) {
	return 0;
    }
    if (*ilo == *ihi) {
	i__1 = *ilo;
	i__2 = *ilo + *ilo * h_dim1;
	W(*ilo).r = H(*ilo,*ilo).r, W(*ilo).i = H(*ilo,*ilo).i;
	return 0;
    }

    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 = 1;
	i2 = *n;
    }

/*     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;
L10:
    if (i < *ilo) {
	goto L130;
    }

/*     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;
    i__1 = itn;
    for (its = 0; its <= itn; ++its) {

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

	i__2 = l + 1;
	for (k = i; k >= l+1; --k) {
	    i__3 = k - 1 + (k - 1) * h_dim1;
	    i__4 = k + k * h_dim1;
	    tst1 = (d__1 = H(k-1,k-1).r, abs(d__1)) + (d__2 = d_imag(&H(k-1,k-1)), abs(d__2)) + ((d__3 = H(k,k).r, abs(
		    d__3)) + (d__4 = d_imag(&H(k,k)), abs(d__4)));
	    if (tst1 == 0.) {
		i__3 = i - l + 1;
		tst1 = zlanhs_("1", &i__3, &H(l,l), ldh, rwork)
			;
	    }
	    i__3 = k + (k - 1) * h_dim1;
/* Computing MAX */
	    d__2 = ulp * tst1;
	    if ((d__1 = H(k,k-1).r, abs(d__1)) <= max(d__2,smlnum)) {
		goto L30;
	    }
/* L20: */
	}
L30:
	l = k;
	if (l > *ilo) {

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

	    i__2 = l + (l - 1) * h_dim1;
	    H(l,l-1).r = 0., H(l,l-1).i = 0.;
	}

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

	if (l >= i) {
	    goto L120;
	}

/*        Now the active submatrix is in rows and columns L to I. If 
  
          eigenvalues only are being computed, only the active submatr
ix   
          need be transformed. */

	if (! (*wantt)) {
	    i1 = l;
	    i2 = i;
	}

	if (its == 10 || its == 20) {

/*           Exceptional shift. */

	    i__2 = i + (i - 1) * h_dim1;
	    i__3 = i - 1 + (i - 2) * h_dim1;
	    d__3 = (d__1 = H(i,i-1).r, abs(d__1)) + (d__2 = H(i-1,i-2).r, abs(
		    d__2));
	    t.r = d__3, t.i = 0.;
	} else {

/*           Wilkinson's shift. */

	    i__2 = i + i * h_dim1;
	    t.r = H(i,i).r, t.i = H(i,i).i;
	    i__2 = i - 1 + i * h_dim1;
	    i__3 = i + (i - 1) * h_dim1;
	    d__1 = H(i,i-1).r;
	    z__1.r = d__1 * H(i-1,i).r, z__1.i = d__1 * H(i-1,i).i;
	    u.r = z__1.r, u.i = z__1.i;
	    if (u.r != 0. || u.i != 0.) {
		i__2 = i - 1 + (i - 1) * h_dim1;
		z__2.r = H(i-1,i-1).r - t.r, z__2.i = H(i-1,i-1).i - t.i;
		z__1.r = z__2.r * .5, z__1.i = z__2.i * .5;