dlabcm.c

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

L50:
        if (rq <= bound[((j + 1) << 1) + 2]) {
            goto L100;
        }
/*<             IF(RQ + ERRB .LT. BOUND(1,J+2)) GO TO 100 >*/
        if (rq + errb < bound[((j + 2) << 1) + 1]) {
            goto L100;
        }

/*  SAVE THE REJECTED VECTOR IF INDICATED */

/*<             IF(RQ - ERRB .LE. BOUND(2,J+1)) GO TO 200 >*/
        if (rq - errb <= bound[((j + 1) << 1) + 2]) {
            goto L200;
        }
/*<             DO 60 I = J, NVAL >*/
        i__2 = nval;
        for (i__ = j; i__ <= i__2; ++i__) {
/*<                IF(BOUND(2,I+1) .GT. RQ) GO TO 70 >*/
            if (bound[((i__ + 1) << 1) + 2] > rq) {
                goto L70;
            }
/*<    60       CONTINUE >*/
/* L60: */
        }
/*<             GO TO 80 >*/
        goto L80;

/*<    70       CALL DCOPY(N, EIGVEC(1,J), 1, EIGVEC(1,I), 1) >*/
L70:
        dcopy_(n, &eigvec[j * eigvec_dim1 + 1], &c__1, &eigvec[i__ * 
                eigvec_dim1 + 1], &c__1);

/*<    80       CALL DLARAN(N, EIGVEC(1,J)) >*/
L80:
        dlaran_(n, &eigvec[j * eigvec_dim1 + 1]);
/*<             GO TO 200 >*/
        goto L200;

/*  PERTURB RQ TOWARD THE MIDDLE */

/*<   100       IF(SIGMA .LT. RQ) SIGMA = DMAX1(SIGMA, RQ-ERRB) >*/
L100:
        if (sigma < rq) {
/* Computing MAX */
            d__1 = sigma, d__2 = rq - errb;
            sigma = max(d__1,d__2);
        }
/*<             IF(SIGMA .GT. RQ) SIGMA = DMIN1(SIGMA, RQ+ERRB) >*/
        if (sigma > rq) {
/* Computing MIN */
            d__1 = sigma, d__2 = rq + errb;
            sigma = min(d__1,d__2);
        }

/*  FACTOR AND SOLVE */

/*<   200       DO 210 I = J, NVAL >*/
L200:
        i__2 = nval;
        for (i__ = j; i__ <= i__2; ++i__) {
/*<                IF(SIGMA .LT. BOUND(1,I+1)) GO TO 220 >*/
            if (sigma < bound[((i__ + 1) << 1) + 1]) {
                goto L220;
            }
/*<   210       CONTINUE >*/
/* L210: */
        }
/*<             I = NVAL + 1 >*/
        i__ = nval + 1;
/*<   220       NUMVEC = I - J >*/
L220:
        numvec = i__ - j;
/*<             NUMVEC = MIN0(NUMVEC, NBAND + 2) >*/
/* Computing MIN */
        i__2 = numvec, i__3 = *nband + 2;
        numvec = min(i__2,i__3);
/*<             IF(RESID .LT. ARTOL) NUMVEC = MIN0(1,NUMVEC) >*/
        if (resid < *artol) {
            numvec = min(1,numvec);
        }
/*<             CALL DCOPY(N, EIGVEC(1,J), 1, VTEMP, 1) >*/
        dcopy_(n, &eigvec[j * eigvec_dim1 + 1], &c__1, &vtemp[1], &c__1);
/*<    >*/
        i__2 = (*nband << 1) - 1;
        dlabfc_(n, nband, &a[a_offset], &sigma, &numvec, lde, &eigvec[j * 
                eigvec_dim1 + 1], &numl, &i__2, &atemp[1], &d__[1], atol);

/*  PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW */

/*<             IF(NUMVEC .EQ. 1) GO TO 227 >*/
        if (numvec == 1) {
            goto L227;
        }
/*<             L = NUMVEC - 1  >*/
        l = numvec - 1;
/*<             DO 225 I = 1,L >*/
        i__2 = l;
        for (i__ = 1; i__ <= i__2; ++i__) {
/*<                M = J + I >*/
            m = j + i__;
/*<                CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,M), 1) >*/
            d__1 = 1. / vnorm;
            dscal_(n, &d__1, &eigvec[m * eigvec_dim1 + 1], &c__1);
/*<   225       CONTINUE >*/
/* L225: */
        }

/*  UPDATE INTERVALS */

/*<   227       NUML = NUML - NL + 1 >*/
L227:
        numl = numl - *nl + 1;
/*<             IF(NUML .GE. 0) BOUND(2,1) = DMIN1(BOUND(2,1), SIGMA) >*/
        if (numl >= 0) {
            bound[4] = min(bound[4],sigma);
        }
/*<             DO 230 I = J, NVAL >*/
        i__2 = nval;
        for (i__ = j; i__ <= i__2; ++i__) {
/*<                IF(SIGMA .LT. BOUND(1,I+1)) GO TO 20 >*/
            if (sigma < bound[((i__ + 1) << 1) + 1]) {
                goto L20;
            }
/*<                IF(NUML .LT. I) BOUND(1,I+1) = SIGMA >*/
            if (numl < i__) {
                bound[((i__ + 1) << 1) + 1] = sigma;
            }
/*<                IF(NUML .GE. I) BOUND(2,I+1) = SIGMA >*/
            if (numl >= i__) {
                bound[((i__ + 1) << 1) + 2] = sigma;
            }
/*<   230       CONTINUE >*/
/* L230: */
        }
/*<    >*/
        if (numl < nval + 1) {
/* Computing MAX */
            d__1 = sigma, d__2 = bound[((nval + 2) << 1) + 1];
            bound[((nval + 2) << 1) + 1] = max(d__1,d__2);
        }
/*<             GO TO 20 >*/
        goto L20;

/*  ACCEPT AN EIGENPAIR */

/*<   300    CALL DLARAN(N, EIGVEC(1,J)) >*/
L300:
        dlaran_(n, &eigvec[j * eigvec_dim1 + 1]);
/*<          FLAG = .TRUE. >*/
        flag__ = TRUE_;
/*<          GO TO 310 >*/
        goto L310;

/*<   305    FLAG = .FALSE. >*/
L305:
        flag__ = FALSE_;
/*<          RQ = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1)) >*/
        rq = (bound[((j + 1) << 1) + 1] + bound[((j + 1) << 1) + 2]) * .5;
/*<    >*/
        i__2 = (*nband << 1) - 1;
        dlabfc_(n, nband, &a[a_offset], &rq, &numvec, lde, &eigvec[j * 
                eigvec_dim1 + 1], &numl, &i__2, &atemp[1], &d__[1], atol);
/*<          VNORM = DNRM2(N, EIGVEC(1,J), 1) >*/
        vnorm = dnrm2_(n, &eigvec[j * eigvec_dim1 + 1], &c__1);
/*<          IF(VNORM .NE. 0.0) CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1) >*/
        if (vnorm != (float)0.) {
            d__1 = 1. / vnorm;
            dscal_(n, &d__1, &eigvec[j * eigvec_dim1 + 1], &c__1);
        }

/*  ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES */

/*<   310    EIGVAL(J) = RQ >*/
L310:
        eigval[j] = rq;
/*<          IF(J .EQ. 1) GO TO 330 >*/
        if (j == 1) {
            goto L330;
        }
/*<          M = J - 1 >*/
        m = j - 1;
/*<          DO 320 I = 1, M >*/
        i__2 = m;
        for (i__ = 1; i__ <= i__2; ++i__) {
/*<    >*/
            d__1 = -ddot_(n, &eigvec[i__ * eigvec_dim1 + 1], &c__1, &eigvec[j 
                    * eigvec_dim1 + 1], &c__1);
            daxpy_(n, &d__1, &eigvec[i__ * eigvec_dim1 + 1], &c__1, &eigvec[j 
                    * eigvec_dim1 + 1], &c__1);
/*<   320    CONTINUE >*/
/* L320: */
        }
/*<   330    VNORM = DNRM2(N, EIGVEC(1,J), 1) >*/
L330:
        vnorm = dnrm2_(n, &eigvec[j * eigvec_dim1 + 1], &c__1);
/*<          IF(VNORM .EQ. 0.0D0) GO TO 305 >*/
        if (vnorm == 0.) {
            goto L305;
        }
/*<          CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1) >*/
        d__1 = 1. / vnorm;
        dscal_(n, &d__1, &eigvec[j * eigvec_dim1 + 1], &c__1);

/*   ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE */

/*<          IF(FLAG) GO TO 305 >*/
        if (flag__) {
            goto L305;
        }
/*<          IF(J .EQ. NVAL) RETURN >*/
        if (j == nval) {
            return 0;
        }
/*<          M = J + 1 >*/
        m = j + 1;
/*<          DO 340 I = M, NVAL >*/
        i__2 = nval;
        for (i__ = m; i__ <= i__2; ++i__) {
/*<    >*/
            d__1 = -ddot_(n, &eigvec[j * eigvec_dim1 + 1], &c__1, &eigvec[i__ 
                    * eigvec_dim1 + 1], &c__1);
            daxpy_(n, &d__1, &eigvec[j * eigvec_dim1 + 1], &c__1, &eigvec[i__ 
                    * eigvec_dim1 + 1], &c__1);
/*<   340    CONTINUE >*/
/* L340: */
        }
/*<   400 CONTINUE >*/
/* L400: */
    }
/*<       RETURN >*/
    return 0;

/*<   500 CONTINUE >*/
/* L500: */
/*<       END >*/
    return 0;
} /* dlabcm_ */

#ifdef __cplusplus
        }
#endif