dnlaso.c

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

            goto L450;
        }
/*<          IND(I) = -1 >*/
        ind[i__] = -1;
/*<          GO TO 460 >*/
        goto L460;

/*<   450    ENOUGH = .FALSE. >*/
L450:
        enough = FALSE_;
/*<          IF (.NOT.TEST) GO TO 470 >*/
        if (! test) {
            goto L470;
        }
/*<          IND(I) = 1 >*/
        ind[i__] = 1;
/*<          NSTART = NSTART + 1 >*/
        ++nstart;
/*<   460 CONTINUE >*/
L460:
        ;
    }

/*  COPY VALUES OF IND INTO VTEMP */

/*<       DO 465 I = 1,NTHETA >*/
    i__1 = ntheta;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          VTEMP(I) = DBLE(FLOAT(IND(I))) >*/
        vtemp[i__] = (doublereal) ((real) ind[i__]);
/*<   465 CONTINUE >*/
/* L465: */
    }
/*<       GO TO 500 >*/
    goto L500;

/* THIS CHECKS FOR NEW GOOD VECTORS. */

/*<   470 NG = 0 >*/
L470:
    ng = 0;
/*<       DO 490 I=1,NTHETA >*/
    i__1 = ntheta;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          IF (VTEMP(I).GT.TOLG) GO TO 480 >*/
        if (vtemp[i__] > tolg) {
            goto L480;
        }
/*<          NG = NG + 1 >*/
        ++ng;
/*<          VTEMP(I) = -1 >*/
        vtemp[i__] = -1.;
/*<          GO TO 490 >*/
        goto L490;

/*<   480    VTEMP(I) = 1 >*/
L480:
        vtemp[i__] = 1.;
/*<   490 CONTINUE >*/
L490:
        ;
    }

/*<       IF (NG.LE.NGOOD) GO TO 80 >*/
    if (ng <= ngood) {
        goto L80;
    }
/*<       NSTART = NTHETA - NG >*/
    nstart = ntheta - ng;

/* ------------------------------------------------------------------ */

/* THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS. */
/* IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED. */

/*<   500 TEST = TEST .AND. .NOT.ENOUGH >*/
L500:
    test = test && ! enough;
/*<       NGOOD = NTHETA - NSTART >*/
    ngood = ntheta - nstart;
/*<       NSTART = NSTART + 1 >*/
    ++nstart;
/*<       NTHETA = NTHETA + 1 >*/
    ++ntheta;

/* THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND */
/* EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR */
/* FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE */
/* EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS */
/* PAUSE. */

/*<       CALL DCOPY(NTHETA, VAL(NUMBER+1), 1, VAL(NPERM+1), 1) >*/
    dcopy_(&ntheta, &val[number + 1], &c__1, &val[*nperm + 1], &c__1);
/*<       IF (NSTART.EQ.0) GO TO 580 >*/
    if (nstart == 0) {
        goto L580;
    }
/*<       IF (NSTART.EQ.NTHETA) GO TO 530 >*/
    if (nstart == ntheta) {
        goto L530;
    }
/*<    >*/
    dvsort_(&ntheta, &vtemp[1], &atemp[1], &c__1, &val[*nperm + 1], maxj, &j, 
            &s[s_offset]);

/* THES ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING */
/* VECTORS. */

/*<   530 IF (.NOT.TEST) NSTART = 0 >*/
L530:
    if (! test) {
        nstart = 0;
    }
/*<       IF (.NOT.TEST) GO TO 580 >*/
    if (! test) {
        goto L580;
    }

/*  FIND MINIMUM ATEMP VALUE TO AVOID POSSIBLE OVERFLOW */

/*<       TEMP = ATEMP(1) >*/
    temp = atemp[1];
/*<       DO 535 I = 1, NSTART >*/
    i__1 = nstart;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          TEMP = DMIN1(TEMP, ATEMP(I)) >*/
/* Computing MIN */
        d__1 = temp, d__2 = atemp[i__];
        temp = min(d__1,d__2);
/*<   535 CONTINUE >*/
/* L535: */
    }
/*<       M = NGOOD + 1 >*/
    m = ngood + 1;
/*<       L = NGOOD + MIN0(NSTART,NBLOCK) >*/
    l = ngood + min(nstart,*nblock);
/*<       DO 540 I=M,L >*/
    i__1 = l;
    for (i__ = m; i__ <= i__1; ++i__) {
/*<          CALL DSCAL(J, TEMP/ATEMP(I), S(1,I), 1) >*/
        d__1 = temp / atemp[i__];
        dscal_(&j, &d__1, &s[i__ * s_dim1 + 1], &c__1);
/*<   540 CONTINUE >*/
/* L540: */
    }
/*<       M = (NSTART-1)/NBLOCK >*/
    m = (nstart - 1) / *nblock;
/*<       IF (M.EQ.0) GO TO 570 >*/
    if (m == 0) {
        goto L570;
    }
/*<       L = NGOOD + NBLOCK >*/
    l = ngood + *nblock;
/*<       DO 560 I=1,M >*/
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          DO 550 K=1,NBLOCK >*/
        i__2 = *nblock;
        for (k = 1; k <= i__2; ++k) {
/*<             L = L + 1 >*/
            ++l;
/*<             IF (L.GT.NTHETA) GO TO 570 >*/
            if (l > ntheta) {
                goto L570;
            }
/*<             I1 = NGOOD + K >*/
            i1 = ngood + k;
/*<             CALL DAXPY(J, TEMP/ATEMP(L), S(1,L), 1, S(1,I1), 1) >*/
            d__1 = temp / atemp[l];
            daxpy_(&j, &d__1, &s[l * s_dim1 + 1], &c__1, &s[i1 * s_dim1 + 1], 
                    &c__1);
/*<   550    CONTINUE >*/
/* L550: */
        }
/*<   560 CONTINUE >*/
/* L560: */
    }
/*<   570 NSTART = MIN0(NSTART,NBLOCK) >*/
L570:
    nstart = min(nstart,*nblock);

/* THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS. */

/*<   580 IF (NGOOD.EQ.0 .OR. .NOT.(TEST .OR. ENOUGH)) GO TO 600 >*/
L580:
    if (ngood == 0 || ! (test || enough)) {
        goto L600;
    }
/*<       DO 590 I=1,NGOOD >*/
    i__1 = ngood;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          M = NPERM + I >*/
        m = *nperm + i__;
/*<          RES(M) = ATEMP(I) >*/
        res[m] = atemp[i__];
/*<   590 CONTINUE >*/
/* L590: */
    }

/* THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE */
/* LANCZOS VECTORS. */

/*<   600 NUMBER = NPERM + NGOOD >*/
L600:
    number = *nperm + ngood;
/*<       IF (TEST .OR. ENOUGH) CALL DCOPY(N*NBLOCK, DZERO, 0, P1, 1) >*/
    if (test || enough) {
        i__1 = *n * *nblock;
        dcopy_(&i__1, dzero, &c__0, &p1[p1_offset], &c__1);
    }
/*<       IF (NGOOD.EQ.0) GO TO 620 >*/
    if (ngood == 0) {
        goto L620;
    }
/*<       M = NPERM + 1 >*/
    m = *nperm + 1;
/*<       DO 610 I=M,NUMBER >*/
    i__1 = number;
    for (i__ = m; i__ <= i__1; ++i__) {
/*<          CALL DCOPY(N, DZERO, 0, VEC(1,I), 1) >*/
        dcopy_(n, dzero, &c__0, &vec[i__ * vec_dim1 + 1], &c__1);
/*<   610 CONTINUE >*/
/* L610: */
    }
/*<   620 DO 670 I=NBLOCK,J,NBLOCK >*/
L620:
    i__1 = j;
    i__2 = *nblock;
    for (i__ = *nblock; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
/*<          CALL IOVECT(N, NBLOCK, P2, I, 1) >*/
        (*iovect)(n, nblock, &p2[p2_offset], &i__, &c__1);
/*<          DO 660 K=1,NBLOCK >*/
        i__3 = *nblock;
        for (k = 1; k <= i__3; ++k) {
/*<             M = I - NBLOCK + K >*/
            m = i__ - *nblock + k;
/*<             IF (NSTART.EQ.0) GO TO 640 >*/
            if (nstart == 0) {
                goto L640;
            }
/*<             DO 630 L=1,NSTART >*/
            i__4 = nstart;
            for (l = 1; l <= i__4; ++l) {
/*<                I1 = NGOOD + L >*/
                i1 = ngood + l;
/*<                CALL DAXPY(N, S(M,I1), P2(1,K), 1, P1(1,L), 1) >*/
                daxpy_(n, &s[m + i1 * s_dim1], &p2[k * p2_dim1 + 1], &c__1, &
                        p1[l * p1_dim1 + 1], &c__1);
/*<   630       CONTINUE >*/
/* L630: */
            }
/*<   640       IF (NGOOD.EQ.0) GO TO 660 >*/
L640:
            if (ngood == 0) {
                goto L660;
            }
/*<             DO 650 L=1,NGOOD >*/
            i__4 = ngood;
            for (l = 1; l <= i__4; ++l) {
/*<                I1 = L + NPERM >*/
                i1 = l + *nperm;
/*<                CALL DAXPY(N, S(M,L), P2(1,K), 1, VEC(1,I1), 1) >*/
                daxpy_(n, &s[m + l * s_dim1], &p2[k * p2_dim1 + 1], &c__1, &
                        vec[i1 * vec_dim1 + 1], &c__1);
/*<   650       CONTINUE >*/
/* L650: */
            }
/*<   660    CONTINUE >*/
L660:
            ;
        }
/*<   670 CONTINUE >*/
/* L670: */
    }
/*<       IF (TEST .OR. ENOUGH) GO TO 690 >*/
    if (test || enough) {
        goto L690;
    }

/* THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE */
/* TAU RECURRENCE. */

/*<       M = NPERM + 1 >*/
    m = *nperm + 1;
/*<       DO 680 I=M,NUMBER >*/
    i__2 = number;
    for (i__ = m; i__ <= i__2; ++i__) {
/*<          TEMP = 1.0D0/DNRM2(N,VEC(1,I),1) >*/
        temp = 1. / dnrm2_(n, &vec[i__ * vec_dim1 + 1], &c__1);
/*<          CALL DSCAL(N, TEMP, VEC(1,I), 1) >*/
        dscal_(n, &temp, &vec[i__ * vec_dim1 + 1], &c__1);
/*<          TAU(I) = 1.0D0 >*/
        tau[i__] = 1.;
/*<          OTAU(I) = 1.0D0 >*/
        otau[i__] = 1.;
/*<   680 CONTINUE >*/
/* L680: */
    }

/*  SHIFT S VECTORS TO ALIGN FOR LATER CALL TO DLAEIG */

/*<       CALL DCOPY(NTHETA, VAL(NPERM+1), 1, VTEMP, 1) >*/
    dcopy_(&ntheta, &val[*nperm + 1], &c__1, &vtemp[1], &c__1);
/*<       CALL DVSORT(NTHETA, VTEMP, ATEMP, 0, TARR, MAXJ, J, S) >*/
    dvsort_(&ntheta, &vtemp[1], &atemp[1], &c__0, tarr, maxj, &j, &s[s_offset]
            );
/*<       GO TO 80 >*/
    goto L80;

/* ------------------------------------------------------------------ */

/* THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE */
/* PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING */
/* THE PERMANENT VECTORS. */

/*<   690 IF (NGOOD.EQ.0 .AND. NOP.GE.MAXOP) GO TO 810 >*/
L690:
    if (ngood == 0 && *nop >= *maxop) {
        goto L810;
    }
/*<       IF (NGOOD.EQ.0) GO TO 30 >*/
    if (ngood == 0) {
        goto L30;
    }

/* THIS ORTHONORMALIZES THE VECTORS */

/*<       CALL DORTQR(NMVEC, N, NPERM+NGOOD, VEC, S) >*/
    i__2 = *nperm + ngood;
    dortqr_(nmvec, n, &i__2, &vec[vec_offset], &s[s_offset]);

/* THIS SORTS THE VALUES AND VECTORS. */

/*<    >*/
    if (*nperm != 0) {
        i__2 = *nperm + ngood;
        dvsort_(&i__2, &val[1], &res[1], &c__0, &temp, nmvec, n, &vec[
                vec_offset]);
    }
/*<       NPERM = NPERM + NGOOD >*/
    *nperm += ngood;
/*<       NLEFT = NLEFT - NGOOD >*/
    nleft -= ngood;
/*<       RNORM = DMAX1(-VAL(1),VAL(NPERM)) >*/
/* Computing MAX */
    d__1 = -val[1], d__2 = val[*nperm];
    rnorm = max(d__1,d__2);

/* THIS DECIDES WHERE TO GO NEXT. */

/*<       IF (NOP.GE.MAXOP .AND. NLEFT.NE.0) GO TO 810 >*/
    if (*nop >= *maxop && nleft != 0) {
        goto L810;
    }
/*<       IF (NLEFT.NE.0) GO TO 30 >*/
    if (nleft != 0) {
        goto L30;
    }
/*<       IF (VAL(NVAL)-VAL(1).LT.TOLA) GO TO 790 >*/
    if (val[*nval] - val[1] < tola) {
        goto L790;
    }

/* THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED */
/* TO LOOK FOR UNDISCLOSED MULTIPLICITIES. */

/*<       M = NPERM - NBLOCK + 1 >*/
    m = *nperm - *nblock + 1;
/*<       IF (M.LE.0) RETURN >*/
    if (m <= 0) {
        return 0;
    }
/*<       DO 780 I=1,M >*/
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__) {
/*<          L = I + NBLOCK - 1 >*/
        l = i__ + *nblock - 1;
/*<          IF (VAL(L)-VAL(I).LT.TOLA) GO TO 30 >*/
        if (val[l] - val[i__] < tola) {
            goto L30;
        }
/*<   780 CONTINUE >*/
/* L780: */
    }

/* THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ */
/* PROCEDURE IS NEEDED. */

/*<   790 M = NPERM - NBLOCK >*/
L790:
    m = *nperm - *nblock;
/*<       IF (M.LE.0) RETURN >*/
    if (m <= 0) {
        return 0;
    }
/*<       DO 800 I=1,M >*/
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__) {
/*<          L = I + NBLOCK >*/
        l = i__ + *nblock;
/*<          IF (VAL(L)-VAL(I).GE.TOLA) GO TO 800 >*/
        if (val[l] - val[i__] >= tola) {
            goto L800;
        }
/*<          RARITZ = .TRUE. >*/
        *raritz = TRUE_;
/*<          RETURN >*/
        return 0;
/*<   800 CONTINUE >*/
L800:
        ;
    }

/*<       RETURN >*/
    return 0;

/* ------------------------------------------------------------------ */

/* THIS REPORTS THAT MAXOP WAS EXCEEDED. */

/*<   810 IERR = -2 >*/
L810:
    *ierr = -2;
/*<       GO TO 790 >*/
    goto L790;

/*<       END >*/
} /* dnwla_ */

#ifdef __cplusplus
        }
#endif