zhseqr.c

提供矩阵类的函数库

	    ftnlen)2);
/* Writing concatenation */
    i__4[0] = 1, a__1[0] = job;
    i__4[1] = 1, a__1[1] = compz;
    s_cat(ch__1, a__1, i__4, &c__2, (ftnlen)2);
    maxb = ilaenv_(&c__8, "ZHSEQR", ch__1, n, ilo, ihi, &c_n1, (ftnlen)6, (
	    ftnlen)2);
    if (ns <= 1 || ns > nh || maxb >= nh) {

/*        Use the standard double-shift algorithm */

	zlahqr_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &w[1], ilo, 
		ihi, &z__[z_offset], ldz, info);
	return 0;
    }
    maxb = max(2,maxb);
/* Computing MIN */
    i__1 = min(ns,maxb);
    ns = min(i__1,15);

/*     Now 1 < NS <= MAXB < NH.   

       Set machine-dependent constants for the stopping criterion.   
       If norm(H) <= sqrt(OVFL), overflow should not occur. */

    unfl = dlamch_("Safe minimum");
    ovfl = 1. / unfl;
    dlabad_(&unfl, &ovfl);
    ulp = dlamch_("Precision");
    smlnum = unfl * (nh / ulp);

/*     ITN is the total number of multiple-shift QR iterations allowed. */

    itn = nh * 30;

/*     The main loop begins here. I is the loop index and decreases from   
       IHI to ILO in steps of at most MAXB. Each iteration of the loop   
       works with the active submatrix in rows and columns L to I.   
       Eigenvalues I+1 to IHI have already converged. Either L = ILO, or   
       H(L,L-1) is negligible so that the matrix splits. */

    i__ = *ihi;
L60:
    if (i__ < *ilo) {
	goto L180;
    }

/*     Perform multiple-shift QR iterations on rows and columns ILO to I   
       until a submatrix of order at most MAXB splits off at the bottom   
       because a subdiagonal element has become negligible. */

    l = *ilo;
    i__1 = itn;
    for (its = 0; its <= i__1; ++its) {

/*        Look for a single small subdiagonal element. */

	i__2 = l + 1;
	for (k = i__; k >= i__2; --k) {
	    i__3 = h___subscr(k - 1, k - 1);
	    i__5 = h___subscr(k, k);
	    tst1 = (d__1 = h__[i__3].r, abs(d__1)) + (d__2 = d_imag(&h___ref(
		    k - 1, k - 1)), abs(d__2)) + ((d__3 = h__[i__5].r, abs(
		    d__3)) + (d__4 = d_imag(&h___ref(k, k)), abs(d__4)));
	    if (tst1 == 0.) {
		i__3 = i__ - l + 1;
		tst1 = zlanhs_("1", &i__3, &h___ref(l, l), ldh, rwork);
/* **   
                Increment op count */
		latime_1.ops += (i__ - l + 1) * 5 * (i__ - l) / 2;
/* ** */
	    }
	    i__3 = h___subscr(k, k - 1);
/* Computing MAX */
	    d__2 = ulp * tst1;
	    if ((d__1 = h__[i__3].r, abs(d__1)) <= max(d__2,smlnum)) {
		goto L80;
	    }
/* L70: */
	}
L80:
	l = k;
/* **   
          Increment op count */
	opst += (i__ - l + 1) * 5;
/* ** */
	if (l > *ilo) {

/*           H(L,L-1) is negligible. */

	    i__2 = h___subscr(l, l - 1);
	    h__[i__2].r = 0., h__[i__2].i = 0.;
	}

/*        Exit from loop if a submatrix of order <= MAXB has split off. */

	if (l >= i__ - maxb + 1) {
	    goto L170;
	}

/*        Now the active submatrix is in rows and columns L to I. If   
          eigenvalues only are being computed, only the active submatrix   
          need be transformed. */

	if (! wantt) {
	    i1 = l;
	    i2 = i__;
	}

	if (its == 20 || its == 30) {

/*           Exceptional shifts. */

	    i__2 = i__;
	    for (ii = i__ - ns + 1; ii <= i__2; ++ii) {
		i__3 = ii;
		i__5 = h___subscr(ii, ii - 1);
		i__6 = h___subscr(ii, ii);
		d__3 = ((d__1 = h__[i__5].r, abs(d__1)) + (d__2 = h__[i__6].r,
			 abs(d__2))) * 1.5;
		w[i__3].r = d__3, w[i__3].i = 0.;
/* L90: */
	    }
/* **   
             Increment op count */
	    opst += ns << 1;
/* ** */
	} else {

/*           Use eigenvalues of trailing submatrix of order NS as shifts. */

	    zlacpy_("Full", &ns, &ns, &h___ref(i__ - ns + 1, i__ - ns + 1), 
		    ldh, s, &c__15);
	    zlahqr_(&c_false, &c_false, &ns, &c__1, &ns, s, &c__15, &w[i__ - 
		    ns + 1], &c__1, &ns, &z__[z_offset], ldz, &ierr);
	    if (ierr > 0) {

/*              If ZLAHQR failed to compute all NS eigenvalues, use the   
                unconverged diagonal elements as the remaining shifts. */

		i__2 = ierr;
		for (ii = 1; ii <= i__2; ++ii) {
		    i__3 = i__ - ns + ii;
		    i__5 = s_subscr(ii, ii);
		    w[i__3].r = s[i__5].r, w[i__3].i = s[i__5].i;
/* L100: */
		}
	    }
	}

/*        Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))   
          where G is the Hessenberg submatrix H(L:I,L:I) and w is   
          the vector of shifts (stored in W). The result is   
          stored in the local array V. */

	v[0].r = 1., v[0].i = 0.;
	i__2 = ns + 1;
	for (ii = 2; ii <= i__2; ++ii) {
	    i__3 = ii - 1;
	    v[i__3].r = 0., v[i__3].i = 0.;
/* L110: */
	}
	nv = 1;
	i__2 = i__;
	for (j = i__ - ns + 1; j <= i__2; ++j) {
	    i__3 = nv + 1;
	    zcopy_(&i__3, v, &c__1, vv, &c__1);
	    i__3 = nv + 1;
	    i__5 = j;
	    z__1.r = -w[i__5].r, z__1.i = -w[i__5].i;
	    zgemv_("No transpose", &i__3, &nv, &c_b2, &h___ref(l, l), ldh, vv,
		     &c__1, &z__1, v, &c__1);
	    ++nv;
/* **   
             Increment op count */
	    opst = opst + (nv << 3) * (*n + 1) + (nv + 1) * 6;
/* **   

             Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,   
             reset it to the unit vector. */

	    itemp = izamax_(&nv, v, &c__1);
/* **   
             Increment op count */
	    opst += nv << 1;
/* ** */
	    i__3 = itemp - 1;
	    rtemp = (d__1 = v[i__3].r, abs(d__1)) + (d__2 = d_imag(&v[itemp - 
		    1]), abs(d__2));
	    if (rtemp == 0.) {
		v[0].r = 1., v[0].i = 0.;
		i__3 = nv;
		for (ii = 2; ii <= i__3; ++ii) {
		    i__5 = ii - 1;
		    v[i__5].r = 0., v[i__5].i = 0.;
/* L120: */
		}
	    } else {
		rtemp = max(rtemp,smlnum);
		d__1 = 1. / rtemp;
		zdscal_(&nv, &d__1, v, &c__1);
/* **   
                Increment op count */
		opst += nv << 1;
/* ** */
	    }
/* L130: */
	}

/*        Multiple-shift QR step */

	i__2 = i__ - 1;
	for (k = l; 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 = ns + 1, i__5 = i__ - k + 1;
	    nr = min(i__3,i__5);
	    if (k > l) {
		zcopy_(&nr, &h___ref(k, k - 1), &c__1, v, &c__1);
	    }
	    zlarfg_(&nr, v, &v[1], &c__1, &tau);
/* **   
             Increment op count */
	    opst = opst + nr * 10 + 12;
/* ** */
	    if (k > l) {
		i__3 = h___subscr(k, k - 1);
		h__[i__3].r = v[0].r, h__[i__3].i = v[0].i;
		i__3 = i__;
		for (ii = k + 1; ii <= i__3; ++ii) {
		    i__5 = h___subscr(ii, k - 1);
		    h__[i__5].r = 0., h__[i__5].i = 0.;
/* L140: */
		}
	    }
	    v[0].r = 1., v[0].i = 0.;

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

	    i__3 = i2 - k + 1;
	    d_cnjg(&z__1, &tau);
	    zlarfx_("Left", &nr, &i__3, v, &z__1, &h___ref(k, k), ldh, &work[
		    1]);

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

   Computing MIN */
	    i__5 = k + nr;
	    i__3 = min(i__5,i__) - i1 + 1;
	    zlarfx_("Right", &i__3, &nr, v, &tau, &h___ref(i1, k), ldh, &work[
		    1]);
/* **   
             Increment op count   
   Computing MIN */
	    i__3 = nr, i__5 = i__ - k;
	    latime_1.ops += ((nr << 2) - 2 << 2) * (i2 - i1 + 2 + min(i__3,
		    i__5));
/* ** */

	    if (wantz) {

/*              Accumulate transformations in the matrix Z */

		zlarfx_("Right", &nh, &nr, v, &tau, &z___ref(*ilo, k), ldz, &
			work[1]);
/* **   
                Increment op count */
		latime_1.ops += ((nr << 2) - 2 << 2) * nh;
/* ** */
	    }
/* L150: */
	}

/*        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.) {
	    d__1 = temp.r;
	    d__2 = d_imag(&temp);
	    rtemp = dlapy2_(&d__1, &d__2);
	    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_(&nh, &temp, &z___ref(*ilo, i__), &c__1);
/* **   
                Increment op count */
		opst += nh * 6;
/* ** */
	    }
	}

/* L160: */
    }

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

    *info = i__;
    return 0;

L170:

/*     A submatrix of order <= MAXB in rows and columns L to I has split   
       off. Use the double-shift QR algorithm to handle it. */

    zlahqr_(&wantt, &wantz, n, &l, &i__, &h__[h_offset], ldh, &w[1], ilo, ihi,
	     &z__[z_offset], ldz, info);
    if (*info > 0) {
	return 0;
    }

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

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

L180:
/* **   
       Compute final op count */
    latime_1.ops += opst;
/* ** */
    i__1 = max(1,*n);
    work[1].r = (doublereal) i__1, work[1].i = 0.;
    return 0;

/*     End of ZHSEQR */

} /* zhseqr_ */

#undef z___ref
#undef z___subscr
#undef s_ref
#undef s_subscr
#undef h___ref
#undef h___subscr