hqr2.c

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

/*<          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, n >*/
        i__2 = *n;
        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 = 1, j >*/
        i__2 = j;
        for (i__ = 1; 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: */
        }
/*     .......... accumulate transformations .......... */
/*<          do 220 i = low, igh >*/
        i__2 = *igh;
        for (i__ = *low; i__ <= i__2; ++i__) {
/*<             p = x * z(i,k) + y * z(i,k+1) >*/
            p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1];
/*<             z(i,k) = z(i,k) - p >*/
            z__[i__ + k * z_dim1] -= p;
/*<             z(i,k+1) = z(i,k+1) - p * q >*/
            z__[i__ + (k + 1) * z_dim1] -= p * q;
/*<   220    continue >*/
/* L220: */
        }
/*<          go to 255 >*/
        goto L255;
/*<   225    continue >*/
L225:
/*     .......... row modification .......... */
/*<          do 230 j = k, n >*/
        i__2 = *n;
        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 = 1, j >*/
        i__2 = j;
        for (i__ = 1; 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: */
        }
/*     .......... accumulate transformations .......... */
/*<          do 250 i = low, igh >*/
        i__2 = *igh;
        for (i__ = *low; i__ <= i__2; ++i__) {
/*<             p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) >*/
            p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1] + 
                    zz * z__[i__ + (k + 2) * z_dim1];
/*<             z(i,k) = z(i,k) - p >*/
            z__[i__ + k * z_dim1] -= p;
/*<             z(i,k+1) = z(i,k+1) - p * q >*/
            z__[i__ + (k + 1) * z_dim1] -= p * q;
/*<             z(i,k+2) = z(i,k+2) - p * r >*/
            z__[i__ + (k + 2) * z_dim1] -= p * r__;
/*<   250    continue >*/
/* L250: */
        }
/*<   255    continue >*/
L255:

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

/*<       go to 70 >*/
    goto L70;
/*     .......... one root found .......... */
/*<   270 h(en,en) = x + t >*/
L270:
    h__[en + en * h_dim1] = x + t;
/*<       wr(en) = h(en,en) >*/
    wr[en] = h__[en + en * h_dim1];
/*<       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)));
/*<       h(en,en) = x + t >*/
    h__[en + en * h_dim1] = x + t;
/*<       x = h(en,en) >*/
    x = h__[en + en * h_dim1];
/*<       h(na,na) = y + t >*/
    h__[na + na * h_dim1] = y + 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.;
/*<       x = h(en,na) >*/
    x = h__[en + na * h_dim1];
/*<       s = dabs(x) + dabs(zz) >*/
    s = abs(x) + abs(zz);
/*<       p = x / s >*/
    p = x / s;
/*<       q = zz / s >*/
    q = zz / s;
/*<       r = dsqrt(p*p+q*q) >*/
    r__ = sqrt(p * p + q * q);
/*<       p = p / r >*/
    p /= r__;
/*<       q = q / r >*/
    q /= r__;
/*     .......... row modification .......... */
/*<       do 290 j = na, n >*/
    i__1 = *n;
    for (j = na; j <= i__1; ++j) {
/*<          zz = h(na,j) >*/
        zz = h__[na + j * h_dim1];
/*<          h(na,j) = q * zz + p * h(en,j) >*/
        h__[na + j * h_dim1] = q * zz + p * h__[en + j * h_dim1];
/*<          h(en,j) = q * h(en,j) - p * zz >*/
        h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * zz;
/*<   290 continue >*/
/* L290: */
    }
/*     .......... column modification .......... */
/*<       do 300 i = 1, en >*/
    i__1 = en;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          zz = h(i,na) >*/
        zz = h__[i__ + na * h_dim1];
/*<          h(i,na) = q * zz + p * h(i,en) >*/
        h__[i__ + na * h_dim1] = q * zz + p * h__[i__ + en * h_dim1];
/*<          h(i,en) = q * h(i,en) - p * zz >*/
        h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * zz;
/*<   300 continue >*/
/* L300: */
    }
/*     .......... accumulate transformations .......... */
/*<       do 310 i = low, igh >*/
    i__1 = *igh;
    for (i__ = *low; i__ <= i__1; ++i__) {
/*<          zz = z(i,na) >*/
        zz = z__[i__ + na * z_dim1];
/*<          z(i,na) = q * zz + p * z(i,en) >*/
        z__[i__ + na * z_dim1] = q * zz + p * z__[i__ + en * z_dim1];
/*<          z(i,en) = q * z(i,en) - p * zz >*/
        z__[i__ + en * z_dim1] = q * z__[i__ + en * z_dim1] - p * zz;
/*<   310 continue >*/
/* L310: */
    }

/*<       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;
/*     .......... all roots found.  backsubstitute to find */
/*                vectors of upper triangular form .......... */
/*<   340 if (norm .eq. 0.0d0) go to 1001 >*/
L340:
    if (norm == 0.) {
        goto L1001;
    }
/*     .......... for en=n step -1 until 1 do -- .......... */
/*<       do 800 nn = 1, n >*/
    i__1 = *n;
    for (nn = 1; nn <= i__1; ++nn) {
/*<          en = n + 1 - nn >*/
        en = *n + 1 - nn;
/*<          p = wr(en) >*/
        p = wr[en];
/*<          q = wi(en) >*/
        q = wi[en];
/*<          na = en - 1 >*/
        na = en - 1;
/*<          if (q) 710, 600, 800 >*/
        if (q < 0.) {
            goto L710;
        } else if (q == 0) {
            goto L600;
        } else {
            goto L800;
        }
/*     .......... real vector .......... */
/*<   600    m = en >*/
L600:
        m = en;
/*<          h(en,en) = 1.0d0 >*/
        h__[en + en * h_dim1] = 1.;
/*<          if (na .eq. 0) go to 800 >*/
        if (na == 0) {
            goto L800;
        }
/*     .......... for i=en-1 step -1 until 1 do -- .......... */
/*<          do 700 ii = 1, na >*/
        i__2 = na;
        for (ii = 1; ii <= i__2; ++ii) {
/*<             i = en - ii >*/
            i__ = en - ii;
/*<             w = h(i,i) - p >*/
            w = h__[i__ + i__ * h_dim1] - p;
/*<             r = 0.0d0 >*/
            r__ = 0.;

/*<             do 610 j = m, en >*/
            i__3 = en;
            for (j = m; j <= i__3; ++j) {
/*<   610       r = r + h(i,j) * h(j,en) >*/
/* L610: */
                r__ += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
            }

/*<             if (wi(i) .ge. 0.0d0) go to 630 >*/
            if (wi[i__] >= 0.) {
                goto L630;
            }