zlahqr.c

提供矩阵类的函数库

		}
		z__3.r = x.r + y.r, z__3.i = x.i + y.i;
		zladiv_(&z__2, &u, &z__3);
		z__1.r = t.r - z__2.r, z__1.i = t.i - z__2.i;
		t.r = z__1.r, t.i = z__1.i;
/* **   
                Increment op count */
		opst += 20;
/* ** */
	    }
	}

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

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

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

	    i__3 = h___subscr(m, m);
	    h11.r = h__[i__3].r, h11.i = h__[i__3].i;
	    i__3 = h___subscr(m + 1, m + 1);
	    h22.r = h__[i__3].r, h22.i = h__[i__3].i;
	    z__1.r = h11.r - t.r, z__1.i = h11.i - t.i;
	    h11s.r = z__1.r, h11s.i = z__1.i;
	    i__3 = h___subscr(m + 1, m);
	    h21 = h__[i__3].r;
	    s = (d__1 = h11s.r, abs(d__1)) + (d__2 = d_imag(&h11s), abs(d__2))
		     + abs(h21);
	    z__1.r = h11s.r / s, z__1.i = h11s.i / s;
	    h11s.r = z__1.r, h11s.i = z__1.i;
	    h21 /= s;
	    v[0].r = h11s.r, v[0].i = h11s.i;
	    v[1].r = h21, v[1].i = 0.;
	    i__3 = h___subscr(m, m - 1);
	    h10 = h__[i__3].r;
	    tst1 = ((d__1 = h11s.r, abs(d__1)) + (d__2 = d_imag(&h11s), abs(
		    d__2))) * ((d__3 = h11.r, abs(d__3)) + (d__4 = d_imag(&
		    h11), abs(d__4)) + ((d__5 = h22.r, abs(d__5)) + (d__6 = 
		    d_imag(&h22), abs(d__6))));
	    if ((d__1 = h10 * h21, abs(d__1)) <= ulp * tst1) {
		goto L50;
	    }
/* L40: */
	}
	i__2 = h___subscr(l, l);
	h11.r = h__[i__2].r, h11.i = h__[i__2].i;
	i__2 = h___subscr(l + 1, l + 1);
	h22.r = h__[i__2].r, h22.i = h__[i__2].i;
	z__1.r = h11.r - t.r, z__1.i = h11.i - t.i;
	h11s.r = z__1.r, h11s.i = z__1.i;
	i__2 = h___subscr(l + 1, l);
	h21 = h__[i__2].r;
	s = (d__1 = h11s.r, abs(d__1)) + (d__2 = d_imag(&h11s), abs(d__2)) + 
		abs(h21);
	z__1.r = h11s.r / s, z__1.i = h11s.i / s;
	h11s.r = z__1.r, h11s.i = z__1.i;
	h21 /= s;
	v[0].r = h11s.r, v[0].i = h11s.i;
	v[1].r = h21, v[1].i = 0.;
L50:
/* **   
          Increment op count */
	opst += (i__ - m) * 14;
/* **   

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

             V(2) is always real before the call to ZLARFG, and hence   
             after the call T2 ( = T1*V(2) ) is also real. */

	    if (k > m) {
		zcopy_(&c__2, &h___ref(k, k - 1), &c__1, v, &c__1);
	    }
	    zlarfg_(&c__2, v, &v[1], &c__1, &t1);
/* **   
             Increment op count */
	    opst += 38;
/* ** */
	    if (k > m) {
		i__3 = h___subscr(k, k - 1);
		h__[i__3].r = v[0].r, h__[i__3].i = v[0].i;
		i__3 = h___subscr(k + 1, k - 1);
		h__[i__3].r = 0., h__[i__3].i = 0.;
	    }
	    v2.r = v[1].r, v2.i = v[1].i;
	    z__1.r = t1.r * v2.r - t1.i * v2.i, z__1.i = t1.r * v2.i + t1.i * 
		    v2.r;
	    t2 = z__1.r;

/*           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) {
		d_cnjg(&z__3, &t1);
		i__4 = h___subscr(k, j);
		z__2.r = z__3.r * h__[i__4].r - z__3.i * h__[i__4].i, z__2.i =
			 z__3.r * h__[i__4].i + z__3.i * h__[i__4].r;
		i__5 = h___subscr(k + 1, j);
		z__4.r = t2 * h__[i__5].r, z__4.i = t2 * h__[i__5].i;
		z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
		sum.r = z__1.r, sum.i = z__1.i;
		i__4 = h___subscr(k, j);
		i__5 = h___subscr(k, j);
		z__1.r = h__[i__5].r - sum.r, z__1.i = h__[i__5].i - sum.i;
		h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
		i__4 = h___subscr(k + 1, j);
		i__5 = h___subscr(k + 1, j);
		z__2.r = sum.r * v2.r - sum.i * v2.i, z__2.i = sum.r * v2.i + 
			sum.i * v2.r;
		z__1.r = h__[i__5].r - z__2.r, z__1.i = h__[i__5].i - z__2.i;
		h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
/* L60: */
	    }

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

   Computing MIN */
	    i__4 = k + 2;
	    i__3 = min(i__4,i__);
	    for (j = i1; j <= i__3; ++j) {
		i__4 = h___subscr(j, k);
		z__2.r = t1.r * h__[i__4].r - t1.i * h__[i__4].i, z__2.i = 
			t1.r * h__[i__4].i + t1.i * h__[i__4].r;
		i__5 = h___subscr(j, k + 1);
		z__3.r = t2 * h__[i__5].r, z__3.i = t2 * h__[i__5].i;
		z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
		sum.r = z__1.r, sum.i = z__1.i;
		i__4 = h___subscr(j, k);
		i__5 = h___subscr(j, k);
		z__1.r = h__[i__5].r - sum.r, z__1.i = h__[i__5].i - sum.i;
		h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
		i__4 = h___subscr(j, k + 1);
		i__5 = h___subscr(j, k + 1);
		d_cnjg(&z__3, &v2);
		z__2.r = sum.r * z__3.r - sum.i * z__3.i, z__2.i = sum.r * 
			z__3.i + sum.i * z__3.r;
		z__1.r = h__[i__5].r - z__2.r, z__1.i = h__[i__5].i - z__2.i;
		h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
/* L70: */
	    }
/* **   
             Increment op count   
   Computing MIN */
	    i__3 = 2, i__4 = i__ - k;
	    latime_1.ops += (i2 - i1 + 2 + min(i__3,i__4)) * 20;
/* ** */

	    if (*wantz) {

/*              Accumulate transformations in the matrix Z */

		i__3 = *ihiz;
		for (j = *iloz; j <= i__3; ++j) {
		    i__4 = z___subscr(j, k);
		    z__2.r = t1.r * z__[i__4].r - t1.i * z__[i__4].i, z__2.i =
			     t1.r * z__[i__4].i + t1.i * z__[i__4].r;
		    i__5 = z___subscr(j, k + 1);
		    z__3.r = t2 * z__[i__5].r, z__3.i = t2 * z__[i__5].i;
		    z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
		    sum.r = z__1.r, sum.i = z__1.i;
		    i__4 = z___subscr(j, k);
		    i__5 = z___subscr(j, k);
		    z__1.r = z__[i__5].r - sum.r, z__1.i = z__[i__5].i - 
			    sum.i;
		    z__[i__4].r = z__1.r, z__[i__4].i = z__1.i;
		    i__4 = z___subscr(j, k + 1);
		    i__5 = z___subscr(j, k + 1);
		    d_cnjg(&z__3, &v2);
		    z__2.r = sum.r * z__3.r - sum.i * z__3.i, z__2.i = sum.r *
			     z__3.i + sum.i * z__3.r;
		    z__1.r = z__[i__5].r - z__2.r, z__1.i = z__[i__5].i - 
			    z__2.i;
		    z__[i__4].r = z__1.r, z__[i__4].i = z__1.i;
/* L80: */
		}
/* **   
                Increment op count */
		latime_1.ops += nz * 20;
/* ** */
	    }

	    if (k == m && m > l) {

/*              If the QR step was started at row M > L because two   
                consecutive small subdiagonals were found, then extra   
                scaling must be performed to ensure that H(M,M-1) remains   
                real. */

		z__1.r = 1. - t1.r, z__1.i = 0. - t1.i;
		temp.r = z__1.r, temp.i = z__1.i;
		d__1 = z_abs(&temp);
		z__1.r = temp.r / d__1, z__1.i = temp.i / d__1;
		temp.r = z__1.r, temp.i = z__1.i;
		i__3 = h___subscr(m + 1, m);
		i__4 = h___subscr(m + 1, m);
		d_cnjg(&z__2, &temp);
		z__1.r = h__[i__4].r * z__2.r - h__[i__4].i * z__2.i, z__1.i =
			 h__[i__4].r * z__2.i + h__[i__4].i * z__2.r;
		h__[i__3].r = z__1.r, h__[i__3].i = z__1.i;
		if (m + 2 <= i__) {
		    i__3 = h___subscr(m + 2, m + 1);
		    i__4 = h___subscr(m + 2, m + 1);
		    z__1.r = h__[i__4].r * temp.r - h__[i__4].i * temp.i, 
			    z__1.i = h__[i__4].r * temp.i + h__[i__4].i * 
			    temp.r;
		    h__[i__3].r = z__1.r, h__[i__3].i = z__1.i;
		}
		i__3 = i__;
		for (j = m; j <= i__3; ++j) {
		    if (j != m + 1) {
			if (i2 > j) {
			    i__4 = i2 - j;
			    zscal_(&i__4, &temp, &h___ref(j, j + 1), ldh);
			}
			i__4 = j - i1;
			d_cnjg(&z__1, &temp);
			zscal_(&i__4, &z__1, &h___ref(i1, j), &c__1);
/* **   
                      Increment op count */
			opst += (i2 - i1 + 3) * 6;
/* ** */
			if (*wantz) {
			    d_cnjg(&z__1, &temp);
			    zscal_(&nz, &z__1, &z___ref(*iloz, j), &c__1);
/* **   
                         Increment op count */
			    opst += nz * 6;
/* ** */
			}
		    }
/* L90: */
		}
	    }
/* L100: */
	}

/*        Ensure that H(I,I-1) is real. */

	i__2 = h___subscr(i__, i__ - 1);
	temp.r = h__[i__2].r, temp.i = h__[i__2].i;
	if (d_imag(&temp) != 0.) {
	    rtemp = z_abs(&temp);
	    i__2 = h___subscr(i__, i__ - 1);
	    h__[i__2].r = rtemp, h__[i__2].i = 0.;
	    z__1.r = temp.r / rtemp, z__1.i = temp.i / rtemp;
	    temp.r = z__1.r, temp.i = z__1.i;
	    if (i2 > i__) {
		i__2 = i2 - i__;
		d_cnjg(&z__1, &temp);
		zscal_(&i__2, &z__1, &h___ref(i__, i__ + 1), ldh);
	    }
	    i__2 = i__ - i1;
	    zscal_(&i__2, &temp, &h___ref(i1, i__), &c__1);
/* **   
             Increment op count */
	    opst += (i2 - i1 + 1) * 6;
/* ** */
	    if (*wantz) {
		zscal_(&nz, &temp, &z___ref(*iloz, i__), &c__1);
/* **   
                Increment op count */
		opst += nz * 6;
/* ** */
	    }
	}

/* L110: */
    }

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

    *info = i__;
    return 0;

L120:

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

    i__1 = i__;
    i__2 = h___subscr(i__, i__);
    w[i__1].r = h__[i__2].r, w[i__1].i = h__[i__2].i;

/*     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;

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

/*     End of ZLAHQR */

} /* zlahqr_ */

#undef z___ref
#undef z___subscr
#undef h___ref
#undef h___subscr