clahqr.c

提供矩阵类的函数库

		cladiv_(&q__2, &u, &q__3);
		q__1.r = t.r - q__2.r, q__1.i = t.i - q__2.i;
		t.r = q__1.r, t.i = q__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;
	    q__1.r = h11.r - t.r, q__1.i = h11.i - t.i;
	    h11s.r = q__1.r, h11s.i = q__1.i;
	    i__3 = h___subscr(m + 1, m);
	    h21 = h__[i__3].r;
	    s = (r__1 = h11s.r, dabs(r__1)) + (r__2 = r_imag(&h11s), dabs(
		    r__2)) + dabs(h21);
	    q__1.r = h11s.r / s, q__1.i = h11s.i / s;
	    h11s.r = q__1.r, h11s.i = q__1.i;
	    h21 /= s;
	    v[0].r = h11s.r, v[0].i = h11s.i;
	    v[1].r = h21, v[1].i = 0.f;
	    i__3 = h___subscr(m, m - 1);
	    h10 = h__[i__3].r;
	    tst1 = ((r__1 = h11s.r, dabs(r__1)) + (r__2 = r_imag(&h11s), dabs(
		    r__2))) * ((r__3 = h11.r, dabs(r__3)) + (r__4 = r_imag(&
		    h11), dabs(r__4)) + ((r__5 = h22.r, dabs(r__5)) + (r__6 = 
		    r_imag(&h22), dabs(r__6))));
	    if ((r__1 = h10 * h21, dabs(r__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;
	q__1.r = h11.r - t.r, q__1.i = h11.i - t.i;
	h11s.r = q__1.r, h11s.i = q__1.i;
	i__2 = h___subscr(l + 1, l);
	h21 = h__[i__2].r;
	s = (r__1 = h11s.r, dabs(r__1)) + (r__2 = r_imag(&h11s), dabs(r__2)) 
		+ dabs(h21);
	q__1.r = h11s.r / s, q__1.i = h11s.i / s;
	h11s.r = q__1.r, h11s.i = q__1.i;
	h21 /= s;
	v[0].r = h11s.r, v[0].i = h11s.i;
	v[1].r = h21, v[1].i = 0.f;
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 CLARFG, and hence   
             after the call T2 ( = T1*V(2) ) is also real. */

	    if (k > m) {
		ccopy_(&c__2, &h___ref(k, k - 1), &c__1, v, &c__1);
	    }
	    clarfg_(&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.f, h__[i__3].i = 0.f;
	    }
	    v2.r = v[1].r, v2.i = v[1].i;
	    q__1.r = t1.r * v2.r - t1.i * v2.i, q__1.i = t1.r * v2.i + t1.i * 
		    v2.r;
	    t2 = q__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) {
		r_cnjg(&q__3, &t1);
		i__4 = h___subscr(k, j);
		q__2.r = q__3.r * h__[i__4].r - q__3.i * h__[i__4].i, q__2.i =
			 q__3.r * h__[i__4].i + q__3.i * h__[i__4].r;
		i__5 = h___subscr(k + 1, j);
		q__4.r = t2 * h__[i__5].r, q__4.i = t2 * h__[i__5].i;
		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
		sum.r = q__1.r, sum.i = q__1.i;
		i__4 = h___subscr(k, j);
		i__5 = h___subscr(k, j);
		q__1.r = h__[i__5].r - sum.r, q__1.i = h__[i__5].i - sum.i;
		h__[i__4].r = q__1.r, h__[i__4].i = q__1.i;
		i__4 = h___subscr(k + 1, j);
		i__5 = h___subscr(k + 1, j);
		q__2.r = sum.r * v2.r - sum.i * v2.i, q__2.i = sum.r * v2.i + 
			sum.i * v2.r;
		q__1.r = h__[i__5].r - q__2.r, q__1.i = h__[i__5].i - q__2.i;
		h__[i__4].r = q__1.r, h__[i__4].i = q__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);
		q__2.r = t1.r * h__[i__4].r - t1.i * h__[i__4].i, q__2.i = 
			t1.r * h__[i__4].i + t1.i * h__[i__4].r;
		i__5 = h___subscr(j, k + 1);
		q__3.r = t2 * h__[i__5].r, q__3.i = t2 * h__[i__5].i;
		q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
		sum.r = q__1.r, sum.i = q__1.i;
		i__4 = h___subscr(j, k);
		i__5 = h___subscr(j, k);
		q__1.r = h__[i__5].r - sum.r, q__1.i = h__[i__5].i - sum.i;
		h__[i__4].r = q__1.r, h__[i__4].i = q__1.i;
		i__4 = h___subscr(j, k + 1);
		i__5 = h___subscr(j, k + 1);
		r_cnjg(&q__3, &v2);
		q__2.r = sum.r * q__3.r - sum.i * q__3.i, q__2.i = sum.r * 
			q__3.i + sum.i * q__3.r;
		q__1.r = h__[i__5].r - q__2.r, q__1.i = h__[i__5].i - q__2.i;
		h__[i__4].r = q__1.r, h__[i__4].i = q__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);
		    q__2.r = t1.r * z__[i__4].r - t1.i * z__[i__4].i, q__2.i =
			     t1.r * z__[i__4].i + t1.i * z__[i__4].r;
		    i__5 = z___subscr(j, k + 1);
		    q__3.r = t2 * z__[i__5].r, q__3.i = t2 * z__[i__5].i;
		    q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
		    sum.r = q__1.r, sum.i = q__1.i;
		    i__4 = z___subscr(j, k);
		    i__5 = z___subscr(j, k);
		    q__1.r = z__[i__5].r - sum.r, q__1.i = z__[i__5].i - 
			    sum.i;
		    z__[i__4].r = q__1.r, z__[i__4].i = q__1.i;
		    i__4 = z___subscr(j, k + 1);
		    i__5 = z___subscr(j, k + 1);
		    r_cnjg(&q__3, &v2);
		    q__2.r = sum.r * q__3.r - sum.i * q__3.i, q__2.i = sum.r *
			     q__3.i + sum.i * q__3.r;
		    q__1.r = z__[i__5].r - q__2.r, q__1.i = z__[i__5].i - 
			    q__2.i;
		    z__[i__4].r = q__1.r, z__[i__4].i = q__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. */

		q__1.r = 1.f - t1.r, q__1.i = 0.f - t1.i;
		temp.r = q__1.r, temp.i = q__1.i;
		r__1 = c_abs(&temp);
		q__1.r = temp.r / r__1, q__1.i = temp.i / r__1;
		temp.r = q__1.r, temp.i = q__1.i;
		i__3 = h___subscr(m + 1, m);
		i__4 = h___subscr(m + 1, m);
		r_cnjg(&q__2, &temp);
		q__1.r = h__[i__4].r * q__2.r - h__[i__4].i * q__2.i, q__1.i =
			 h__[i__4].r * q__2.i + h__[i__4].i * q__2.r;
		h__[i__3].r = q__1.r, h__[i__3].i = q__1.i;
		if (m + 2 <= i__) {
		    i__3 = h___subscr(m + 2, m + 1);
		    i__4 = h___subscr(m + 2, m + 1);
		    q__1.r = h__[i__4].r * temp.r - h__[i__4].i * temp.i, 
			    q__1.i = h__[i__4].r * temp.i + h__[i__4].i * 
			    temp.r;
		    h__[i__3].r = q__1.r, h__[i__3].i = q__1.i;
		}
		i__3 = i__;
		for (j = m; j <= i__3; ++j) {
		    if (j != m + 1) {
			if (i2 > j) {
			    i__4 = i2 - j;
			    cscal_(&i__4, &temp, &h___ref(j, j + 1), ldh);
			}
			i__4 = j - i1;
			r_cnjg(&q__1, &temp);
			cscal_(&i__4, &q__1, &h___ref(i1, j), &c__1);
/* **   
                      Increment op count */
			opst += (i2 - i1 + 3) * 6;
/* ** */
			if (*wantz) {
			    r_cnjg(&q__1, &temp);
			    cscal_(&nz, &q__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 (r_imag(&temp) != 0.f) {
	    rtemp = c_abs(&temp);
	    i__2 = h___subscr(i__, i__ - 1);
	    h__[i__2].r = rtemp, h__[i__2].i = 0.f;
	    q__1.r = temp.r / rtemp, q__1.i = temp.i / rtemp;
	    temp.r = q__1.r, temp.i = q__1.i;
	    if (i2 > i__) {
		i__2 = i2 - i__;
		r_cnjg(&q__1, &temp);
		cscal_(&i__2, &q__1, &h___ref(i__, i__ + 1), ldh);
	    }
	    i__2 = i__ - i1;
	    cscal_(&i__2, &temp, &h___ref(i1, i__), &c__1);
/* **   
             Increment op count */
	    opst += (i2 - i1 + 1) * 6;
/* ** */
	    if (*wantz) {
		cscal_(&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 CLAHQR */

} /* clahqr_ */

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