hqr.c

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

        q /= s;
/*<          r = r / s >*/
        r__ /= s;
/*<          if (m .eq. l) go to 150 >*/
        if (m == l) {
            goto L150;
        }
/*<          tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1))) >*/
        tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) + 
                abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
/*<          tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r)) >*/
        tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q) 
                + abs(r__));
/*<          if (tst2 .eq. tst1) go to 150 >*/
        if (tst2 == tst1) {
            goto L150;
        }
/*<   140 continue >*/
/* L140: */
    }

/*<   150 mp2 = m + 2 >*/
L150:
    mp2 = m + 2;

/*<       do 160 i = mp2, en >*/
    i__1 = en;
    for (i__ = mp2; i__ <= i__1; ++i__) {
/*<          h(i,i-2) = 0.0d0 >*/
        h__[i__ + (i__ - 2) * h_dim1] = 0.;
/*<          if (i .eq. mp2) go to 160 >*/
        if (i__ == mp2) {
            goto L160;
        }
/*<          h(i,i-3) = 0.0d0 >*/
        h__[i__ + (i__ - 3) * h_dim1] = 0.;
/*<   160 continue >*/
L160:
        ;
    }
/*     .......... double qr step involving rows l to en and */
/*                columns m to en .......... */
/*<       do 260 k = m, na >*/
    i__1 = na;
    for (k = m; k <= i__1; ++k) {
/*<          notlas = k .ne. na >*/
        notlas = k != na;
/*<          if (k .eq. m) go to 170 >*/
        if (k == m) {
            goto L170;
        }
/*<          p = h(k,k-1) >*/
        p = h__[k + (k - 1) * h_dim1];
/*<          q = h(k+1,k-1) >*/
        q = h__[k + 1 + (k - 1) * h_dim1];
/*<          r = 0.0d0 >*/
        r__ = 0.;
/*<          if (notlas) r = h(k+2,k-1) >*/
        if (notlas) {
            r__ = h__[k + 2 + (k - 1) * h_dim1];
        }
/*<          x = dabs(p) + dabs(q) + dabs(r) >*/
        x = abs(p) + abs(q) + abs(r__);
/*<          if (x .eq. 0.0d0) go to 260 >*/
        if (x == 0.) {
            goto L260;
        }
/*<          p = p / x >*/
        p /= x;
/*<          q = q / x >*/
        q /= x;
/*<          r = r / x >*/
        r__ /= x;
/*<   170    s = dsign(dsqrt(p*p+q*q+r*r),p) >*/
L170:
        d__1 = sqrt(p * p + q * q + r__ * r__);
        s = d_sign(&d__1, &p);
/*<          if (k .eq. m) go to 180 >*/
        if (k == m) {
            goto L180;
        }
/*<          h(k,k-1) = -s * x >*/
        h__[k + (k - 1) * h_dim1] = -s * x;
/*<          go to 190 >*/
        goto L190;
/*<   180    if (l .ne. m) h(k,k-1) = -h(k,k-1) >*/
L180:
        if (l != m) {
            h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
        }
/*<   190    p = p + s >*/
L190:
        p += s;
/*<          x = p / s >*/
        x = p / s;
/*<          y = q / s >*/
        y = q / s;
/*<          zz = r / s >*/
        zz = r__ / s;
/*<          q = q / p >*/
        q /= p;
/*<          r = r / p >*/
        r__ /= p;
/*<          if (notlas) go to 225 >*/
        if (notlas) {
            goto L225;
        }
/*     .......... row modification .......... */
/*<          do 200 j = k, EN >*/
        i__2 = en;
        for (j = k; j <= i__2; ++j) {
/*<             p = h(k,j) + q * h(k+1,j) >*/
            p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
/*<             h(k,j) = h(k,j) - p * x >*/
            h__[k + j * h_dim1] -= p * x;
/*<             h(k+1,j) = h(k+1,j) - p * y >*/
            h__[k + 1 + j * h_dim1] -= p * y;
/*<   200    continue >*/
/* L200: */
        }

/*<          j = min0(en,k+3) >*/
/* Computing MIN */
        i__2 = en, i__3 = k + 3;
        j = min(i__2,i__3);
/*     .......... column modification .......... */
/*<          do 210 i = L, j >*/
        i__2 = j;
        for (i__ = l; i__ <= i__2; ++i__) {
/*<             p = x * h(i,k) + y * h(i,k+1) >*/
            p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
/*<             h(i,k) = h(i,k) - p >*/
            h__[i__ + k * h_dim1] -= p;
/*<             h(i,k+1) = h(i,k+1) - p * q >*/
            h__[i__ + (k + 1) * h_dim1] -= p * q;
/*<   210    continue >*/
/* L210: */
        }
/*<          go to 255 >*/
        goto L255;
/*<   225    continue >*/
L225:
/*     .......... row modification .......... */
/*<          do 230 j = k, EN >*/
        i__2 = en;
        for (j = k; j <= i__2; ++j) {
/*<             p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) >*/
            p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
                    k + 2 + j * h_dim1];
/*<             h(k,j) = h(k,j) - p * x >*/
            h__[k + j * h_dim1] -= p * x;
/*<             h(k+1,j) = h(k+1,j) - p * y >*/
            h__[k + 1 + j * h_dim1] -= p * y;
/*<             h(k+2,j) = h(k+2,j) - p * zz >*/
            h__[k + 2 + j * h_dim1] -= p * zz;
/*<   230    continue >*/
/* L230: */
        }

/*<          j = min0(en,k+3) >*/
/* Computing MIN */
        i__2 = en, i__3 = k + 3;
        j = min(i__2,i__3);
/*     .......... column modification .......... */
/*<          do 240 i = L, j >*/
        i__2 = j;
        for (i__ = l; i__ <= i__2; ++i__) {
/*<             p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) >*/
            p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] + 
                    zz * h__[i__ + (k + 2) * h_dim1];
/*<             h(i,k) = h(i,k) - p >*/
            h__[i__ + k * h_dim1] -= p;
/*<             h(i,k+1) = h(i,k+1) - p * q >*/
            h__[i__ + (k + 1) * h_dim1] -= p * q;
/*<             h(i,k+2) = h(i,k+2) - p * r >*/
            h__[i__ + (k + 2) * h_dim1] -= p * r__;
/*<   240    continue >*/
/* L240: */
        }
/*<   255    continue >*/
L255:

/*<   260 continue >*/
L260:
        ;
    }

/*<       go to 70 >*/
    goto L70;
/*     .......... one root found .......... */
/*<   270 wr(en) = x + t >*/
L270:
    wr[en] = x + t;
/*<       wi(en) = 0.0d0 >*/
    wi[en] = 0.;
/*<       en = na >*/
    en = na;
/*<       go to 60 >*/
    goto L60;
/*     .......... two roots found .......... */
/*<   280 p = (y - x) / 2.0d0 >*/
L280:
    p = (y - x) / 2.;
/*<       q = p * p + w >*/
    q = p * p + w;
/*<       zz = dsqrt(dabs(q)) >*/
    zz = sqrt((abs(q)));
/*<       x = x + t >*/
    x += t;
/*<       if (q .lt. 0.0d0) go to 320 >*/
    if (q < 0.) {
        goto L320;
    }
/*     .......... real pair .......... */
/*<       zz = p + dsign(zz,p) >*/
    zz = p + d_sign(&zz, &p);
/*<       wr(na) = x + zz >*/
    wr[na] = x + zz;
/*<       wr(en) = wr(na) >*/
    wr[en] = wr[na];
/*<       if (zz .ne. 0.0d0) wr(en) = x - w / zz >*/
    if (zz != 0.) {
        wr[en] = x - w / zz;
    }
/*<       wi(na) = 0.0d0 >*/
    wi[na] = 0.;
/*<       wi(en) = 0.0d0 >*/
    wi[en] = 0.;
/*<       go to 330 >*/
    goto L330;
/*     .......... complex pair .......... */
/*<   320 wr(na) = x + p >*/
L320:
    wr[na] = x + p;
/*<       wr(en) = x + p >*/
    wr[en] = x + p;
/*<       wi(na) = zz >*/
    wi[na] = zz;
/*<       wi(en) = -zz >*/
    wi[en] = -zz;
/*<   330 en = enm2 >*/
L330:
    en = enm2;
/*<       go to 60 >*/
    goto L60;
/*     .......... set error -- all eigenvalues have not */
/*                converged after 30*n iterations .......... */
/*<  1000 ierr = en >*/
L1000:
    *ierr = en;
/*<  1001 return >*/
L1001:
    return 0;
/*<       end >*/
} /* hqr_ */

#ifdef __cplusplus
        }
#endif