📄 slaebz.c
字号:
/* L40: */
}
/* Increment opcount for initializing C. */
latime_1.ops += *minp << 1;
}
/* Iteration loop */
i__1 = *nitmax;
for (jit = 1; jit <= i__1; ++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 <= i__2; ++ji) {
/* Compute N(c), the number of eigenvalues less than 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 = -(*pivmin);
work[ji] = dmin(r__1,r__2);
}
i__3 = *n;
for (j = 2; j <= i__3; ++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 = -(*pivmin);
work[ji] = dmin(r__1,r__2);
}
/* L50: */
}
/* L60: */
}
/* Increment iteration counter. */
latime_1.itcnt = latime_1.itcnt + kl - kf + 1;
/* Increment opcount for evaluating Sturm sequences on
each interval. */
latime_1.ops += (kl - kf + 1) * (*n - 1) * 3;
if (*ijob <= 2) {
/* IJOB=2: Choose all intervals containing eigenvalues. */
klnew = kl;
i__2 = kl;
for (ji = kf; ji <= i__2; ++ji) {
/* Insure that N(w) is monotone
Computing MIN
Computing MAX */
i__5 = nab_ref(ji, 1), i__6 = iwork[ji];
i__3 = nab_ref(ji, 2), i__4 = max(i__5,i__6);
iwork[ji] = min(i__3,i__4);
/* Update the Queue -- add intervals if both halves
contain eigenvalues. */
if (iwork[ji] == nab_ref(ji, 2)) {
/* No eigenvalue in the upper interval:
just use the lower interval. */
ab_ref(ji, 2) = c__[ji];
} else if (iwork[ji] == nab_ref(ji, 1)) {
/* No eigenvalue in the lower interval:
just use the upper interval. */
ab_ref(ji, 1) = c__[ji];
} else {
++klnew;
if (klnew <= *mmax) {
/* Eigenvalue in both intervals -- add upper to
queue. */
ab_ref(klnew, 2) = ab_ref(ji, 2);
nab_ref(klnew, 2) = nab_ref(ji, 2);
ab_ref(klnew, 1) = c__[ji];
nab_ref(klnew, 1) = iwork[ji];
ab_ref(ji, 2) = c__[ji];
nab_ref(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 <= i__2; ++ji) {
if (iwork[ji] <= nval[ji]) {
ab_ref(ji, 1) = c__[ji];
nab_ref(ji, 1) = iwork[ji];
}
if (iwork[ji] >= nval[ji]) {
ab_ref(ji, 2) = c__[ji];
nab_ref(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 <= i__2; ++ji) {
/* Compute N(w), the number of eigenvalues less than w */
tmp1 = c__[ji];
tmp2 = d__[1] - tmp1;
itmp1 = 0;
if (tmp2 <= *pivmin) {
itmp1 = 1;
/* Computing MIN */
r__1 = tmp2, r__2 = -(*pivmin);
tmp2 = dmin(r__1,r__2);
}
/* A series of compiler directives to defeat vectorization
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 <= i__3; ++j) {
tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1;
if (tmp2 <= *pivmin) {
++itmp1;
/* Computing MIN */
r__1 = tmp2, r__2 = -(*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_ref(ji, 1);
i__3 = nab_ref(ji, 2), i__4 = max(i__5,itmp1);
itmp1 = min(i__3,i__4);
/* Update the Queue -- add intervals if both halves
contain eigenvalues. */
if (itmp1 == nab_ref(ji, 2)) {
/* No eigenvalue in the upper interval:
just use the lower interval. */
ab_ref(ji, 2) = tmp1;
} else if (itmp1 == nab_ref(ji, 1)) {
/* No eigenvalue in the lower interval:
just use the upper interval. */
ab_ref(ji, 1) = tmp1;
} else if (klnew < *mmax) {
/* Eigenvalue in both intervals -- add upper to queue. */
++klnew;
ab_ref(klnew, 2) = ab_ref(ji, 2);
nab_ref(klnew, 2) = nab_ref(ji, 2);
ab_ref(klnew, 1) = tmp1;
nab_ref(klnew, 1) = itmp1;
ab_ref(ji, 2) = tmp1;
nab_ref(ji, 2) = itmp1;
} else {
*info = *mmax + 1;
return 0;
}
} else {
/* IJOB=3: Binary search. Keep only the interval
containing w s.t. N(w) = NVAL */
if (itmp1 <= nval[ji]) {
ab_ref(ji, 1) = tmp1;
nab_ref(ji, 1) = itmp1;
}
if (itmp1 >= nval[ji]) {
ab_ref(ji, 2) = tmp1;
nab_ref(ji, 2) = itmp1;
}
}
/* L100: */
}
/* Increment iteration counter. */
latime_1.itcnt = latime_1.itcnt + kl - kf + 1;
/* Increment opcount for evaluating Sturm sequences on
each interval. */
latime_1.ops += (kl - kf + 1) * (*n - 1) * 3;
kl = klnew;
/* End of Serial Version of the loop */
}
/* Check for convergence */
kfnew = kf;
i__2 = kl;
for (ji = kf; ji <= i__2; ++ji) {
tmp1 = (r__1 = ab_ref(ji, 2) - ab_ref(ji, 1), dabs(r__1));
/* Computing MAX */
r__3 = (r__1 = ab_ref(ji, 2), dabs(r__1)), r__4 = (r__2 = ab_ref(
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_ref(ji, 1) >= nab_ref(ji, 2)) {
/* Converged -- Swap with position KFNEW,
then increment KFNEW */
if (ji > kfnew) {
tmp1 = ab_ref(ji, 1);
tmp2 = ab_ref(ji, 2);
itmp1 = nab_ref(ji, 1);
itmp2 = nab_ref(ji, 2);
ab_ref(ji, 1) = ab_ref(kfnew, 1);
ab_ref(ji, 2) = ab_ref(kfnew, 2);
nab_ref(ji, 1) = nab_ref(kfnew, 1);
nab_ref(ji, 2) = nab_ref(kfnew, 2);
ab_ref(kfnew, 1) = tmp1;
ab_ref(kfnew, 2) = tmp2;
nab_ref(kfnew, 1) = itmp1;
nab_ref(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 <= i__2; ++ji) {
c__[ji] = (ab_ref(ji, 1) + ab_ref(ji, 2)) * .5f;
/* L120: */
}
/* Increment opcount for convergence check and choosing midpoints. */
latime_1.ops += kl - kf + 1 << 2;
/* 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_ */
#undef nab_ref
#undef ab_ref
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -