hqr2.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*<             zz = w >*/
            zz = w;
/*<             s = r >*/
            s = r__;
/*<             go to 700 >*/
            goto L700;
/*<   630       m = i >*/
L630:
            m = i__;
/*<             if (wi(i) .ne. 0.0d0) go to 640 >*/
            if (wi[i__] != 0.) {
                goto L640;
            }
/*<             t = w >*/
            t = w;
/*<             if (t .ne. 0.0d0) go to 635 >*/
            if (t != 0.) {
                goto L635;
            }
/*<                tst1 = norm >*/
            tst1 = norm;
/*<                t = tst1 >*/
            t = tst1;
/*<   632          t = 0.01d0 * t >*/
L632:
            t *= .01;
/*<                tst2 = norm + t >*/
            tst2 = norm + t;
/*<                if (tst2 .gt. tst1) go to 632 >*/
            if (tst2 > tst1) {
                goto L632;
            }
/*<   635       h(i,en) = -r / t >*/
L635:
            h__[i__ + en * h_dim1] = -r__ / t;
/*<             go to 680 >*/
            goto L680;
/*     .......... solve real equations .......... */
/*<   640       x = h(i,i+1) >*/
L640:
            x = h__[i__ + (i__ + 1) * h_dim1];
/*<             y = h(i+1,i) >*/
            y = h__[i__ + 1 + i__ * h_dim1];
/*<             q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) >*/
            q = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__];
/*<             t = (x * s - zz * r) / q >*/
            t = (x * s - zz * r__) / q;
/*<             h(i,en) = t >*/
            h__[i__ + en * h_dim1] = t;
/*<             if (dabs(x) .le. dabs(zz)) go to 650 >*/
            if (abs(x) <= abs(zz)) {
                goto L650;
            }
/*<             h(i+1,en) = (-r - w * t) / x >*/
            h__[i__ + 1 + en * h_dim1] = (-r__ - w * t) / x;
/*<             go to 680 >*/
            goto L680;
/*<   650       h(i+1,en) = (-s - y * t) / zz >*/
L650:
            h__[i__ + 1 + en * h_dim1] = (-s - y * t) / zz;

/*     .......... overflow control .......... */
/*<   680       t = dabs(h(i,en)) >*/
L680:
            t = (d__1 = h__[i__ + en * h_dim1], abs(d__1));
/*<             if (t .eq. 0.0d0) go to 700 >*/
            if (t == 0.) {
                goto L700;
            }
/*<             tst1 = t >*/
            tst1 = t;
/*<             tst2 = tst1 + 1.0d0/tst1 >*/
            tst2 = tst1 + 1. / tst1;
/*<             if (tst2 .gt. tst1) go to 700 >*/
            if (tst2 > tst1) {
                goto L700;
            }
/*<             do 690 j = i, en >*/
            i__3 = en;
            for (j = i__; j <= i__3; ++j) {
/*<                h(j,en) = h(j,en)/t >*/
                h__[j + en * h_dim1] /= t;
/*<   690       continue >*/
/* L690: */
            }

/*<   700    continue >*/
L700:
            ;
        }
/*     .......... end real vector .......... */
/*<          go to 800 >*/
        goto L800;
/*     .......... complex vector .......... */
/*<   710    m = na >*/
L710:
        m = na;
/*     .......... last vector component chosen imaginary so that */
/*                eigenvector matrix is triangular .......... */
/*<          if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720 >*/
        if ((d__1 = h__[en + na * h_dim1], abs(d__1)) <= (d__2 = h__[na + en *
                 h_dim1], abs(d__2))) {
            goto L720;
        }
/*<          h(na,na) = q / h(en,na) >*/
        h__[na + na * h_dim1] = q / h__[en + na * h_dim1];
/*<          h(na,en) = -(h(en,en) - p) / h(en,na) >*/
        h__[na + en * h_dim1] = -(h__[en + en * h_dim1] - p) / h__[en + na * 
                h_dim1];
/*<          go to 730 >*/
        goto L730;
/*<   720    call cdiv(0.0d0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) >*/
L720:
        d__1 = -h__[na + en * h_dim1];
        d__2 = h__[na + na * h_dim1] - p;
        cdiv_(&c_b49, &d__1, &d__2, &q, &h__[na + na * h_dim1], &h__[na + en *
                 h_dim1]);
/*<   730    h(en,na) = 0.0d0 >*/
L730:
        h__[en + na * h_dim1] = 0.;
/*<          h(en,en) = 1.0d0 >*/
        h__[en + en * h_dim1] = 1.;
/*<          enm2 = na - 1 >*/
        enm2 = na - 1;
/*<          if (enm2 .eq. 0) go to 800 >*/
        if (enm2 == 0) {
            goto L800;
        }
/*     .......... for i=en-2 step -1 until 1 do -- .......... */
/*<          do 795 ii = 1, enm2 >*/
        i__2 = enm2;
        for (ii = 1; ii <= i__2; ++ii) {
/*<             i = na - ii >*/
            i__ = na - ii;
/*<             w = h(i,i) - p >*/
            w = h__[i__ + i__ * h_dim1] - p;
/*<             ra = 0.0d0 >*/
            ra = 0.;
/*<             sa = 0.0d0 >*/
            sa = 0.;

/*<             do 760 j = m, en >*/
            i__3 = en;
            for (j = m; j <= i__3; ++j) {
/*<                ra = ra + h(i,j) * h(j,na) >*/
                ra += h__[i__ + j * h_dim1] * h__[j + na * h_dim1];
/*<                sa = sa + h(i,j) * h(j,en) >*/
                sa += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
/*<   760       continue >*/
/* L760: */
            }

/*<             if (wi(i) .ge. 0.0d0) go to 770 >*/
            if (wi[i__] >= 0.) {
                goto L770;
            }
/*<             zz = w >*/
            zz = w;
/*<             r = ra >*/
            r__ = ra;
/*<             s = sa >*/
            s = sa;
/*<             go to 795 >*/
            goto L795;
/*<   770       m = i >*/
L770:
            m = i__;
/*<             if (wi(i) .ne. 0.0d0) go to 780 >*/
            if (wi[i__] != 0.) {
                goto L780;
            }
/*<             call cdiv(-ra,-sa,w,q,h(i,na),h(i,en)) >*/
            d__1 = -ra;
            d__2 = -sa;
            cdiv_(&d__1, &d__2, &w, &q, &h__[i__ + na * h_dim1], &h__[i__ + 
                    en * h_dim1]);
/*<             go to 790 >*/
            goto L790;
/*     .......... solve complex equations .......... */
/*<   780       x = h(i,i+1) >*/
L780:
            x = h__[i__ + (i__ + 1) * h_dim1];
/*<             y = h(i+1,i) >*/
            y = h__[i__ + 1 + i__ * h_dim1];
/*<             vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q >*/
            vr = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__] - q * q;
/*<             vi = (wr(i) - p) * 2.0d0 * q >*/
            vi = (wr[i__] - p) * 2. * q;
/*<             if (vr .ne. 0.0d0 .or. vi .ne. 0.0d0) go to 784 >*/
            if (vr != 0. || vi != 0.) {
                goto L784;
            }
/*<    >*/
            tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz));
/*<                vr = tst1 >*/
            vr = tst1;
/*<   783          vr = 0.01d0 * vr >*/
L783:
            vr *= .01;
/*<                tst2 = tst1 + vr >*/
            tst2 = tst1 + vr;
/*<                if (tst2 .gt. tst1) go to 783 >*/
            if (tst2 > tst1) {
                goto L783;
            }
/*<    >*/
L784:
            d__1 = x * r__ - zz * ra + q * sa;
            d__2 = x * s - zz * sa - q * ra;
            cdiv_(&d__1, &d__2, &vr, &vi, &h__[i__ + na * h_dim1], &h__[i__ + 
                    en * h_dim1]);
/*<             if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785 >*/
            if (abs(x) <= abs(zz) + abs(q)) {
                goto L785;
            }
/*<             h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x >*/
            h__[i__ + 1 + na * h_dim1] = (-ra - w * h__[i__ + na * h_dim1] + 
                    q * h__[i__ + en * h_dim1]) / x;
/*<             h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x >*/
            h__[i__ + 1 + en * h_dim1] = (-sa - w * h__[i__ + en * h_dim1] - 
                    q * h__[i__ + na * h_dim1]) / x;
/*<             go to 790 >*/
            goto L790;
/*<    >*/
L785:
            d__1 = -r__ - y * h__[i__ + na * h_dim1];
            d__2 = -s - y * h__[i__ + en * h_dim1];
            cdiv_(&d__1, &d__2, &zz, &q, &h__[i__ + 1 + na * h_dim1], &h__[
                    i__ + 1 + en * h_dim1]);

/*     .......... overflow control .......... */
/*<   790       t = dmax1(dabs(h(i,na)), dabs(h(i,en))) >*/
L790:
/* Computing MAX */
            d__3 = (d__1 = h__[i__ + na * h_dim1], abs(d__1)), d__4 = (d__2 = 
                    h__[i__ + en * h_dim1], abs(d__2));
            t = max(d__3,d__4);
/*<             if (t .eq. 0.0d0) go to 795 >*/
            if (t == 0.) {
                goto L795;
            }
/*<             tst1 = t >*/
            tst1 = t;
/*<             tst2 = tst1 + 1.0d0/tst1 >*/
            tst2 = tst1 + 1. / tst1;
/*<             if (tst2 .gt. tst1) go to 795 >*/
            if (tst2 > tst1) {
                goto L795;
            }
/*<             do 792 j = i, en >*/
            i__3 = en;
            for (j = i__; j <= i__3; ++j) {
/*<                h(j,na) = h(j,na)/t >*/
                h__[j + na * h_dim1] /= t;
/*<                h(j,en) = h(j,en)/t >*/
                h__[j + en * h_dim1] /= t;
/*<   792       continue >*/
/* L792: */
            }

/*<   795    continue >*/
L795:
            ;
        }
/*     .......... end complex vector .......... */
/*<   800 continue >*/
L800:
        ;
    }
/*     .......... end back substitution. */
/*                vectors of isolated roots .......... */
/*<       do 840 i = 1, n >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          if (i .ge. low .and. i .le. igh) go to 840 >*/
        if (i__ >= *low && i__ <= *igh) {
            goto L840;
        }

/*<          do 820 j = i, n >*/
        i__2 = *n;
        for (j = i__; j <= i__2; ++j) {
/*<   820    z(i,j) = h(i,j) >*/
/* L820: */
            z__[i__ + j * z_dim1] = h__[i__ + j * h_dim1];
        }

/*<   840 continue >*/
L840:
        ;
    }
/*     .......... multiply by transformation matrix to give */
/*                vectors of original full matrix. */
/*                for j=n step -1 until low do -- .......... */
/*<       do 880 jj = low, n >*/
    i__1 = *n;
    for (jj = *low; jj <= i__1; ++jj) {
/*<          j = n + low - jj >*/
        j = *n + *low - jj;
/*<          m = min0(j,igh) >*/
        m = min(j,*igh);

/*<          do 880 i = low, igh >*/
        i__2 = *igh;
        for (i__ = *low; i__ <= i__2; ++i__) {
/*<             zz = 0.0d0 >*/
            zz = 0.;

/*<             do 860 k = low, m >*/
            i__3 = m;
            for (k = *low; k <= i__3; ++k) {
/*<   860       zz = zz + z(i,k) * h(k,j) >*/
/* L860: */
                zz += z__[i__ + k * z_dim1] * h__[k + j * h_dim1];
            }

/*<             z(i,j) = zz >*/
            z__[i__ + j * z_dim1] = zz;
/*<   880 continue >*/
/* L880: */
        }
    }

/*<       go to 1001 >*/
    goto L1001;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
/*<  1000 ierr = en >*/
L1000:
    *ierr = en;
/*<  1001 return >*/
L1001:
    return 0;
/*<       end >*/
} /* hqr2_ */

#ifdef __cplusplus
        }
#endif