dhgeqz.c

提供矩阵类的函数库

	    ifrstm = ifirst;
	}

/*        Compute single shifts.   

          At this point, IFIRST < ILAST, and the diagonal elements of   
          B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in   
          magnitude) */

	if (iiter / 10 * 10 == iiter) {

/*           Exceptional shift.  Chosen for no particularly good reason.   
             (Single shift only.) */

	    if ((doublereal) maxit * safmin * (d__1 = a_ref(ilast - 1, ilast),
		     abs(d__1)) < (d__2 = b_ref(ilast - 1, ilast - 1), abs(
		    d__2))) {
		eshift += a_ref(ilast - 1, ilast) / b_ref(ilast - 1, ilast - 
			1);
	    } else {
		eshift += 1. / (safmin * (doublereal) maxit);
	    }
	    s1 = 1.;
	    wr = eshift;

/*           ------------------- Begin Timing Code ---------------------- */
	    opst += 4.;
/*           -------------------- End Timing Code ----------------------- */

	} else {

/*           Shifts based on the generalized eigenvalues of the   
             bottom-right 2x2 block of A and B. The first eigenvalue   
             returned by DLAG2 is the Wilkinson shift (AEP p.512), */

	    d__1 = safmin * 100.;
	    dlag2_(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast 
		    - 1), ldb, &d__1, &s1, &s2, &wr, &wr2, &wi);

/* Computing MAX   
   Computing MAX */
	    d__3 = 1., d__4 = abs(wr), d__3 = max(d__3,d__4), d__4 = abs(wi);
	    d__1 = s1, d__2 = safmin * max(d__3,d__4);
	    temp = max(d__1,d__2);

/*           ------------------- Begin Timing Code ---------------------- */
	    opst += 57.;
/*           -------------------- End Timing Code ----------------------- */

	    if (wi != 0.) {
		goto L200;
	    }
	}

/*        Fiddle with shift to avoid overflow */

	temp = min(ascale,1.) * (safmax * .5);
	if (s1 > temp) {
	    scale = temp / s1;
	} else {
	    scale = 1.;
	}

	temp = min(bscale,1.) * (safmax * .5);
	if (abs(wr) > temp) {
/* Computing MIN */
	    d__1 = scale, d__2 = temp / abs(wr);
	    scale = min(d__1,d__2);
	}
	s1 = scale * s1;
	wr = scale * wr;

/*        Now check for two consecutive small subdiagonals. */

	i__2 = ifirst + 1;
	for (j = ilast - 1; j >= i__2; --j) {
	    istart = j;
	    temp = (d__1 = s1 * a_ref(j, j - 1), abs(d__1));
	    temp2 = (d__1 = s1 * a_ref(j, j) - wr * b_ref(j, j), abs(d__1));
	    tempr = max(temp,temp2);
	    if (tempr < 1. && tempr != 0.) {
		temp /= tempr;
		temp2 /= tempr;
	    }
	    if ((d__1 = ascale * a_ref(j + 1, j) * temp, abs(d__1)) <= ascale 
		    * atol * temp2) {
		goto L130;
	    }
/* L120: */
	}

	istart = ifirst;
L130:

/*        Do an implicit single-shift QZ sweep.   

          Initial Q */

	temp = s1 * a_ref(istart, istart) - wr * b_ref(istart, istart);
	temp2 = s1 * a_ref(istart + 1, istart);
	dlartg_(&temp, &temp2, &c__, &s, &tempr);

/*        Sweep */

	i__2 = ilast - 1;
	for (j = istart; j <= i__2; ++j) {
	    if (j > istart) {
		temp = a_ref(j, j - 1);
		dlartg_(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 
			1));
		a_ref(j + 1, j - 1) = 0.;
	    }

	    i__3 = ilastm;
	    for (jc = j; jc <= i__3; ++jc) {
		temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
		a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
		a_ref(j, jc) = temp;
		temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
		b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
		b_ref(j, jc) = temp2;
/* L140: */
	    }
	    if (ilq) {
		i__3 = *n;
		for (jr = 1; jr <= i__3; ++jr) {
		    temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
		    q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j 
			    + 1);
		    q_ref(jr, j) = temp;
/* L150: */
		}
	    }

	    temp = b_ref(j + 1, j + 1);
	    dlartg_(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
	    b_ref(j + 1, j) = 0.;

/* Computing MIN */
	    i__4 = j + 2;
	    i__3 = min(i__4,ilast);
	    for (jr = ifrstm; jr <= i__3; ++jr) {
		temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
		a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
		a_ref(jr, j + 1) = temp;
/* L160: */
	    }
	    i__3 = j;
	    for (jr = ifrstm; jr <= i__3; ++jr) {
		temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
		b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
		b_ref(jr, j + 1) = temp;
/* L170: */
	    }
	    if (ilz) {
		i__3 = *n;
		for (jr = 1; jr <= i__3; ++jr) {
		    temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
		    z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
			    jr, j);
		    z___ref(jr, j + 1) = temp;
/* L180: */
		}
	    }
/* L190: */
	}

/*        --------------------- Begin Timing Code ----------------------- */
	opst += (doublereal) ((ilast - istart) * ((ilastm - ifrstm) * 12 + 58 
		+ (nq + nz) * 6) + 6);
/*        ---------------------- End Timing Code ------------------------ */

	goto L350;

/*        Use Francis double-shift   

          Note: the Francis double-shift should work with real shifts,   
                but only if the block is at least 3x3.   
                This code may break if this point is reached with   
                a 2x2 block with real eigenvalues. */

L200:
	if (ifirst + 1 == ilast) {

/*           Special case -- 2x2 block with complex eigenvectors   

             Step 1: Standardize, that is, rotate so that   

                         ( B11  0  )   
                     B = (         )  with B11 non-negative.   
                         (  0  B22 ) */

	    dlasv2_(&b_ref(ilast - 1, ilast - 1), &b_ref(ilast - 1, ilast), &
		    b_ref(ilast, ilast), &b22, &b11, &sr, &cr, &sl, &cl);

	    if (b11 < 0.) {
		cr = -cr;
		sr = -sr;
		b11 = -b11;
		b22 = -b22;
	    }

	    i__2 = ilastm + 1 - ifirst;
	    drot_(&i__2, &a_ref(ilast - 1, ilast - 1), lda, &a_ref(ilast, 
		    ilast - 1), lda, &cl, &sl);
	    i__2 = ilast + 1 - ifrstm;
	    drot_(&i__2, &a_ref(ifrstm, ilast - 1), &c__1, &a_ref(ifrstm, 
		    ilast), &c__1, &cr, &sr);

	    if (ilast < ilastm) {
		i__2 = ilastm - ilast;
		drot_(&i__2, &b_ref(ilast - 1, ilast + 1), ldb, &b_ref(ilast, 
			ilast + 1), lda, &cl, &sl);
	    }
	    if (ifrstm < ilast - 1) {
		i__2 = ifirst - ifrstm;
		drot_(&i__2, &b_ref(ifrstm, ilast - 1), &c__1, &b_ref(ifrstm, 
			ilast), &c__1, &cr, &sr);
	    }

	    if (ilq) {
		drot_(n, &q_ref(1, ilast - 1), &c__1, &q_ref(1, ilast), &c__1,
			 &cl, &sl);
	    }
	    if (ilz) {
		drot_(n, &z___ref(1, ilast - 1), &c__1, &z___ref(1, ilast), &
			c__1, &cr, &sr);
	    }

	    b_ref(ilast - 1, ilast - 1) = b11;
	    b_ref(ilast - 1, ilast) = 0.;
	    b_ref(ilast, ilast - 1) = 0.;
	    b_ref(ilast, ilast) = b22;

/*           If B22 is negative, negate column ILAST */

	    if (b22 < 0.) {
		i__2 = ilast;
		for (j = ifrstm; j <= i__2; ++j) {
		    a_ref(j, ilast) = -a_ref(j, ilast);
		    b_ref(j, ilast) = -b_ref(j, ilast);
/* L210: */
		}

		if (ilz) {
		    i__2 = *n;
		    for (j = 1; j <= i__2; ++j) {
			z___ref(j, ilast) = -z___ref(j, ilast);
/* L220: */
		    }
		}
	    }

/*           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)   

             Recompute shift */

	    d__1 = safmin * 100.;
	    dlag2_(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast 
		    - 1), ldb, &d__1, &s1, &temp, &wr, &temp2, &wi);

/*           ------------------- Begin Timing Code ---------------------- */
	    opst += (doublereal) ((ilastm + ilast - ifirst - ifrstm) * 12 + 
		    103 + (nq + nz) * 6);
/*           -------------------- End Timing Code -----------------------   

             If standardization has perturbed the shift onto real line,   
             do another (real single-shift) QR step. */

	    if (wi == 0.) {
		goto L350;
	    }
	    s1inv = 1. / s1;

/*           Do EISPACK (QZVAL) computation of alpha and beta */

	    a11 = a_ref(ilast - 1, ilast - 1);
	    a21 = a_ref(ilast, ilast - 1);
	    a12 = a_ref(ilast - 1, ilast);
	    a22 = a_ref(ilast, ilast);

/*           Compute complex Givens rotation on right   
             (Assume some element of C = (sA - wB) > unfl )   
                              __   
             (sA - wB) ( CZ   -SZ )   
                       ( SZ    CZ ) */

	    c11r = s1 * a11 - wr * b11;
	    c11i = -wi * b11;
	    c12 = s1 * a12;
	    c21 = s1 * a21;
	    c22r = s1 * a22 - wr * b22;
	    c22i = -wi * b22;

	    if (abs(c11r) + abs(c11i) + abs(c12) > abs(c21) + abs(c22r) + abs(
		    c22i)) {
		t = dlapy3_(&c12, &c11r, &c11i);
		cz = c12 / t;
		szr = -c11r / t;
		szi = -c11i / t;
	    } else {
		cz = dlapy2_(&c22r, &c22i);
		if (cz <= safmin) {
		    cz = 0.;
		    szr = 1.;
		    szi = 0.;
		} else {
		    tempr = c22r / cz;
		    tempi = c22i / cz;
		    t = dlapy2_(&cz, &c21);
		    cz /= t;
		    szr = -c21 * tempr / t;
		    szi = c21 * tempi / t;
		}
	    }

/*           Compute Givens rotation on left   

             (  CQ   SQ )   
             (  __      )  A or B   
             ( -SQ   CQ ) */

	    an = abs(a11) + abs(a12) + abs(a21) + abs(a22);
	    bn = abs(b11) + abs(b22);
	    wabs = abs(wr) + abs(wi);
	    if (s1 * an > wabs * bn) {
		cq = cz * b11;
		sqr = szr * b22;
		sqi = -szi * b22;
	    } else {
		a1r = cz * a11 + szr * a12;
		a1i = szi * a12;
		a2r = cz * a21 + szr * a22;
		a2i = szi * a22;
		cq = dlapy2_(&a1r, &a1i);
		if (cq <= safmin) {
		    cq = 0.;
		    sqr = 1.;
		    sqi = 0.;
		} else {
		    tempr = a1r / cq;
		    tempi = a1i / cq;
		    sqr = tempr * a2r + tempi * a2i;
		    sqi = tempi * a2r - tempr * a2i;
		}
	    }
	    t = dlapy3_(&cq, &sqr, &sqi);
	    cq /= t;
	    sqr /= t;
	    sqi /= t;

/*           Compute diagonal elements of QBZ */

	    tempr = sqr * szr - sqi * szi;
	    tempi = sqr * szi + sqi * szr;
	    b1r = cq * cz * b11 + tempr * b22;
	    b1i = tempi * b22;
	    b1a = dlapy2_(&b1r, &b1i);
	    b2r = cq * cz * b22 + tempr * b11;
	    b2i = -tempi * b11;
	    b2a = dlapy2_(&b2r, &b2i);

/*           Normalize so beta > 0, and Im( alpha1 ) > 0 */

	    beta[ilast - 1] = b1a;
	    beta[ilast] = b2a;
	    alphar[ilast - 1] = wr * b1a * s1inv;
	    alphai[ilast - 1] = wi * b1a * s1inv;
	    alphar[ilast] = wr * b2a * s1inv;
	    alphai[ilast] = -(wi * b2a) * s1inv;

/*           ------------------- Begin Timing Code ---------------------- */
	    opst += 93.;
/*           -------------------- End Timing Code -----------------------   

             Step 3: Go to next block -- exit if finished. */

	    ilast = ifirst - 1;
	    if (ilast < *ilo) {
		goto L380;
	    }

/*           Reset counters */

	    iiter = 0;
	    eshift = 0.;
	    if (! ilschr) {
		ilastm = ilast;
		if (ifrstm > ilast) {
		    ifrstm = *ilo;