slaebz.c

来自「NIST Handwriting OCR Testbed」· C语言 代码 · 共 648 行 · 第 1/2 页

    if (*ijob == 2) {
	i__1 = *minp;
	for (ji = 1; ji <= *minp; ++ji) {
	    C(ji) = (AB(ji,1) + AB(ji,2)) * .5f;
/* L40: */
	}
    }

/*     Iteration loop */

    i__1 = *nitmax;
    for (jit = 1; jit <= *nitmax; ++jit) {

/*        Loop over intervals */

	if (kl - kf + 1 >= *nbmin && *nbmin > 0) {

/*           Begin of Parallel Version of the loop */

	    i__2 = kl;
	    for (ji = kf; ji <= kl; ++ji) {

/*              Compute N(c), the number of eigenvalues less t
han c */

		WORK(ji) = D(1) - C(ji);
		IWORK(ji) = 0;
		if (WORK(ji) <= *pivmin) {
		    IWORK(ji) = 1;
/* Computing MIN */
		    r__1 = WORK(ji), r__2 = -(doublereal)(*pivmin);
		    WORK(ji) = dmin(r__1,r__2);
		}

		i__3 = *n;
		for (j = 2; j <= *n; ++j) {
		    WORK(ji) = D(j) - E2(j - 1) / WORK(ji) - C(ji);
		    if (WORK(ji) <= *pivmin) {
			++IWORK(ji);
/* Computing MIN */
			r__1 = WORK(ji), r__2 = -(doublereal)(*pivmin);
			WORK(ji) = dmin(r__1,r__2);
		    }
/* L50: */
		}
/* L60: */
	    }

	    if (*ijob <= 2) {

/*              IJOB=2: Choose all intervals containing eigenv
alues. */

		klnew = kl;
		i__2 = kl;
		for (ji = kf; ji <= kl; ++ji) {

/*                 Insure that N(w) is monotone   

   Computing MIN   
   Computing MAX */
		    i__5 = NAB(ji,1), i__6 = IWORK(ji);
		    i__3 = NAB(ji,2), i__4 = max(i__5,i__6);
		    IWORK(ji) = min(i__3,i__4);

/*                 Update the Queue -- add intervals if bo
th halves   
                   contain eigenvalues. */

		    if (IWORK(ji) == NAB(ji,2)) {

/*                    No eigenvalue in the upper inter
val:   
                      just use the lower interval. */

			AB(ji,2) = C(ji);

		    } else if (IWORK(ji) == NAB(ji,1)) {

/*                    No eigenvalue in the lower inter
val:   
                      just use the upper interval. */

			AB(ji,1) = C(ji);
		    } else {
			++klnew;
			if (klnew <= *mmax) {

/*                       Eigenvalue in both interv
als -- add upper to   
                         queue. */

			    AB(klnew,2) = AB(ji,2);
			    NAB(klnew,2) = NAB(ji,2);
			    AB(klnew,1) = C(ji);
			    NAB(klnew,1) = IWORK(ji);
			    AB(ji,2) = C(ji);
			    NAB(ji,2) = IWORK(ji);
			} else {
			    *info = *mmax + 1;
			}
		    }
/* L70: */
		}
		if (*info != 0) {
		    return 0;
		}
		kl = klnew;
	    } else {

/*              IJOB=3: Binary search.  Keep only the interval
 containing   
                        w   s.t. N(w) = NVAL */

		i__2 = kl;
		for (ji = kf; ji <= kl; ++ji) {
		    if (IWORK(ji) <= NVAL(ji)) {
			AB(ji,1) = C(ji);
			NAB(ji,1) = IWORK(ji);
		    }
		    if (IWORK(ji) >= NVAL(ji)) {
			AB(ji,2) = C(ji);
			NAB(ji,2) = IWORK(ji);
		    }
/* L80: */
		}
	    }

	} else {

/*           End of Parallel Version of the loop   

             Begin of Serial Version of the loop */

	    klnew = kl;
	    i__2 = kl;
	    for (ji = kf; ji <= kl; ++ji) {

/*              Compute N(w), the number of eigenvalues less t
han w */

		tmp1 = C(ji);
		tmp2 = D(1) - tmp1;
		itmp1 = 0;
		if (tmp2 <= *pivmin) {
		    itmp1 = 1;
/* Computing MIN */
		    r__1 = tmp2, r__2 = -(doublereal)(*pivmin);
		    tmp2 = dmin(r__1,r__2);
		}

/*              A series of compiler directives to defeat vect
orization   
                for the next loop   

   $PL$ CMCHAR=' '   
   DIR$          NEXTSCALAR   
   $DIR          SCALAR   
   DIR$          NEXT SCALAR   
   VD$L          NOVECTOR   
   DEC$          NOVECTOR   
   VD$           NOVECTOR   
   VDIR          NOVECTOR   
   VOCL          LOOP,SCALAR   
   IBM           PREFER SCALAR   
   $PL$ CMCHAR='*' */

		i__3 = *n;
		for (j = 2; j <= *n; ++j) {
		    tmp2 = D(j) - E2(j - 1) / tmp2 - tmp1;
		    if (tmp2 <= *pivmin) {
			++itmp1;
/* Computing MIN */
			r__1 = tmp2, r__2 = -(doublereal)(*pivmin);
			tmp2 = dmin(r__1,r__2);
		    }
/* L90: */
		}

		if (*ijob <= 2) {

/*                 IJOB=2: Choose all intervals containing
 eigenvalues.   

                   Insure that N(w) is monotone   

   Computing MIN   
   Computing MAX */
		    i__5 = NAB(ji,1);
		    i__3 = NAB(ji,2), i__4 = max(i__5,itmp1);
		    itmp1 = min(i__3,i__4);

/*                 Update the Queue -- add intervals if bo
th halves   
                   contain eigenvalues. */

		    if (itmp1 == NAB(ji,2)) {

/*                    No eigenvalue in the upper inter
val:   
                      just use the lower interval. */

			AB(ji,2) = tmp1;

		    } else if (itmp1 == NAB(ji,1)) {

/*                    No eigenvalue in the lower inter
val:   
                      just use the upper interval. */

			AB(ji,1) = tmp1;
		    } else if (klnew < *mmax) {

/*                    Eigenvalue in both intervals -- 
add upper to queue. */

			++klnew;
			AB(klnew,2) = AB(ji,2);
			NAB(klnew,2) = NAB(ji,2);
			AB(klnew,1) = tmp1;
			NAB(klnew,1) = itmp1;
			AB(ji,2) = tmp1;
			NAB(ji,2) = itmp1;
		    } else {
			*info = *mmax + 1;
			return 0;
		    }
		} else {

/*                 IJOB=3: Binary search.  Keep only the i
nterval   
                           containing  w  s.t. N(w) = NVAL
 */

		    if (itmp1 <= NVAL(ji)) {
			AB(ji,1) = tmp1;
			NAB(ji,1) = itmp1;
		    }
		    if (itmp1 >= NVAL(ji)) {
			AB(ji,2) = tmp1;
			NAB(ji,2) = itmp1;
		    }
		}
/* L100: */
	    }
	    kl = klnew;

/*           End of Serial Version of the loop */

	}

/*        Check for convergence */

	kfnew = kf;
	i__2 = kl;
	for (ji = kf; ji <= kl; ++ji) {
	    tmp1 = (r__1 = AB(ji,2) - AB(ji,1), dabs(
		    r__1));
/* Computing MAX */
	    r__3 = (r__1 = AB(ji,2), dabs(r__1)), r__4 = (r__2 
		    = AB(ji,1), dabs(r__2));
	    tmp2 = dmax(r__3,r__4);
/* Computing MAX */
	    r__1 = max(*abstol,*pivmin), r__2 = *reltol * tmp2;
	    if (tmp1 < dmax(r__1,r__2) || NAB(ji,1) >= NAB(ji,2)) {

/*              Converged -- Swap with position KFNEW,   
                             then increment KFNEW */

		if (ji > kfnew) {
		    tmp1 = AB(ji,1);
		    tmp2 = AB(ji,2);
		    itmp1 = NAB(ji,1);
		    itmp2 = NAB(ji,2);
		    AB(ji,1) = AB(kfnew,1);
		    AB(ji,2) = AB(kfnew,2);
		    NAB(ji,1) = NAB(kfnew,1);
		    NAB(ji,2) = NAB(kfnew,2);
		    AB(kfnew,1) = tmp1;
		    AB(kfnew,2) = tmp2;
		    NAB(kfnew,1) = itmp1;
		    NAB(kfnew,2) = itmp2;
		    if (*ijob == 3) {
			itmp1 = NVAL(ji);
			NVAL(ji) = NVAL(kfnew);
			NVAL(kfnew) = itmp1;
		    }
		}
		++kfnew;
	    }
/* L110: */
	}
	kf = kfnew;

/*        Choose Midpoints */

	i__2 = kl;
	for (ji = kf; ji <= kl; ++ji) {
	    C(ji) = (AB(ji,1) + AB(ji,2)) * .5f;
/* L120: */
	}

/*        If no more intervals to refine, quit. */

	if (kf > kl) {
	    goto L140;
	}
/* L130: */
    }

/*     Converged */

L140:
/* Computing MAX */
    i__1 = kl + 1 - kf;
    *info = max(i__1,0);
    *mout = kl;

    return 0;

/*     End of SLAEBZ */

} /* slaebz_ */