zlahqr.c

来自「算断裂的」· C语言 代码 · 共 576 行 · 第 1/2 页

		x.r = z__1.r, x.i = z__1.i;
		z__3.r = x.r * x.r - x.i * x.i, z__3.i = x.r * x.i + x.i * 
			x.r;
		z__2.r = z__3.r + u.r, z__2.i = z__3.i + u.i;
		z_sqrt(&z__1, &z__2);
		y.r = z__1.r, y.i = z__1.i;
		if (x.r * y.r + d_imag(&x) * d_imag(&y) < 0.) {
		    z__1.r = -y.r, z__1.i = -y.i;
		    y.r = z__1.r, y.i = z__1.i;
		}
		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;
	    }
	}

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

	i__2 = l;
	for (m = i - 1; m >= l; --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 = m + m * h_dim1;
	    h11.r = H(m,m).r, h11.i = H(m,m).i;
	    i__3 = m + 1 + (m + 1) * h_dim1;
	    h22.r = H(m+1,m+1).r, h22.i = H(m+1,m+1).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 = m + 1 + m * h_dim1;
	    h21 = H(m+1,m).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.;
	    if (m == l) {
		goto L50;
	    }
	    i__3 = m + (m - 1) * h_dim1;
	    h10 = H(m,m-1).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: */
	}
L50:

/*        Single-shift QR step */

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

/*           The first iteration of this loop determines a reflect
ion 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 t
o   
             restore the Hessenberg form in the (K-1)th column, an
d thus   
             chases the bulge one step toward the bottom of the ac
tive   
             submatrix.   

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

	    if (k > m) {
		zcopy_(&c__2, &H(k,k-1), &c__1, v, &c__1);
	    }
	    zlarfg_(&c__2, v, &V(1), &c__1, &t1);
	    if (k > m) {
		i__3 = k + (k - 1) * h_dim1;
		H(k,k-1).r = V(0).r, H(k,k-1).i = V(0).i;
		i__3 = k + 1 + (k - 1) * h_dim1;
		H(k+1,k-1).r = 0., H(k+1,k-1).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 ma
trix   
             in columns K to I2. */

	    i__3 = i2;
	    for (j = k; j <= i2; ++j) {
		d_cnjg(&z__3, &t1);
		i__4 = k + j * h_dim1;
		z__2.r = z__3.r * H(k,j).r - z__3.i * H(k,j).i, z__2.i = 
			z__3.r * H(k,j).i + z__3.i * H(k,j).r;
		i__5 = k + 1 + j * h_dim1;
		z__4.r = t2 * H(k+1,j).r, z__4.i = t2 * H(k+1,j).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 = k + j * h_dim1;
		i__5 = k + j * h_dim1;
		z__1.r = H(k,j).r - sum.r, z__1.i = H(k,j).i - sum.i;
		H(k,j).r = z__1.r, H(k,j).i = z__1.i;
		i__4 = k + 1 + j * h_dim1;
		i__5 = k + 1 + j * h_dim1;
		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(k+1,j).r - z__2.r, z__1.i = H(k+1,j).i - z__2.i;
		H(k+1,j).r = z__1.r, H(k+1,j).i = z__1.i;
/* L60: */
	    }

/*           Apply G from the right to transform the columns of th
e   
             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 <= min(k+2,i); ++j) {
		i__4 = j + k * h_dim1;
		z__2.r = t1.r * H(j,k).r - t1.i * H(j,k).i, z__2.i = t1.r * 
			H(j,k).i + t1.i * H(j,k).r;
		i__5 = j + (k + 1) * h_dim1;
		z__3.r = t2 * H(j,k+1).r, z__3.i = t2 * H(j,k+1).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 = j + k * h_dim1;
		i__5 = j + k * h_dim1;
		z__1.r = H(j,k).r - sum.r, z__1.i = H(j,k).i - sum.i;
		H(j,k).r = z__1.r, H(j,k).i = z__1.i;
		i__4 = j + (k + 1) * h_dim1;
		i__5 = j + (k + 1) * h_dim1;
		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(j,k+1).r - z__2.r, z__1.i = H(j,k+1).i - z__2.i;
		H(j,k+1).r = z__1.r, H(j,k+1).i = z__1.i;
/* L70: */
	    }

	    if (*wantz) {

/*              Accumulate transformations in the matrix Z */

		i__3 = *ihiz;
		for (j = *iloz; j <= *ihiz; ++j) {
		    i__4 = j + k * z_dim1;
		    z__2.r = t1.r * Z(j,k).r - t1.i * Z(j,k).i, z__2.i = 
			    t1.r * Z(j,k).i + t1.i * Z(j,k).r;
		    i__5 = j + (k + 1) * z_dim1;
		    z__3.r = t2 * Z(j,k+1).r, z__3.i = t2 * Z(j,k+1).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 = j + k * z_dim1;
		    i__5 = j + k * z_dim1;
		    z__1.r = Z(j,k).r - sum.r, z__1.i = Z(j,k).i - sum.i;
		    Z(j,k).r = z__1.r, Z(j,k).i = z__1.i;
		    i__4 = j + (k + 1) * z_dim1;
		    i__5 = j + (k + 1) * z_dim1;
		    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(j,k+1).r - z__2.r, z__1.i = Z(j,k+1).i - z__2.i;
		    Z(j,k+1).r = z__1.r, Z(j,k+1).i = z__1.i;
/* L80: */
		}
	    }

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

/*              If the QR step was started at row M > L becaus
e two   
                consecutive small subdiagonals were found, the
n 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__2 = temp.r;
		d__3 = d_imag(&temp);
		d__1 = dlapy2_(&d__2, &d__3);
		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 = m + 1 + m * h_dim1;
		i__4 = m + 1 + m * h_dim1;
		d_cnjg(&z__2, &temp);
		z__1.r = H(m+1,m).r * z__2.r - H(m+1,m).i * z__2.i, z__1.i = H(m+1,m).r * z__2.i + H(m+1,m).i * z__2.r;
		H(m+1,m).r = z__1.r, H(m+1,m).i = z__1.i;
		if (m + 2 <= i) {
		    i__3 = m + 2 + (m + 1) * h_dim1;
		    i__4 = m + 2 + (m + 1) * h_dim1;
		    z__1.r = H(m+2,m+1).r * temp.r - H(m+2,m+1).i * temp.i, z__1.i =
			     H(m+2,m+1).r * temp.i + H(m+2,m+1).i * temp.r;
		    H(m+2,m+1).r = z__1.r, H(m+2,m+1).i = z__1.i;
		}
		i__3 = i;
		for (j = m; j <= i; ++j) {
		    if (j != m + 1) {
			if (i2 > j) {
			    i__4 = i2 - j;
			    zscal_(&i__4, &temp, &H(j,j+1), 
				    ldh);
			}
			i__4 = j - i1;
			d_cnjg(&z__1, &temp);
			zscal_(&i__4, &z__1, &H(i1,j), &c__1);
			if (*wantz) {
			    d_cnjg(&z__1, &temp);
			    zscal_(&nz, &z__1, &Z(*iloz,j), &c__1);
			}
		    }
/* L90: */
		}
	    }
/* L100: */
	}

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

	i__2 = i + (i - 1) * h_dim1;
	temp.r = H(i,i-1).r, temp.i = H(i,i-1).i;
	if (d_imag(&temp) != 0.) {
	    d__1 = temp.r;
	    d__2 = d_imag(&temp);
	    rtemp = dlapy2_(&d__1, &d__2);
	    i__2 = i + (i - 1) * h_dim1;
	    H(i,i-1).r = rtemp, H(i,i-1).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(i,i+1), ldh);
	    }
	    i__2 = i - i1;
	    zscal_(&i__2, &temp, &H(i1,i), &c__1);
	    if (*wantz) {
		zscal_(&nz, &temp, &Z(*iloz,i), &c__1);
	    }
	}

/* 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 = i + i * h_dim1;
    W(i).r = H(i,i).r, W(i).i = H(i,i).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:
    return 0;

/*     End of ZLAHQR */

} /* zlahqr_ */