dnlaso.c

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

    }

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

/* THIS SECTION CALLS DNWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM */
/* WITH SELECTIVE ORTHOGONALIZATION. */

/*<   130 NBAND = NBLOCK + 1 >*/
L130:
    nband = *nblock + 1;
/*<       I1 = 1 + N*NBLOCK >*/
    i1 = *n * *nblock + 1;
/*<       I2 = I1 + N*NBLOCK >*/
    i2 = i1 + *n * *nblock;
/*<       I3 = I2 + NV >*/
    i3 = i2 + nv;
/*<       I4 = I3 + NV >*/
    i4 = i3 + nv;
/*<       I5 = I4 + NV >*/
    i5 = i4 + nv;
/*<       I6 = I5 + MAXJ*NBAND >*/
    i6 = i5 + *maxj * nband;
/*<       I7 = I6 + NBLOCK*NBLOCK >*/
    i7 = i6 + *nblock * *nblock;
/*<       I8 = I7 + NBLOCK*NBLOCK >*/
    i8 = i7 + *nblock * *nblock;
/*<       I9 = I8 + MAXJ*(NV+1) >*/
    i9 = i8 + *maxj * (nv + 1);
/*<       I10 = I9 + NBLOCK >*/
    i10 = i9 + *nblock;
/*<       I11 = I10 + 2*NV + 6 >*/
    i11 = i10 + (nv << 1) + 6;
/*<       I12 = I11 + MAXJ*(2*NBLOCK+1) >*/
    i12 = i11 + *maxj * ((*nblock << 1) + 1);
/*<       I13 = I12 + MAXJ >*/
    i13 = i12 + *maxj;
/*<    >*/
    dnwla_(op, iovect, n, &nband, &nv, nfig, nperm, &val[
            val_offset], nmvec, &vec[vec_offset], nblock, maxop, maxj, &nop, &
            work[1], &work[i1], &work[i2], &work[i3], &work[i4], &work[i5], &
            work[i6], &work[i7], &work[i8], &work[i9], &work[i10], &work[i11],
             &work[i12], &work[i13], &ind[1], &small, &raritz, &delta, &eps, 
            ierr);

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

/* THIS SECTION CALLS DNPPLA (THE POST PROCESSOR). */

/*<       IF (NPERM.EQ.0) GO TO 140 >*/
    if (*nperm == 0) {
        goto L140;
    }
/*<       I1 = N*NBLOCK + 1 >*/
    i1 = *n * *nblock + 1;
/*<       I2 = I1 + NPERM*NPERM >*/
    i2 = i1 + *nperm * *nperm;
/*<       I3 = I2 + NPERM*NPERM >*/
    i3 = i2 + *nperm * *nperm;
/*<       I4 = I3 + MAX0(N*NBLOCK,2*NPERM*NPERM) >*/
/* Computing MAX */
    i__1 = *n * *nblock, i__2 = (*nperm << 1) * *nperm;
    i4 = i3 + max(i__1,i__2);
/*<       I5 = I4 + N*NBLOCK >*/
    i5 = i4 + *n * *nblock;
/*<       I6 = I5 + 2*NPERM + 4 >*/
    i6 = i5 + (*nperm << 1) + 4;
/*<    >*/
    dnppla_(op, iovect, n, nperm, &nop, nmval, &val[val_offset], 
            nmvec, &vec[vec_offset], nblock, &work[i1], &work[i2], &work[i3], 
            &work[i4], &work[i5], &work[i6], &delta, &small, &raritz, &eps);

/*<   140 IND(1) = NOP >*/
L140:
    ind[1] = nop;
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* dnlaso_ */


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

/*<    >*/
/* Subroutine */ int dnwla_(
        void (*op)(integer*,integer*,doublereal*,doublereal*),
        void (*iovect)(integer*,integer*,doublereal*,integer*,integer*),
        integer *n, integer *nband, 
        integer *nval, integer *nfig, integer *nperm, doublereal *val, 
        integer *nmvec, doublereal *vec, integer *nblock, integer *maxop, 
        integer *maxj, integer *nop, doublereal *p1, doublereal *p0, 
        doublereal *res, doublereal *tau, doublereal *otau, doublereal *t, 
        doublereal *alp, doublereal *bet, doublereal *s, doublereal *p2, 
        doublereal *bound, doublereal *atemp, doublereal *vtemp, doublereal *
        d__, integer *ind, logical *small, logical *raritz, doublereal *delta,
         doublereal *eps, integer *ierr)
{
    /* System generated locals */
    integer vec_dim1, vec_offset, p0_dim1, p0_offset, p1_dim1, p1_offset, 
            p2_dim1, p2_offset, t_dim1, t_offset, alp_dim1, alp_offset, 
            bet_dim1, bet_offset, s_dim1, s_offset, i__1, i__2, i__3, i__4;
    doublereal d__1, d__2, d__3;

    /* Builtin functions */
    double sqrt(doublereal), pow_dd(doublereal *, doublereal *);

    /* Local variables */
    integer i__, j, k, l=0, m, i1, ii, ng;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
            integer *);
    doublereal tola=0, temp, tolg=0, tmin, tmax, tarr[1];
    logical test;
    doublereal utol;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
            integer *);
    integer ngood, nleft;
    doublereal anorm=0;
    extern /* Subroutine */ int dmvpc_(integer *, doublereal *, integer *, 
            integer *, doublereal *, integer *, doublereal *, doublereal *, 
            doublereal *), dcopy_(integer *, doublereal *, integer *, 
            doublereal *, integer *);
    integer mtemp;
    doublereal dzero[1];
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
            integer *, doublereal *, integer *);
    doublereal pnorm, epsrt, rnorm;
    extern /* Subroutine */ int dlaeig_(integer *, integer *, integer *, 
            integer *, doublereal *, doublereal *, integer *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *), dlager_(integer *, 
            integer *, integer *, doublereal *, doublereal *, doublereal *), 
            dlaran_(integer *, doublereal *);
    doublereal betmin, alpmin=0, betmax, alpmax=0;
    integer ntheta;
    logical enough;
    integer number, nstart;
    extern /* Subroutine */ int dortqr_(integer *, integer *, integer *, 
            doublereal *, doublereal *), dvsort_(integer *, doublereal *, 
            doublereal *, integer *, doublereal *, integer *, integer *, 
            doublereal *);


/*<    >*/
/*<       LOGICAL RARITZ, SMALL >*/
/*<    >*/
/*<       EXTERNAL OP, IOVECT >*/

/* DNWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE */
/* ORTHOGONALIZATION. */

/*   NBAND  NBLOCK + 1  THE BAND WIDTH OF T. */

/*   NVAL   THE NUMBER OF DESIRED EIGENVALUES. */

/*   NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS */
/*     INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE */
/*     ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL */
/*     TO THE CURRENT KRYLOV SUBSPACE. */

/*   NOP   THE NUMBER OF CALLS TO OP. */

/*   P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS. */

/*   RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS. */

/*   TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION. */

/*   T   THE BAND MATRIX. */

/*   ALP   THE CURRENT DIAGONAL BLOCK. */

/*   BET   THE CURRENT OFF DIAGONAL BLOCK. */

/*   BOUND, ATEMP, VTEMP, D  TEMPORARY STORAGE USED BY THE BAND */
/*      EIGENVALUE SOLVER DLAEIG. */

/*   S   EIGENVECTORS OF T. */

/*   SMALL   .TRUE.  IF THE SMALL EIGENVALUES ARE DESIRED. */

/*   RARITZ  RETURNED AS  .TRUE.  IF A FINAL RAYLEIGH-RITZ PROCEDURE */
/*     IS TO BE DONE. */

/*   DELTA   RETURNED AS THE VALUE OF THE (NVAL+1)TH EIGENVALUE */
/*     OF THE MATRIX.  USED IN ESTIMATING THE ACCURACY OF THE */
/*     COMPUTED EIGENVALUES. */


/* INTERNAL VARIABLES USED */

/*<    >*/
/*<       LOGICAL ENOUGH, TEST >*/
/*<    >*/
/*<    >*/

/* J   THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT */
/*    KRYLOV SUBSPACE. */

/* NGOOD   THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS */
/*    LIE IN THE CURRENT KRYLOV SUBSPACE). */

/* NLEFT   THE NUMBER OF VALUES WHICH REMAIN TO BE DETERMINED, */
/*    I.E.  NLEFT = NVAL - NPERM. */

/* NUMBER = NPERM + NGOOD. */

/* ANORM   AN ESTIMATE OF THE NORM OF THE MATRIX. */

/* EPS   THE RELATIVE MACHINE PRECISION. */

/* UTOL   THE USER TOLERANCE. */

/* TARR  AN ARRAY OF LENGTH ONE USED TO INSURE TYPE CONSISTENCY IN CALLS TO */
/*       DLAEIG */

/* DZERO AN ARRAY OF LENGTH ONE CONTAINING DZERO, USED TO INSURE TYPE */
/*       CONSISTENCY IN CALLS TO DCOPY */

/*<       DZERO(1) = 0.0D0 >*/
    /* Parameter adjustments */
    p2_dim1 = *n;
    p2_offset = 1 + p2_dim1;
    p2 -= p2_offset;
    p0_dim1 = *n;
    p0_offset = 1 + p0_dim1;
    p0 -= p0_offset;
    p1_dim1 = *n;
    p1_offset = 1 + p1_dim1;
    p1 -= p1_offset;
    t_dim1 = *nband;
    t_offset = 1 + t_dim1;
    t -= t_offset;
    --val;
    vec_dim1 = *nmvec;
    vec_offset = 1 + vec_dim1;
    vec -= vec_offset;
    bet_dim1 = *nblock;
    bet_offset = 1 + bet_dim1;
    bet -= bet_offset;
    alp_dim1 = *nblock;
    alp_offset = 1 + alp_dim1;
    alp -= alp_offset;
    s_dim1 = *maxj;
    s_offset = 1 + s_dim1;
    s -= s_offset;
    --res;
    --tau;
    --otau;
    --bound;
    --atemp;
    --vtemp;
    --d__;
    --ind;

    /* Function Body */
    dzero[0] = 0.;
/*<       RNORM = 0.0D0 >*/
    rnorm = 0.;
/*<       IF (NPERM.NE.0) RNORM = DMAX1(-VAL(1),VAL(NPERM)) >*/
    if (*nperm != 0) {
/* Computing MAX */
        d__1 = -val[1], d__2 = val[*nperm];
        rnorm = max(d__1,d__2);
    }
/*<       PNORM = RNORM >*/
    pnorm = rnorm;
/*<       DELTA = 10.D30 >*/
    *delta = 1e31;
/*<       EPSRT = DSQRT(EPS) >*/
    epsrt = sqrt(*eps);
/*<       NLEFT = NVAL - NPERM >*/
    nleft = *nval - *nperm;
/*<       NOP = 0 >*/
    *nop = 0;
/*<       NUMBER = NPERM >*/
    number = *nperm;
/*<       RARITZ = .FALSE. >*/
    *raritz = FALSE_;
/*<       UTOL = DMAX1(DBLE(FLOAT(N))*EPS,10.0D0**DBLE((-FLOAT(NFIG)))) >*/
/* Computing MAX */
    d__3 = (doublereal) (-((real) (*nfig)));
    d__1 = (doublereal) ((real) (*n)) * *eps, d__2 = pow_dd(&c_b15, &d__3);
    utol = max(d__1,d__2);
/*<       J = MAXJ >*/
    j = *maxj;

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

/* ANY ITERATION OF THE ALGORITHM BEGINS HERE. */

/*<    30 DO 50 I=1,NBLOCK >*/
L30:
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          TEMP = DNRM2(N,P1(1,I),1) >*/
        temp = dnrm2_(n, &p1[i__ * p1_dim1 + 1], &c__1);
/*<          IF (TEMP.EQ.0D0) CALL DLARAN(N, P1(1,I)) >*/
        if (temp == 0.) {
            dlaran_(n, &p1[i__ * p1_dim1 + 1]);
        }
/*<    50 CONTINUE >*/
/* L50: */
    }
/*<       IF (NPERM.EQ.0) GO TO 70 >*/
    if (*nperm == 0) {
        goto L70;
    }
/*<       DO 60 I=1,NPERM >*/
    i__1 = *nperm;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          TAU(I) = 1.0D0 >*/
        tau[i__] = 1.;
/*<          OTAU(I) = 0.0D0 >*/
        otau[i__] = 0.;
/*<    60 CONTINUE >*/
/* L60: */
    }
/*<    70 CALL DCOPY(N*NBLOCK, DZERO, 0, P0, 1) >*/
L70:
    i__1 = *n * *nblock;
    dcopy_(&i__1, dzero, &c__0, &p0[p0_offset], &c__1);
/*<       CALL DCOPY(NBLOCK*NBLOCK, DZERO, 0, BET, 1) >*/
    i__1 = *nblock * *nblock;
    dcopy_(&i__1, dzero, &c__0, &bet[bet_offset], &c__1);
/*<       CALL DCOPY(J*NBAND, DZERO, 0, T, 1) >*/
    i__1 = j * *nband;
    dcopy_(&i__1, dzero, &c__0, &t[t_offset], &c__1);
/*<       MTEMP = NVAL + 1 >*/
    mtemp = *nval + 1;
/*<       DO 75 I = 1, MTEMP >*/
    i__1 = mtemp;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          CALL DCOPY(J, DZERO, 0, S(1,I), 1) >*/
        dcopy_(&j, dzero, &c__0, &s[i__ * s_dim1 + 1], &c__1);
/*<    75 CONTINUE >*/
/* L75: */
    }
/*<       NGOOD = 0 >*/
    ngood = 0;
/*<       TMIN = 1.0D30 >*/
    tmin = 1e30;
/*<       TMAX = -1.0D30 >*/
    tmax = -1e30;
/*<       TEST = .TRUE. >*/
    test = TRUE_;
/*<       ENOUGH = .FALSE. >*/
    enough = FALSE_;
/*<       BETMAX = 0.0D0 >*/
    betmax = 0.;
/*<       J = 0 >*/
    j = 0;

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

/* THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. */

/*<    80 J = J + NBLOCK >*/
L80:
    j += *nblock;

/* THIS IS THE SELECTIVE ORTHOGONALIZATION. */

/*<       IF (NUMBER.EQ.0) GO TO 110 >*/
    if (number == 0) {
        goto L110;
    }
/*<       DO 100 I=1,NUMBER >*/
    i__1 = number;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          IF (TAU(I).LT.EPSRT) GO TO 100 >*/
        if (tau[i__] < epsrt) {
            goto L100;
        }
/*<          TEST = .TRUE. >*/
        test = TRUE_;
/*<          TAU(I) = 0.0D0 >*/
        tau[i__] = 0.;
/*<          IF (OTAU(I).NE.0.0D0) OTAU(I) = 1.0D0 >*/
        if (otau[i__] != 0.) {
            otau[i__] = 1.;
        }
/*<          DO 90 K=1,NBLOCK >*/
        i__2 = *nblock;
        for (k = 1; k <= i__2; ++k) {
/*<             TEMP = -DDOT(N,VEC(1,I),1,P1(1,K),1) >*/
            temp = -ddot_(n, &vec[i__ * vec_dim1 + 1], &c__1, &p1[k * p1_dim1 
                    + 1], &c__1);
/*<             CALL DAXPY(N, TEMP, VEC(1,I), 1, P1(1,K), 1) >*/
            daxpy_(n, &temp, &vec[i__ * vec_dim1 + 1], &c__1, &p1[k * p1_dim1 
                    + 1], &c__1);

/* THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A */
/* NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS */
/* TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. */

/*<    >*/
            if ((d__1 = temp * bet[k + k * bet_dim1], abs(d__1)) > (
                    doublereal) ((real) (*n)) * epsrt * anorm && i__ > *nperm)
                     {
                goto L380;
            }
/*<    90    CONTINUE >*/
/* L90: */
        }
/*<   100 CONTINUE >*/
L100:
        ;
    }

/* IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. */

/*<   110 IF(.NOT. TEST) GO TO 160 >*/
L110:
    if (! test) {
        goto L160;
    }
/*<       CALL DORTQR(N, N, NBLOCK, P1, ALP) >*/
    dortqr_(n, n, nblock, &p1[p1_offset], &alp[alp_offset]);
/*<       TEST = .FALSE. >*/
    test = FALSE_;
/*<       IF(J .EQ. NBLOCK) GO TO 160 >*/
    if (j == *nblock) {
        goto L160;
    }
/*<       DO 130 I = 1,NBLOCK >*/
    i__1 = *nblock;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          IF(ALP(I,I) .GT. 0.0D0) GO TO 130 >*/
        if (alp[i__ + i__ * alp_dim1] > 0.) {
            goto L130;
        }
/*<          M = J - 2*NBLOCK + I >*/
        m = j - (*nblock << 1) + i__;