dlahqr.c

著名的LAPACK矩阵计算软件包, 是比较新的版本, 一般用到矩阵分解的朋友也许会用到

	    h33 = h___ref(i__ - 1, i__ - 1);
	    h43h34 = h___ref(i__, i__ - 1) * h___ref(i__ - 1, i__);
	    s = h___ref(i__ - 1, i__ - 2) * h___ref(i__ - 1, i__ - 2);
	    disc = (h33 - h44) * .5;
	    disc = disc * disc + h43h34;
/* **   
             Increment op count */
	    opst += 6;
/* ** */
	    if (disc > 0.) {

/*              Real roots: use Wilkinson's shift twice */

		disc = sqrt(disc);
		ave = (h33 + h44) * .5;
/* **   
                Increment op count */
		opst += 2;
/* ** */
		if (abs(h33) - abs(h44) > 0.) {
		    h33 = h33 * h44 - h43h34;
		    h44 = h33 / (d_sign(&disc, &ave) + ave);
/* **   
                   Increment op count */
		    opst += 4;
/* ** */
		} else {
		    h44 = d_sign(&disc, &ave) + ave;
/* **   
                   Increment op count */
		    opst += 1;
/* ** */
		}
		h33 = h44;
		h43h34 = 0.;
	    }
	}

/*        Look for two consecutive small subdiagonal elements. */

	i__2 = l;
	for (m = i__ - 2; m >= i__2; --m) {

/*           Determine the effect of starting the double-shift QR   
             iteration at row M, and see if this would make H(M,M-1)   
             negligible. */

	    h11 = h___ref(m, m);
	    h22 = h___ref(m + 1, m + 1);
	    h21 = h___ref(m + 1, m);
	    h12 = h___ref(m, m + 1);
	    h44s = h44 - h11;
	    h33s = h33 - h11;
	    v1 = (h33s * h44s - h43h34) / h21 + h12;
	    v2 = h22 - h11 - h33s - h44s;
	    v3 = h___ref(m + 2, m + 1);
	    s = abs(v1) + abs(v2) + abs(v3);
	    v1 /= s;
	    v2 /= s;
	    v3 /= s;
	    v[0] = v1;
	    v[1] = v2;
	    v[2] = v3;
	    if (m == l) {
		goto L50;
	    }
	    h00 = h___ref(m - 1, m - 1);
	    h10 = h___ref(m, m - 1);
	    tst1 = abs(v1) * (abs(h00) + abs(h11) + abs(h22));
	    if (abs(h10) * (abs(v2) + abs(v3)) <= ulp * tst1) {
		goto L50;
	    }
/* L40: */
	}
L50:
/* **   
          Increment op count */
	opst += (i__ - m - 1) * 20;
/* **   

          Double-shift QR step */

	i__2 = i__ - 1;
	for (k = m; k <= i__2; ++k) {

/*           The first iteration of this loop determines a reflection G   
             from the vector V and applies it from left and right to H,   
             thus creating a nonzero bulge below the subdiagonal.   

             Each subsequent iteration determines a reflection G to   
             restore the Hessenberg form in the (K-1)th column, and thus   
             chases the bulge one step toward the bottom of the active   
             submatrix. NR is the order of G.   

   Computing MIN */
	    i__3 = 3, i__4 = i__ - k + 1;
	    nr = min(i__3,i__4);
	    if (k > m) {
		dcopy_(&nr, &h___ref(k, k - 1), &c__1, v, &c__1);
	    }
	    dlarfg_(&nr, v, &v[1], &c__1, &t1);
/* **   
             Increment op count */
	    opst = opst + nr * 3 + 9;
/* ** */
	    if (k > m) {
		h___ref(k, k - 1) = v[0];
		h___ref(k + 1, k - 1) = 0.;
		if (k < i__ - 1) {
		    h___ref(k + 2, k - 1) = 0.;
		}
	    } else if (m > l) {
		h___ref(k, k - 1) = -h___ref(k, k - 1);
	    }
	    v2 = v[1];
	    t2 = t1 * v2;
	    if (nr == 3) {
		v3 = v[2];
		t3 = t1 * v3;

/*              Apply G from the left to transform the rows of the matrix   
                in columns K to I2. */

		i__3 = i2;
		for (j = k; j <= i__3; ++j) {
		    sum = h___ref(k, j) + v2 * h___ref(k + 1, j) + v3 * 
			    h___ref(k + 2, j);
		    h___ref(k, j) = h___ref(k, j) - sum * t1;
		    h___ref(k + 1, j) = h___ref(k + 1, j) - sum * t2;
		    h___ref(k + 2, j) = h___ref(k + 2, j) - sum * t3;
/* L60: */
		}

/*              Apply G from the right to transform the columns of the   
                matrix in rows I1 to min(K+3,I).   

   Computing MIN */
		i__4 = k + 3;
		i__3 = min(i__4,i__);
		for (j = i1; j <= i__3; ++j) {
		    sum = h___ref(j, k) + v2 * h___ref(j, k + 1) + v3 * 
			    h___ref(j, k + 2);
		    h___ref(j, k) = h___ref(j, k) - sum * t1;
		    h___ref(j, k + 1) = h___ref(j, k + 1) - sum * t2;
		    h___ref(j, k + 2) = h___ref(j, k + 2) - sum * t3;
/* L70: */
		}
/* **   
                Increment op count   
   Computing MIN */
		i__3 = 3, i__4 = i__ - k;
		latime_1.ops += (i2 - i1 + 2 + min(i__3,i__4)) * 10;
/* ** */

		if (*wantz) {

/*                 Accumulate transformations in the matrix Z */

		    i__3 = *ihiz;
		    for (j = *iloz; j <= i__3; ++j) {
			sum = z___ref(j, k) + v2 * z___ref(j, k + 1) + v3 * 
				z___ref(j, k + 2);
			z___ref(j, k) = z___ref(j, k) - sum * t1;
			z___ref(j, k + 1) = z___ref(j, k + 1) - sum * t2;
			z___ref(j, k + 2) = z___ref(j, k + 2) - sum * t3;
/* L80: */
		    }
/* **   
                   Increment op count */
		    latime_1.ops += nz * 10;
/* ** */
		}
	    } else if (nr == 2) {

/*              Apply G from the left to transform the rows of the matrix   
                in columns K to I2. */

		i__3 = i2;
		for (j = k; j <= i__3; ++j) {
		    sum = h___ref(k, j) + v2 * h___ref(k + 1, j);
		    h___ref(k, j) = h___ref(k, j) - sum * t1;
		    h___ref(k + 1, j) = h___ref(k + 1, j) - sum * t2;
/* L90: */
		}

/*              Apply G from the right to transform the columns of the   
                matrix in rows I1 to min(K+3,I). */

		i__3 = i__;
		for (j = i1; j <= i__3; ++j) {
		    sum = h___ref(j, k) + v2 * h___ref(j, k + 1);
		    h___ref(j, k) = h___ref(j, k) - sum * t1;
		    h___ref(j, k + 1) = h___ref(j, k + 1) - sum * t2;
/* L100: */
		}
/* **   
                Increment op count */
		latime_1.ops += (i2 - i1 + 3) * 6;
/* ** */

		if (*wantz) {

/*                 Accumulate transformations in the matrix Z */

		    i__3 = *ihiz;
		    for (j = *iloz; j <= i__3; ++j) {
			sum = z___ref(j, k) + v2 * z___ref(j, k + 1);
			z___ref(j, k) = z___ref(j, k) - sum * t1;
			z___ref(j, k + 1) = z___ref(j, k + 1) - sum * t2;
/* L110: */
		    }
/* **   
                   Increment op count */
		    latime_1.ops += nz * 6;
/* ** */
		}
	    }
/* L120: */
	}

/* L130: */
    }

/*     Failure to converge in remaining number of iterations */

    *info = i__;
    return 0;

L140:

    if (l == i__) {

/*        H(I,I-1) is negligible: one eigenvalue has converged. */

	wr[i__] = h___ref(i__, i__);
	wi[i__] = 0.;
    } else if (l == i__ - 1) {

/*        H(I-1,I-2) is negligible: a pair of eigenvalues have converged.   

          Transform the 2-by-2 submatrix to standard Schur form,   
          and compute and store the eigenvalues. */

	dlanv2_(&h___ref(i__ - 1, i__ - 1), &h___ref(i__ - 1, i__), &h___ref(
		i__, i__ - 1), &h___ref(i__, i__), &wr[i__ - 1], &wi[i__ - 1],
		 &wr[i__], &wi[i__], &cs, &sn);

	if (*wantt) {

/*           Apply the transformation to the rest of H. */

	    if (i2 > i__) {
		i__1 = i2 - i__;
		drot_(&i__1, &h___ref(i__ - 1, i__ + 1), ldh, &h___ref(i__, 
			i__ + 1), ldh, &cs, &sn);
	    }
	    i__1 = i__ - i1 - 1;
	    drot_(&i__1, &h___ref(i1, i__ - 1), &c__1, &h___ref(i1, i__), &
		    c__1, &cs, &sn);
/* **   
             Increment op count */
	    latime_1.ops += (i2 - i1 - 1) * 6;
/* ** */
	}
	if (*wantz) {

/*           Apply the transformation to Z. */

	    drot_(&nz, &z___ref(*iloz, i__ - 1), &c__1, &z___ref(*iloz, i__), 
		    &c__1, &cs, &sn);
/* **   
             Increment op count */
	    latime_1.ops += nz * 6;
/* ** */
	}
    }

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

    itn -= its;
    i__ = l - 1;
    goto L10;

L150:
/* **   
       Compute final op count */
    latime_1.ops += opst;
/* ** */
    return 0;

/*     End of DLAHQR */

} /* dlahqr_ */

#undef z___ref
#undef h___ref