dnlaso.c

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

/*<          L = NBLOCK + 1 >*/
        l = *nblock + 1;
/*<          DO 120 K = I,NBLOCK >*/
        i__2 = *nblock;
        for (k = i__; k <= i__2; ++k) {
/*<             BET(I,K) = -BET(I,K) >*/
            bet[i__ + k * bet_dim1] = -bet[i__ + k * bet_dim1];
/*<             T(L,M) = -T(L,M) >*/
            t[l + m * t_dim1] = -t[l + m * t_dim1];
/*<             L = L - 1 >*/
            --l;
/*<             M = M + 1 >*/
            ++m;
/*<   120    CONTINUE >*/
/* L120: */
        }
/*<   130 CONTINUE >*/
L130:
        ;
    }

/* THIS IS THE LANCZOS STEP. */

/*<   160 CALL OP(N, NBLOCK, P1, P2) >*/
L160:
    (*op)(n, nblock, &p1[p1_offset], &p2[p2_offset]);
/*<       NOP = NOP + 1 >*/
    ++(*nop);
/*<       CALL IOVECT(N, NBLOCK, P1, J, 0) >*/
    (*iovect)(n, nblock, &p1[p1_offset], &j, &c__0);

/* THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) */

/*<       DO 180 I=1,NBLOCK >*/
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          DO 170 K=I,NBLOCK >*/
        i__2 = *nblock;
        for (k = i__; k <= i__2; ++k) {
/*<             CALL DAXPY(N, -BET(I,K), P0(1,K), 1, P2(1,I), 1) >*/
            d__1 = -bet[i__ + k * bet_dim1];
            daxpy_(n, &d__1, &p0[k * p0_dim1 + 1], &c__1, &p2[i__ * p2_dim1 + 
                    1], &c__1);
/*<   170    CONTINUE >*/
/* L170: */
        }
/*<   180 CONTINUE >*/
/* L180: */
    }

/* THIS COMPUTES ALP AND P2=P2-P1*ALP. */

/*<       DO 200 I=1,NBLOCK >*/
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          DO 190 K=1,I >*/
        i__2 = i__;
        for (k = 1; k <= i__2; ++k) {
/*<             II = I - K + 1 >*/
            ii = i__ - k + 1;
/*<             ALP(II,K) = DDOT(N,P1(1,I),1,P2(1,K),1) >*/
            alp[ii + k * alp_dim1] = ddot_(n, &p1[i__ * p1_dim1 + 1], &c__1, &
                    p2[k * p2_dim1 + 1], &c__1);
/*<             CALL DAXPY(N, -ALP(II,K), P1(1,I), 1, P2(1,K), 1) >*/
            d__1 = -alp[ii + k * alp_dim1];
            daxpy_(n, &d__1, &p1[i__ * p1_dim1 + 1], &c__1, &p2[k * p2_dim1 + 
                    1], &c__1);
/*<    >*/
            if (k != i__) {
                d__1 = -alp[ii + k * alp_dim1];
                daxpy_(n, &d__1, &p1[k * p1_dim1 + 1], &c__1, &p2[i__ * 
                        p2_dim1 + 1], &c__1);
            }
/*<   190   CONTINUE >*/
/* L190: */
        }
/*<   200 CONTINUE >*/
/* L200: */
    }

/*  REORTHOGONALIZATION OF THE SECOND BLOCK */

/*<       IF(J .NE. NBLOCK) GO TO 220 >*/
    if (j != *nblock) {
        goto L220;
    }
/*<       DO 215 I=1,NBLOCK >*/
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          DO 210 K=1,I >*/
        i__2 = i__;
        for (k = 1; k <= i__2; ++k) {
/*<             TEMP = DDOT(N,P1(1,I),1,P2(1,K),1) >*/
            temp = ddot_(n, &p1[i__ * p1_dim1 + 1], &c__1, &p2[k * p2_dim1 + 
                    1], &c__1);
/*<             CALL DAXPY(N, -TEMP, P1(1,I), 1, P2(1,K), 1) >*/
            d__1 = -temp;
            daxpy_(n, &d__1, &p1[i__ * p1_dim1 + 1], &c__1, &p2[k * p2_dim1 + 
                    1], &c__1);
/*<    >*/
            if (k != i__) {
                d__1 = -temp;
                daxpy_(n, &d__1, &p1[k * p1_dim1 + 1], &c__1, &p2[i__ * 
                        p2_dim1 + 1], &c__1);
            }
/*<             II = I - K + 1 >*/
            ii = i__ - k + 1;
/*<             ALP(II,K) = ALP(II,K) + TEMP             >*/
            alp[ii + k * alp_dim1] += temp;
/*<   210   CONTINUE >*/
/* L210: */
        }
/*<   215 CONTINUE >*/
/* L215: */
    }

/* THIS ORTHONORMALIZES THE NEXT BLOCK */

/*<   220 CALL DORTQR(N, N, NBLOCK, P2, BET) >*/
L220:
    dortqr_(n, n, nblock, &p2[p2_offset], &bet[bet_offset]);

/* THIS STORES ALP AND BET IN T. */

/*<       DO 250 I=1,NBLOCK >*/
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          M = J - NBLOCK + I >*/
        m = j - *nblock + i__;
/*<          DO 230 K=I,NBLOCK >*/
        i__2 = *nblock;
        for (k = i__; k <= i__2; ++k) {
/*<             L = K - I + 1 >*/
            l = k - i__ + 1;
/*<             T(L,M) = ALP(L,I) >*/
            t[l + m * t_dim1] = alp[l + i__ * alp_dim1];
/*<   230    CONTINUE >*/
/* L230: */
        }
/*<          DO 240 K=1,I >*/
        i__2 = i__;
        for (k = 1; k <= i__2; ++k) {
/*<             L = NBLOCK - I + K + 1 >*/
            l = *nblock - i__ + k + 1;
/*<             T(L,M) = BET(K,I) >*/
            t[l + m * t_dim1] = bet[k + i__ * bet_dim1];
/*<   240    CONTINUE >*/
/* L240: */
        }
/*<   250 CONTINUE >*/
/* L250: */
    }

/* THIS NEGATES T IF SMALL IS FALSE. */

/*<       IF (SMALL) GO TO 280 >*/
    if (*small) {
        goto L280;
    }
/*<       M = J - NBLOCK + 1 >*/
    m = j - *nblock + 1;
/*<       DO 270 I=M,J >*/
    i__1 = j;
    for (i__ = m; i__ <= i__1; ++i__) {
/*<          DO 260 K=1,L >*/
        i__2 = l;
        for (k = 1; k <= i__2; ++k) {
/*<             T(K,I) = -T(K,I) >*/
            t[k + i__ * t_dim1] = -t[k + i__ * t_dim1];
/*<   260    CONTINUE >*/
/* L260: */
        }
/*<   270 CONTINUE >*/
/* L270: */
    }

/* THIS SHIFTS THE LANCZOS VECTORS */

/*<   280 CALL DCOPY(NBLOCK*N, P1, 1, P0, 1) >*/
L280:
    i__1 = *nblock * *n;
    dcopy_(&i__1, &p1[p1_offset], &c__1, &p0[p0_offset], &c__1);
/*<       CALL DCOPY(NBLOCK*N, P2, 1, P1, 1) >*/
    i__1 = *nblock * *n;
    dcopy_(&i__1, &p2[p2_offset], &c__1, &p1[p1_offset], &c__1);
/*<       CALL DLAGER(J, NBAND, J-NBLOCK+1, T, TMIN, TMAX) >*/
    i__1 = j - *nblock + 1;
    dlager_(&j, nband, &i__1, &t[t_offset], &tmin, &tmax);
/*<       ANORM = DMAX1(RNORM, TMAX, -TMIN) >*/
/* Computing MAX */
    d__1 = max(rnorm,tmax), d__2 = -tmin;
    anorm = max(d__1,d__2);
/*<       IF (NUMBER.EQ.0) GO TO 305 >*/
    if (number == 0) {
        goto L305;
    }

/* THIS COMPUTES THE EXTREME EIGENVALUES OF ALP. */

/*<       CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) >*/
    dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
/*<    >*/
    dlaeig_(nblock, nblock, &c__1, &c__1, &alp[alp_offset], tarr, nblock, &p2[
            p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin, 
            &tmax);
/*<       ALPMIN = TARR(1) >*/
    alpmin = tarr[0];
/*<       CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) >*/
    dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
/*<    >*/
    dlaeig_(nblock, nblock, nblock, nblock, &alp[alp_offset], tarr, nblock, &
            p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &
            tmin, &tmax);
/*<       ALPMAX = TARR(1) >*/
    alpmax = tarr[0];

/* THIS COMPUTES ALP = BET(TRANSPOSE)*BET. */

/*<   305 DO 310 I = 1, NBLOCK >*/
L305:
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          DO 300 K = 1, I >*/
        i__2 = i__;
        for (k = 1; k <= i__2; ++k) {
/*<             L = I - K + 1 >*/
            l = i__ - k + 1;
/*<    >*/
            i__3 = *nblock - i__ + 1;
            alp[l + k * alp_dim1] = ddot_(&i__3, &bet[i__ + i__ * bet_dim1], 
                    nblock, &bet[k + i__ * bet_dim1], nblock);
/*<   300    CONTINUE >*/
/* L300: */
        }
/*<   310 CONTINUE >*/
/* L310: */
    }
/*<       IF(NUMBER .EQ. 0) GO TO 330 >*/
    if (number == 0) {
        goto L330;
    }

/* THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET. */

/*<       CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) >*/
    dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
/*<    >*/
    d__1 = anorm * anorm;
    dlaeig_(nblock, nblock, &c__1, &c__1, &alp[alp_offset], tarr, nblock, &p2[
            p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &c_b88,
             &d__1);
/*<       BETMIN = DSQRT(TARR(1)) >*/
    betmin = sqrt(tarr[0]);

/* THIS UPDATES TAU AND OTAU. */

/*<       DO 320 I=1,NUMBER >*/
    i__1 = number;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<    >*/
/* Computing MAX */
        d__1 = alpmax - val[i__], d__2 = val[i__] - alpmin;
        temp = (tau[i__] * max(d__1,d__2) + otau[i__] * betmax + *eps * anorm)
                 / betmin;
/*<          IF (I.LE.NPERM) TEMP = TEMP + RES(I)/BETMIN >*/
        if (i__ <= *nperm) {
            temp += res[i__] / betmin;
        }
/*<          OTAU(I) = TAU(I) >*/
        otau[i__] = tau[i__];
/*<          TAU(I) = TEMP >*/
        tau[i__] = temp;
/*<   320 CONTINUE >*/
/* L320: */
    }

/* THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET. */

/*<   330 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1) >*/
L330:
    dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
/*<    >*/
    d__1 = anorm * anorm;
    dlaeig_(nblock, nblock, nblock, nblock, &alp[alp_offset], tarr, nblock, &
            p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &
            c_b88, &d__1);
/*<       BETMAX = DSQRT(TARR(1)) >*/
    betmax = sqrt(tarr[0]);
/*<       IF (J.LE.2*NBLOCK) GO TO 80 >*/
    if (j <= *nblock << 1) {
        goto L80;
    }

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

/* THIS SECTION COMPUTES AND EXAMINES THE SMALLEST NONGOOD AND */
/* LARGEST DESIRED EIGENVALUES OF T TO SEE IF A CLOSER LOOK */
/* IS JUSTIFIED. */

/*<       TOLG = EPSRT*ANORM >*/
    tolg = epsrt * anorm;
/*<       TOLA = UTOL*RNORM >*/
    tola = utol * rnorm;
/*<    >*/
    if (*maxj - j < *nblock || (*nop >= *maxop && nleft != 0)) {
        goto L390;
    }
/*<       GO TO 400 >*/
    goto L400;

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

/* THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO */
/* SEE IF FURTHER ACTION IS INDICATED, ENTRY IS AT 380 OR 390 IF AN */
/* ITERATION (OR TERMINATION) IS KNOWN TO BE NEEDED, OTHERWISE ENTRY */
/* IS AT 400. */

/*<   380 J = J - NBLOCK >*/
L380:
    j -= *nblock;
/*<       IERR = -8 >*/
    *ierr = -8;
/*<   390 IF (NLEFT.EQ.0) RETURN >*/
L390:
    if (nleft == 0) {
        return 0;
    }
/*<       TEST = .TRUE. >*/
    test = TRUE_;
/*<   400 NTHETA = MIN0(J/2,NLEFT+1) >*/
L400:
/* Computing MIN */
    i__1 = j / 2, i__2 = nleft + 1;
    ntheta = min(i__1,i__2);
/*<    >*/
    dlaeig_(&j, nband, &c__1, &ntheta, &t[t_offset], &val[number + 1], maxj, &
            s[s_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin,
             &tmax);
/*<       CALL DMVPC(NBLOCK, BET, MAXJ, J, S, NTHETA, ATEMP, VTEMP, D) >*/
    dmvpc_(nblock, &bet[bet_offset], maxj, &j, &s[s_offset], &ntheta, &atemp[
            1], &vtemp[1], &d__[1]);

/* THIS CHECKS FOR TERMINATION OF A CHECK RUN */

/*<       IF(NLEFT .NE. 0 .OR. J .LT. 6*NBLOCK) GO TO 410 >*/
    if (nleft != 0 || j < *nblock * 6) {
        goto L410;
    }
/*<       IF(VAL(NUMBER+1)-ATEMP(1) .GT. VAL(NPERM) - TOLA) GO TO 790  >*/
    if (val[number + 1] - atemp[1] > val[*nperm] - tola) {
        goto L790;
    }

/* THIS UPDATES NLEFT BY EXAMINING THE COMPUTED EIGENVALUES OF T */
/* TO DETERMINE IF SOME PERMANENT VALUES ARE NO LONGER DESIRED. */

/*<  410  IF (NTHETA.LE.NLEFT) GO TO 470 >*/
L410:
    if (ntheta <= nleft) {
        goto L470;
    }
/*<       IF (NPERM.EQ.0) GO TO 430 >*/
    if (*nperm == 0) {
        goto L430;
    }
/*<       M = NUMBER + NLEFT + 1 >*/
    m = number + nleft + 1;
/*<       IF (VAL(M).GE.VAL(NPERM)) GO TO 430 >*/
    if (val[m] >= val[*nperm]) {
        goto L430;
    }
/*<       NPERM = NPERM - 1 >*/
    --(*nperm);
/*<       NGOOD = 0 >*/
    ngood = 0;
/*<       NUMBER = NPERM >*/
    number = *nperm;
/*<       NLEFT = NLEFT + 1 >*/
    ++nleft;
/*<       GO TO 400 >*/
    goto L400;

/* THIS UPDATES DELTA. */

/*<   430 M = NUMBER + NLEFT + 1 >*/
L430:
    m = number + nleft + 1;
/*<       DELTA = DMIN1(DELTA,VAL(M)) >*/
/* Computing MIN */
    d__1 = *delta, d__2 = val[m];
    *delta = min(d__1,d__2);
/*<       ENOUGH = .TRUE. >*/
    enough = TRUE_;
/*<       IF(NLEFT .EQ. 0) GO TO 80 >*/
    if (nleft == 0) {
        goto L80;
    }
/*<       NTHETA = NLEFT >*/
    ntheta = nleft;
/*<       VTEMP(NTHETA+1) = 1 >*/
    vtemp[ntheta + 1] = 1.;

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

/* THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL. */

/* THIS CHECKS FOR ENOUGH ACCEPTABLE VALUES. */

/*<       IF (.NOT.(TEST .OR. ENOUGH)) GO TO 470 >*/
    if (! (test || enough)) {
        goto L470;
    }
/*<       DELTA = DMIN1(DELTA,ANORM) >*/
    *delta = min(*delta,anorm);
/*<       PNORM = DMAX1(RNORM,DMAX1(-VAL(NUMBER+1),DELTA)) >*/
/* Computing MAX */
/* Computing MAX */
    d__3 = -val[number + 1];
    d__1 = rnorm, d__2 = max(d__3,*delta);
    pnorm = max(d__1,d__2);
/*<       TOLA = UTOL*PNORM >*/
    tola = utol * pnorm;
/*<       NSTART = 0 >*/
    nstart = 0;
/*<       DO 460 I=1,NTHETA >*/
    i__1 = ntheta;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          M = NUMBER + I >*/
        m = number + i__;
/*<    >*/
/* Computing MIN */
        d__1 = atemp[i__] * atemp[i__] / (*delta - val[m]), d__2 = atemp[i__];
        if (min(d__1,d__2) > tola) {