dhgeqz.c

提供矩阵类的函数库

		}
	    }
	    goto L350;
	} else {

/*           Usual case: 3x3 or larger block, using Francis implicit   
                         double-shift   

                                      2   
             Eigenvalue equation is  w  - c w + d = 0,   

                                           -1 2        -1   
             so compute 1st column of  (A B  )  - c A B   + d   
             using the formula in QZIT (from EISPACK)   

             We assume that the block is at least 3x3 */

	    ad11 = ascale * a_ref(ilast - 1, ilast - 1) / (bscale * b_ref(
		    ilast - 1, ilast - 1));
	    ad21 = ascale * a_ref(ilast, ilast - 1) / (bscale * b_ref(ilast - 
		    1, ilast - 1));
	    ad12 = ascale * a_ref(ilast - 1, ilast) / (bscale * b_ref(ilast, 
		    ilast));
	    ad22 = ascale * a_ref(ilast, ilast) / (bscale * b_ref(ilast, 
		    ilast));
	    u12 = b_ref(ilast - 1, ilast) / b_ref(ilast, ilast);
	    ad11l = ascale * a_ref(ifirst, ifirst) / (bscale * b_ref(ifirst, 
		    ifirst));
	    ad21l = ascale * a_ref(ifirst + 1, ifirst) / (bscale * b_ref(
		    ifirst, ifirst));
	    ad12l = ascale * a_ref(ifirst, ifirst + 1) / (bscale * b_ref(
		    ifirst + 1, ifirst + 1));
	    ad22l = ascale * a_ref(ifirst + 1, ifirst + 1) / (bscale * b_ref(
		    ifirst + 1, ifirst + 1));
	    ad32l = ascale * a_ref(ifirst + 2, ifirst + 1) / (bscale * b_ref(
		    ifirst + 1, ifirst + 1));
	    u12l = b_ref(ifirst, ifirst + 1) / b_ref(ifirst + 1, ifirst + 1);

	    v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12 
		    * ad11l + (ad12l - ad11l * u12l) * ad21l;
	    v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 - 
		    ad11l) + ad21 * u12) * ad21l;
	    v[2] = ad32l * ad21l;

	    istart = ifirst;

	    dlarfg_(&c__3, v, &v[1], &c__1, &tau);
	    v[0] = 1.;

/*           Sweep */

	    i__2 = ilast - 2;
	    for (j = istart; j <= i__2; ++j) {

/*              All but last elements: use 3x3 Householder transforms.   

                Zero (j-1)st column of A */

		if (j > istart) {
		    v[0] = a_ref(j, j - 1);
		    v[1] = a_ref(j + 1, j - 1);
		    v[2] = a_ref(j + 2, j - 1);

		    dlarfg_(&c__3, &a_ref(j, j - 1), &v[1], &c__1, &tau);
		    v[0] = 1.;
		    a_ref(j + 1, j - 1) = 0.;
		    a_ref(j + 2, j - 1) = 0.;
		}

		i__3 = ilastm;
		for (jc = j; jc <= i__3; ++jc) {
		    temp = tau * (a_ref(j, jc) + v[1] * a_ref(j + 1, jc) + v[
			    2] * a_ref(j + 2, jc));
		    a_ref(j, jc) = a_ref(j, jc) - temp;
		    a_ref(j + 1, jc) = a_ref(j + 1, jc) - temp * v[1];
		    a_ref(j + 2, jc) = a_ref(j + 2, jc) - temp * v[2];
		    temp2 = tau * (b_ref(j, jc) + v[1] * b_ref(j + 1, jc) + v[
			    2] * b_ref(j + 2, jc));
		    b_ref(j, jc) = b_ref(j, jc) - temp2;
		    b_ref(j + 1, jc) = b_ref(j + 1, jc) - temp2 * v[1];
		    b_ref(j + 2, jc) = b_ref(j + 2, jc) - temp2 * v[2];
/* L230: */
		}
		if (ilq) {
		    i__3 = *n;
		    for (jr = 1; jr <= i__3; ++jr) {
			temp = tau * (q_ref(jr, j) + v[1] * q_ref(jr, j + 1) 
				+ v[2] * q_ref(jr, j + 2));
			q_ref(jr, j) = q_ref(jr, j) - temp;
			q_ref(jr, j + 1) = q_ref(jr, j + 1) - temp * v[1];
			q_ref(jr, j + 2) = q_ref(jr, j + 2) - temp * v[2];
/* L240: */
		    }
		}

/*              Zero j-th column of B (see DLAGBC for details)   

                Swap rows to pivot */

		ilpivt = FALSE_;
/* Computing MAX */
		d__3 = (d__1 = b_ref(j + 1, j + 1), abs(d__1)), d__4 = (d__2 =
			 b_ref(j + 1, j + 2), abs(d__2));
		temp = max(d__3,d__4);
/* Computing MAX */
		d__3 = (d__1 = b_ref(j + 2, j + 1), abs(d__1)), d__4 = (d__2 =
			 b_ref(j + 2, j + 2), abs(d__2));
		temp2 = max(d__3,d__4);
		if (max(temp,temp2) < safmin) {
		    scale = 0.;
		    u1 = 1.;
		    u2 = 0.;
		    goto L250;
		} else if (temp >= temp2) {
		    w11 = b_ref(j + 1, j + 1);
		    w21 = b_ref(j + 2, j + 1);
		    w12 = b_ref(j + 1, j + 2);
		    w22 = b_ref(j + 2, j + 2);
		    u1 = b_ref(j + 1, j);
		    u2 = b_ref(j + 2, j);
		} else {
		    w21 = b_ref(j + 1, j + 1);
		    w11 = b_ref(j + 2, j + 1);
		    w22 = b_ref(j + 1, j + 2);
		    w12 = b_ref(j + 2, j + 2);
		    u2 = b_ref(j + 1, j);
		    u1 = b_ref(j + 2, j);
		}

/*              Swap columns if nec. */

		if (abs(w12) > abs(w11)) {
		    ilpivt = TRUE_;
		    temp = w12;
		    temp2 = w22;
		    w12 = w11;
		    w22 = w21;
		    w11 = temp;
		    w21 = temp2;
		}

/*              LU-factor */

		temp = w21 / w11;
		u2 -= temp * u1;
		w22 -= temp * w12;
		w21 = 0.;

/*              Compute SCALE */

		scale = 1.;
		if (abs(w22) < safmin) {
		    scale = 0.;
		    u2 = 1.;
		    u1 = -w12 / w11;
		    goto L250;
		}
		if (abs(w22) < abs(u2)) {
		    scale = (d__1 = w22 / u2, abs(d__1));
		}
		if (abs(w11) < abs(u1)) {
/* Computing MIN */
		    d__2 = scale, d__3 = (d__1 = w11 / u1, abs(d__1));
		    scale = min(d__2,d__3);
		}

/*              Solve */

		u2 = scale * u2 / w22;
		u1 = (scale * u1 - w12 * u2) / w11;

L250:
		if (ilpivt) {
		    temp = u2;
		    u2 = u1;
		    u1 = temp;
		}

/*              Compute Householder Vector   

   Computing 2nd power */
		d__1 = scale;
/* Computing 2nd power */
		d__2 = u1;
/* Computing 2nd power */
		d__3 = u2;
		t = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
		tau = scale / t + 1.;
		vs = -1. / (scale + t);
		v[0] = 1.;
		v[1] = vs * u1;
		v[2] = vs * u2;

/*              Apply transformations from the right.   

   Computing MIN */
		i__4 = j + 3;
		i__3 = min(i__4,ilast);
		for (jr = ifrstm; jr <= i__3; ++jr) {
		    temp = tau * (a_ref(jr, j) + v[1] * a_ref(jr, j + 1) + v[
			    2] * a_ref(jr, j + 2));
		    a_ref(jr, j) = a_ref(jr, j) - temp;
		    a_ref(jr, j + 1) = a_ref(jr, j + 1) - temp * v[1];
		    a_ref(jr, j + 2) = a_ref(jr, j + 2) - temp * v[2];
/* L260: */
		}
		i__3 = j + 2;
		for (jr = ifrstm; jr <= i__3; ++jr) {
		    temp = tau * (b_ref(jr, j) + v[1] * b_ref(jr, j + 1) + v[
			    2] * b_ref(jr, j + 2));
		    b_ref(jr, j) = b_ref(jr, j) - temp;
		    b_ref(jr, j + 1) = b_ref(jr, j + 1) - temp * v[1];
		    b_ref(jr, j + 2) = b_ref(jr, j + 2) - temp * v[2];
/* L270: */
		}
		if (ilz) {
		    i__3 = *n;
		    for (jr = 1; jr <= i__3; ++jr) {
			temp = tau * (z___ref(jr, j) + v[1] * z___ref(jr, j + 
				1) + v[2] * z___ref(jr, j + 2));
			z___ref(jr, j) = z___ref(jr, j) - temp;
			z___ref(jr, j + 1) = z___ref(jr, j + 1) - temp * v[1];
			z___ref(jr, j + 2) = z___ref(jr, j + 2) - temp * v[2];
/* L280: */
		    }
		}
		b_ref(j + 1, j) = 0.;
		b_ref(j + 2, j) = 0.;
/* L290: */
	    }

/*           Last elements: Use Givens rotations   

             Rotations from the left */

	    j = ilast - 1;
	    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__2 = ilastm;
	    for (jc = j; jc <= i__2; ++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;
/* L300: */
	    }
	    if (ilq) {
		i__2 = *n;
		for (jr = 1; jr <= i__2; ++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;
/* L310: */
		}
	    }

/*           Rotations from the right. */

	    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.;

	    i__2 = ilast;
	    for (jr = ifrstm; jr <= i__2; ++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;
/* L320: */
	    }
	    i__2 = ilast - 1;
	    for (jr = ifrstm; jr <= i__2; ++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;
/* L330: */
	    }
	    if (ilz) {
		i__2 = *n;
		for (jr = 1; jr <= i__2; ++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;
/* L340: */
		}
	    }

/*           ------------------- Begin Timing Code ---------------------- */
	    opst += (doublereal) ((ilastm - ifrstm) * 12 + 86 + (nq + nz) * 6)
		     + (doublereal) (ilast - 1 - istart) * (doublereal) ((
		    ilastm - ifrstm) * 20 + 128 + (nq + nz) * 10);
/*           -------------------- End Timing Code -----------------------   

             End of Double-Shift code */

	}

	goto L350;

/*        End of iteration loop */

L350:
/*        --------------------- Begin Timing Code ----------------------- */
	latime_1.ops += opst;
	opst = 0.;
/*        ---------------------- End Timing Code ------------------------   


   L360: */
    }

/*     Drop-through = non-convergence   

   L370:   
       ---------------------- Begin Timing Code ------------------------- */
    latime_1.ops += opst;
    opst = 0.;
/*     ----------------------- End Timing Code -------------------------- */

    *info = ilast;
    goto L420;

/*     Successful completion of all QZ steps */

L380:

/*     Set Eigenvalues 1:ILO-1 */

    i__1 = *ilo - 1;
    for (j = 1; j <= i__1; ++j) {
	if (b_ref(j, j) < 0.) {
	    if (ilschr) {
		i__2 = j;
		for (jr = 1; jr <= i__2; ++jr) {
		    a_ref(jr, j) = -a_ref(jr, j);
		    b_ref(jr, j) = -b_ref(jr, j);
/* L390: */
		}
	    } else {
		a_ref(j, j) = -a_ref(j, j);
		b_ref(j, j) = -b_ref(j, j);
	    }
	    if (ilz) {
		i__2 = *n;
		for (jr = 1; jr <= i__2; ++jr) {
		    z___ref(jr, j) = -z___ref(jr, j);
/* L400: */
		}
	    }
	}
	alphar[j] = a_ref(j, j);
	alphai[j] = 0.;
	beta[j] = b_ref(j, j);
/* L410: */
    }

/*     Normal Termination */

    *info = 0;

/*     Exit (other than argument error) -- return optimal workspace size */

L420:

/*     ---------------------- Begin Timing Code ------------------------- */
    latime_1.ops += opst;
    opst = 0.;
    latime_1.itcnt = (doublereal) jiter;
/*     ----------------------- End Timing Code -------------------------- */

    work[1] = (doublereal) (*n);
    return 0;

/*     End of DHGEQZ */

} /* dhgeqz_ */

#undef z___ref
#undef q_ref
#undef b_ref
#undef a_ref