slarrv.c

著名的LAPACK矩阵计算软件包, 是比较新的版本, 一般用到矩阵分解的朋友也许会用到

	r__2 = (r__1 = work[in], dabs(r__1));
	work[indgap + in] = dmax(r__2,eps);
	ndone = 0;

	ndepth = 0;
	lsbdpt = 1;
	nclus = 1;
	iwork[iindc1 + 1] = 1;
	iwork[iindc1 + 2] = in;

/*        While( NDONE.LT.IM ) do */

L40:
	if (ndone < im) {
	    oldncl = nclus;
	    nclus = 0;
	    lsbdpt = 1 - lsbdpt;
	    i__2 = oldncl;
	    for (i__ = 1; i__ <= i__2; ++i__) {
		if (lsbdpt == 0) {
		    oldcls = iindc1;
		    newcls = iindc2;
		} else {
		    oldcls = iindc2;
		    newcls = iindc1;
		}

/*              If NDEPTH > 1, retrieve the relatively robust   
                representation (RRR) and perform limited bisection   
                (if necessary) to get approximate eigenvalues. */

		j = oldcls + (i__ << 1);
		oldfst = iwork[j - 1];
		oldlst = iwork[j];
		if (ndepth > 0) {
		    j = oldien + oldfst;
		    scopy_(&in, &z___ref(ibegin, j), &c__1, &d__[ibegin], &
			    c__1);
		    scopy_(&in, &z___ref(ibegin, j + 1), &c__1, &l[ibegin], &
			    c__1);
		    sigma = l[iend];
		}
		k = ibegin;
		latime_1.ops += (real) (in - 1 << 1);
		i__3 = in - 1;
		for (j = 1; j <= i__3; ++j) {
		    work[indld + j] = d__[k] * l[k];
		    work[indlld + j] = work[indld + j] * l[k];
		    ++k;
/* L50: */
		}
		if (ndepth > 0) {
		    slarrb_(&in, &d__[ibegin], &l[ibegin], &work[indld + 1], &
			    work[indlld + 1], &oldfst, &oldlst, &sigma, &
			    reltol, &work[1], &work[indgap + 1], &work[inderr]
			    , &work[indwrk], &iwork[iindwk], &iinfo);
		    if (iinfo != 0) {
			*info = 1;
			return 0;
		    }
		}

/*              Classify eigenvalues of the current representation (RRR)   
                as (i) isolated, (ii) loosely clustered or (iii) tightly   
                clustered */

		newfrs = oldfst;
		i__3 = oldlst;
		for (j = oldfst; j <= i__3; ++j) {
		    latime_1.ops += 1.f;
		    if (j == oldlst || work[indgap + j] >= reltol * (r__1 = 
			    work[j], dabs(r__1))) {
			newlst = j;
		    } else {

/*                    continue (to the next loop) */

			latime_1.ops += 1.f;
			relgap = work[indgap + j] / (r__1 = work[j], dabs(
				r__1));
			if (j == newfrs) {
			    minrgp = relgap;
			} else {
			    minrgp = dmin(minrgp,relgap);
			}
			goto L140;
		    }
		    newsiz = newlst - newfrs + 1;
		    maxitr = 10;
		    newftt = oldien + newfrs;
		    if (newsiz > 1) {
			mgscls = newsiz <= 20 && minrgp >= mgstol;
			if (! mgscls) {
			    slarrf_(&in, &d__[ibegin], &l[ibegin], &work[
				    indld + 1], &work[indlld + 1], &newfrs, &
				    newlst, &work[1], &z___ref(ibegin, newftt)
				    , &z___ref(ibegin, newftt + 1), &work[
				    indwrk], &iwork[iindwk], info);
			    if (*info == 0) {
				++nclus;
				k = newcls + (nclus << 1);
				iwork[k - 1] = newfrs;
				iwork[k] = newlst;
			    } else {
				*info = 0;
				if (minrgp >= mgstol) {
				    mgscls = TRUE_;
				} else {

/*                             Call SSTEIN to process this tight cluster.   
                               This happens only if MINRGP <= MGSTOL   
                               and SLARRF returns INFO = 1. The latter   
                               means that a new RRR to "break" the   
                               cluster could not be found. */

				    work[indwrk] = d__[ibegin];
				    latime_1.ops += (real) (in - 1);
				    i__4 = in - 1;
				    for (k = 1; k <= i__4; ++k) {
					work[indwrk + k] = d__[ibegin + k] + 
						work[indlld + k];
/* L60: */
				    }
				    i__4 = newsiz;
				    for (k = 1; k <= i__4; ++k) {
					iwork[iindwk + k - 1] = 1;
/* L70: */
				    }
				    i__4 = newlst;
				    for (k = newfrs; k <= i__4; ++k) {
					isuppz[(ibegin + k << 1) - 3] = 1;
					isuppz[(ibegin + k << 1) - 2] = in;
/* L80: */
				    }
				    temp[0] = in;
				    sstein_(&in, &work[indwrk], &work[indld + 
					    1], &newsiz, &work[newfrs], &
					    iwork[iindwk], temp, &z___ref(
					    ibegin, newftt), ldz, &work[
					    indwrk + in], &iwork[iindwk + in],
					     &iwork[iindwk + (in << 1)], &
					    iinfo);
				    if (iinfo != 0) {
					*info = 2;
					return 0;
				    }
				    ndone += newsiz;
				}
			    }
			}
		    } else {
			mgscls = FALSE_;
		    }
		    if (newsiz == 1 || mgscls) {
			ktot = newftt;
			i__4 = newlst;
			for (k = newfrs; k <= i__4; ++k) {
			    iter = 0;
L90:
			    lambda = work[k];
			    slar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &
				    l[ibegin], &work[indld + 1], &work[indlld 
				    + 1], &gersch[(oldien << 1) + 1], &
				    z___ref(ibegin, ktot), &ztz, &mingma, &
				    iwork[iindr + ktot], &isuppz[(ktot << 1) 
				    - 1], &work[indwrk]);
			    latime_1.ops += 4.f;
			    tmp1 = 1.f / ztz;
			    nrminv = sqrt(tmp1);
			    resid = dabs(mingma) * nrminv;
			    rqcorr = mingma * tmp1;
			    if (k == in) {
				gap = work[indgap + k - 1];
			    } else if (k == 1) {
				gap = work[indgap + k];
			    } else {
/* Computing MIN */
				r__1 = work[indgap + k - 1], r__2 = work[
					indgap + k];
				gap = dmin(r__1,r__2);
			    }
			    ++iter;
			    latime_1.ops += 3.f;
			    if (resid > *tol * gap && dabs(rqcorr) > eps * 
				    4.f * dabs(lambda)) {
				latime_1.ops += 1.f;
				work[k] = lambda + rqcorr;
				if (iter < maxitr) {
				    goto L90;
				}
			    }
			    iwork[ktot] = 1;
			    if (newsiz == 1) {
				++ndone;
			    }
			    latime_1.ops += (real) in;
			    sscal_(&in, &nrminv, &z___ref(ibegin, ktot), &
				    c__1);
			    ++ktot;
/* L100: */
			}
			if (newsiz > 1) {
			    itmp1 = isuppz[(newftt << 1) - 1];
			    itmp2 = isuppz[newftt * 2];
			    ktot = oldien + newlst;
			    i__4 = ktot;
			    for (p = newftt + 1; p <= i__4; ++p) {
				i__5 = p - 1;
				for (q = newftt; q <= i__5; ++q) {
				    latime_1.ops += (real) (in << 2);
				    tmp1 = -sdot_(&in, &z___ref(ibegin, p), &
					    c__1, &z___ref(ibegin, q), &c__1);
				    saxpy_(&in, &tmp1, &z___ref(ibegin, q), &
					    c__1, &z___ref(ibegin, p), &c__1);
/* L110: */
				}
				latime_1.ops += (real) (in * 3 + 1);
				tmp1 = 1.f / snrm2_(&in, &z___ref(ibegin, p), 
					&c__1);
				sscal_(&in, &tmp1, &z___ref(ibegin, p), &c__1)
					;
/* Computing MIN */
				i__5 = itmp1, i__6 = isuppz[(p << 1) - 1];
				itmp1 = min(i__5,i__6);
/* Computing MAX */
				i__5 = itmp2, i__6 = isuppz[p * 2];
				itmp2 = max(i__5,i__6);
/* L120: */
			    }
			    i__4 = ktot;
			    for (p = newftt; p <= i__4; ++p) {
				isuppz[(p << 1) - 1] = itmp1;
				isuppz[p * 2] = itmp2;
/* L130: */
			    }
			    ndone += newsiz;
			}
		    }
		    newfrs = j + 1;
L140:
		    ;
		}
/* L150: */
	    }
	    ++ndepth;
	    goto L40;
	}
	j = ibegin << 1;
	i__2 = iend;
	for (i__ = ibegin; i__ <= i__2; ++i__) {
	    isuppz[j - 1] += oldien;
	    isuppz[j] += oldien;
	    j += 2;
/* L160: */
	}
	ibegin = iend + 1;
L170:
	;
    }

    return 0;

/*     End of SLARRV */

} /* slarrv_ */

#undef z___ref