dhseqr.c

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

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

/*     I1 and I2 are the indices of the first row and last column of H   
       to which transformations must be applied. If eigenvalues only are   
       being computed, I1 and I2 are set inside the main loop. */

    if (wantt) {
	i1 = 1;
	i2 = *n;
    }

/*     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;
L50:
    l = *ilo;
    if (i__ < *ilo) {
	goto L170;
    }

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

    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) {
	    tst1 = (d__1 = h___ref(k - 1, k - 1), abs(d__1)) + (d__2 = 
		    h___ref(k, k), abs(d__2));
	    if (tst1 == 0.) {
		i__4 = i__ - l + 1;
		tst1 = dlanhs_("1", &i__4, &h___ref(l, l), ldh, &work[1]);
/* **   
                Increment op count */
		latime_1.ops += (i__ - l + 1) * (i__ - l + 2) / 2;
/* ** */
	    }
/* Computing MAX */
	    d__2 = ulp * tst1;
	    if ((d__1 = h___ref(k, k - 1), abs(d__1)) <= max(d__2,smlnum)) {
		goto L70;
	    }
/* L60: */
	}
L70:
	l = k;
/* **   
          Increment op count */
	opst += (i__ - l + 1) * 3;
/* ** */
	if (l > *ilo) {

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

	    h___ref(l, l - 1) = 0.;
	}

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

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

/*        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) {
		wr[ii] = ((d__1 = h___ref(ii, ii - 1), abs(d__1)) + (d__2 = 
			h___ref(ii, ii), abs(d__2))) * 1.5;
		wi[ii] = 0.;
/* L80: */
	    }
/* **   
             Increment op count */
	    opst += ns << 1;
/* ** */
	} else {

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

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

/*              If DLAHQR 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) {
		    wr[i__ - ns + ii] = s_ref(ii, ii);
		    wi[i__ - ns + ii] = 0.;
/* L90: */
		}
	    }
	}

/*        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 WR and WI). The result is   
          stored in the local array V. */

	v[0] = 1.;
	i__2 = ns + 1;
	for (ii = 2; ii <= i__2; ++ii) {
	    v[ii - 1] = 0.;
/* L100: */
	}
	nv = 1;
	i__2 = i__;
	for (j = i__ - ns + 1; j <= i__2; ++j) {
	    if (wi[j] >= 0.) {
		if (wi[j] == 0.) {

/*                 real shift */

		    i__4 = nv + 1;
		    dcopy_(&i__4, v, &c__1, vv, &c__1);
		    i__4 = nv + 1;
		    d__1 = -wr[j];
		    dgemv_("No transpose", &i__4, &nv, &c_b10, &h___ref(l, l),
			     ldh, vv, &c__1, &d__1, v, &c__1);
		    ++nv;
/* **   
                   Increment op count */
		    opst = opst + (nv << 1) * (nv + 1) + nv + 1;
/* ** */
		} else if (wi[j] > 0.) {

/*                 complex conjugate pair of shifts */

		    i__4 = nv + 1;
		    dcopy_(&i__4, v, &c__1, vv, &c__1);
		    i__4 = nv + 1;
		    d__1 = wr[j] * -2.;
		    dgemv_("No transpose", &i__4, &nv, &c_b10, &h___ref(l, l),
			     ldh, v, &c__1, &d__1, vv, &c__1);
		    i__4 = nv + 1;
		    itemp = idamax_(&i__4, vv, &c__1);
/* Computing MAX */
		    d__2 = (d__1 = vv[itemp - 1], abs(d__1));
		    temp = 1. / max(d__2,smlnum);
		    i__4 = nv + 1;
		    dscal_(&i__4, &temp, vv, &c__1);
		    absw = dlapy2_(&wr[j], &wi[j]);
		    temp = temp * absw * absw;
		    i__4 = nv + 2;
		    i__5 = nv + 1;
		    dgemv_("No transpose", &i__4, &i__5, &c_b10, &h___ref(l, 
			    l), ldh, vv, &c__1, &temp, v, &c__1);
		    nv += 2;
/* **   
                   Increment op count   
   Computing 2nd power */
		    i__4 = nv + 1;
		    opst = opst + (i__4 * i__4 << 2) + (nv << 2) + 9;
/* ** */
		}

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

		itemp = idamax_(&nv, v, &c__1);
/* **   
                Increment op count */
		opst += nv;
/* ** */
		temp = (d__1 = v[itemp - 1], abs(d__1));
		if (temp == 0.) {
		    v[0] = 1.;
		    i__4 = nv;
		    for (ii = 2; ii <= i__4; ++ii) {
			v[ii - 1] = 0.;
/* L110: */
		    }
		} else {
		    temp = max(temp,smlnum);
		    d__1 = 1. / temp;
		    dscal_(&nv, &d__1, v, &c__1);
/* **   
                   Increment op count */
		    opst += nv;
/* ** */
		}
	    }
/* L120: */
	}

/*        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__4 = ns + 1, i__5 = i__ - k + 1;
	    nr = min(i__4,i__5);
	    if (k > l) {
		dcopy_(&nr, &h___ref(k, k - 1), &c__1, v, &c__1);
	    }
	    dlarfg_(&nr, v, &v[1], &c__1, &tau);
/* **   
             Increment op count */
	    opst = opst + nr * 3 + 9;
/* ** */
	    if (k > l) {
		h___ref(k, k - 1) = v[0];
		i__4 = i__;
		for (ii = k + 1; ii <= i__4; ++ii) {
		    h___ref(ii, k - 1) = 0.;
/* L130: */
		}
	    }
	    v[0] = 1.;

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

	    i__4 = i2 - k + 1;
	    dlarfx_("Left", &nr, &i__4, v, &tau, &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__4 = min(i__5,i__) - i1 + 1;
	    dlarfx_("Right", &i__4, &nr, v, &tau, &h___ref(i1, k), ldh, &work[
		    1]);
/* **   
             Increment op count   
   Computing MIN */
	    i__4 = nr, i__5 = i__ - k;
	    latime_1.ops += ((nr << 2) - 2) * (i2 - i1 + 2 + min(i__4,i__5));
/* ** */

	    if (wantz) {

/*              Accumulate transformations in the matrix Z */

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

/* L150: */
    }

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

    *info = i__;
    return 0;

L160:

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

    dlahqr_(&wantt, &wantz, n, &l, &i__, &h__[h_offset], ldh, &wr[1], &wi[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 L50;

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

/*     End of DHSEQR */

} /* dhseqr_ */

#undef z___ref
#undef s_ref
#undef h___ref