X-Git-Url: http://vcs.maemo.org/git/?a=blobdiff_plain;f=3rdparty%2Flapack%2Fdlaebz.c;fp=3rdparty%2Flapack%2Fdlaebz.c;h=a76f635b2589c30015d7612a95154f93d10c6bff;hb=e4c14cdbdf2fe805e79cd96ded236f57e7b89060;hp=0000000000000000000000000000000000000000;hpb=454138ff8a20f6edb9b65a910101403d8b520643;p=opencv diff --git a/3rdparty/lapack/dlaebz.c b/3rdparty/lapack/dlaebz.c new file mode 100644 index 0000000..a76f635 --- /dev/null +++ b/3rdparty/lapack/dlaebz.c @@ -0,0 +1,627 @@ +#include "clapack.h" + +/* Subroutine */ int dlaebz_(integer *ijob, integer *nitmax, integer *n, + integer *mmax, integer *minp, integer *nbmin, doublereal *abstol, + doublereal *reltol, doublereal *pivmin, doublereal *d__, doublereal * + e, doublereal *e2, integer *nval, doublereal *ab, doublereal *c__, + integer *mout, integer *nab, doublereal *work, integer *iwork, + integer *info) +{ + /* System generated locals */ + integer nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4, + i__5, i__6; + doublereal d__1, d__2, d__3, d__4; + + /* Local variables */ + integer j, kf, ji, kl, jp, jit; + doublereal tmp1, tmp2; + integer itmp1, itmp2, kfnew, klnew; + + +/* -- LAPACK auxiliary routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DLAEBZ contains the iteration loops which compute and use the */ +/* function N(w), which is the count of eigenvalues of a symmetric */ +/* tridiagonal matrix T less than or equal to its argument w. It */ +/* performs a choice of two types of loops: */ + +/* IJOB=1, followed by */ +/* IJOB=2: It takes as input a list of intervals and returns a list of */ +/* sufficiently small intervals whose union contains the same */ +/* eigenvalues as the union of the original intervals. */ +/* The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. */ +/* The output interval (AB(j,1),AB(j,2)] will contain */ +/* eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. */ + +/* IJOB=3: It performs a binary search in each input interval */ +/* (AB(j,1),AB(j,2)] for a point w(j) such that */ +/* N(w(j))=NVAL(j), and uses C(j) as the starting point of */ +/* the search. If such a w(j) is found, then on output */ +/* AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output */ +/* (AB(j,1),AB(j,2)] will be a small interval containing the */ +/* point where N(w) jumps through NVAL(j), unless that point */ +/* lies outside the initial interval. */ + +/* Note that the intervals are in all cases half-open intervals, */ +/* i.e., of the form (a,b] , which includes b but not a . */ + +/* To avoid underflow, the matrix should be scaled so that its largest */ +/* element is no greater than overflow**(1/2) * underflow**(1/4) */ +/* in absolute value. To assure the most accurate computation */ +/* of small eigenvalues, the matrix should be scaled to be */ +/* not much smaller than that, either. */ + +/* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */ +/* Matrix", Report CS41, Computer Science Dept., Stanford */ +/* University, July 21, 1966 */ + +/* Note: the arguments are, in general, *not* checked for unreasonable */ +/* values. */ + +/* Arguments */ +/* ========= */ + +/* IJOB (input) INTEGER */ +/* Specifies what is to be done: */ +/* = 1: Compute NAB for the initial intervals. */ +/* = 2: Perform bisection iteration to find eigenvalues of T. */ +/* = 3: Perform bisection iteration to invert N(w), i.e., */ +/* to find a point which has a specified number of */ +/* eigenvalues of T to its left. */ +/* Other values will cause DLAEBZ to return with INFO=-1. */ + +/* NITMAX (input) INTEGER */ +/* The maximum number of "levels" of bisection to be */ +/* performed, i.e., an interval of width W will not be made */ +/* smaller than 2^(-NITMAX) * W. If not all intervals */ +/* have converged after NITMAX iterations, then INFO is set */ +/* to the number of non-converged intervals. */ + +/* N (input) INTEGER */ +/* The dimension n of the tridiagonal matrix T. It must be at */ +/* least 1. */ + +/* MMAX (input) INTEGER */ +/* The maximum number of intervals. If more than MMAX intervals */ +/* are generated, then DLAEBZ will quit with INFO=MMAX+1. */ + +/* MINP (input) INTEGER */ +/* The initial number of intervals. It may not be greater than */ +/* MMAX. */ + +/* NBMIN (input) INTEGER */ +/* The smallest number of intervals that should be processed */ +/* using a vector loop. If zero, then only the scalar loop */ +/* will be used. */ + +/* ABSTOL (input) DOUBLE PRECISION */ +/* The minimum (absolute) width of an interval. When an */ +/* interval is narrower than ABSTOL, or than RELTOL times the */ +/* larger (in magnitude) endpoint, then it is considered to be */ +/* sufficiently small, i.e., converged. This must be at least */ +/* zero. */ + +/* RELTOL (input) DOUBLE PRECISION */ +/* The minimum relative width of an interval. When an interval */ +/* is narrower than ABSTOL, or than RELTOL times the larger (in */ +/* magnitude) endpoint, then it is considered to be */ +/* sufficiently small, i.e., converged. Note: this should */ +/* always be at least radix*machine epsilon. */ + +/* PIVMIN (input) DOUBLE PRECISION */ +/* The minimum absolute value of a "pivot" in the Sturm */ +/* sequence loop. This *must* be at least max |e(j)**2| * */ +/* safe_min and at least safe_min, where safe_min is at least */ +/* the smallest number that can divide one without overflow. */ + +/* D (input) DOUBLE PRECISION array, dimension (N) */ +/* The diagonal elements of the tridiagonal matrix T. */ + +/* E (input) DOUBLE PRECISION array, dimension (N) */ +/* The offdiagonal elements of the tridiagonal matrix T in */ +/* positions 1 through N-1. E(N) is arbitrary. */ + +/* E2 (input) DOUBLE PRECISION array, dimension (N) */ +/* The squares of the offdiagonal elements of the tridiagonal */ +/* matrix T. E2(N) is ignored. */ + +/* NVAL (input/output) INTEGER array, dimension (MINP) */ +/* If IJOB=1 or 2, not referenced. */ +/* If IJOB=3, the desired values of N(w). The elements of NVAL */ +/* will be reordered to correspond with the intervals in AB. */ +/* Thus, NVAL(j) on output will not, in general be the same as */ +/* NVAL(j) on input, but it will correspond with the interval */ +/* (AB(j,1),AB(j,2)] on output. */ + +/* AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2) */ +/* The endpoints of the intervals. AB(j,1) is a(j), the left */ +/* endpoint of the j-th interval, and AB(j,2) is b(j), the */ +/* right endpoint of the j-th interval. The input intervals */ +/* will, in general, be modified, split, and reordered by the */ +/* calculation. */ + +/* C (input/output) DOUBLE PRECISION array, dimension (MMAX) */ +/* If IJOB=1, ignored. */ +/* If IJOB=2, workspace. */ +/* If IJOB=3, then on input C(j) should be initialized to the */ +/* first search point in the binary search. */ + +/* MOUT (output) INTEGER */ +/* If IJOB=1, the number of eigenvalues in the intervals. */ +/* If IJOB=2 or 3, the number of intervals output. */ +/* If IJOB=3, MOUT will equal MINP. */ + +/* NAB (input/output) INTEGER array, dimension (MMAX,2) */ +/* If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). */ +/* If IJOB=2, then on input, NAB(i,j) should be set. It must */ +/* satisfy the condition: */ +/* N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), */ +/* which means that in interval i only eigenvalues */ +/* NAB(i,1)+1,...,NAB(i,2) will be considered. Usually, */ +/* NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with */ +/* IJOB=1. */ +/* On output, NAB(i,j) will contain */ +/* max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of */ +/* the input interval that the output interval */ +/* (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the */ +/* the input values of NAB(k,1) and NAB(k,2). */ +/* If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), */ +/* unless N(w) > NVAL(i) for all search points w , in which */ +/* case NAB(i,1) will not be modified, i.e., the output */ +/* value will be the same as the input value (modulo */ +/* reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) */ +/* for all search points w , in which case NAB(i,2) will */ +/* not be modified. Normally, NAB should be set to some */ +/* distinctive value(s) before DLAEBZ is called. */ + +/* WORK (workspace) DOUBLE PRECISION array, dimension (MMAX) */ +/* Workspace. */ + +/* IWORK (workspace) INTEGER array, dimension (MMAX) */ +/* Workspace. */ + +/* INFO (output) INTEGER */ +/* = 0: All intervals converged. */ +/* = 1--MMAX: The last INFO intervals did not converge. */ +/* = MMAX+1: More than MMAX intervals were generated. */ + +/* Further Details */ +/* =============== */ + +/* This routine is intended to be called only by other LAPACK */ +/* routines, thus the interface is less user-friendly. It is intended */ +/* for two purposes: */ + +/* (a) finding eigenvalues. In this case, DLAEBZ should have one or */ +/* more initial intervals set up in AB, and DLAEBZ should be called */ +/* with IJOB=1. This sets up NAB, and also counts the eigenvalues. */ +/* Intervals with no eigenvalues would usually be thrown out at */ +/* this point. Also, if not all the eigenvalues in an interval i */ +/* are desired, NAB(i,1) can be increased or NAB(i,2) decreased. */ +/* For example, set NAB(i,1)=NAB(i,2)-1 to get the largest */ +/* eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX */ +/* no smaller than the value of MOUT returned by the call with */ +/* IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1 */ +/* through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the */ +/* tolerance specified by ABSTOL and RELTOL. */ + +/* (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). */ +/* In this case, start with a Gershgorin interval (a,b). Set up */ +/* AB to contain 2 search intervals, both initially (a,b). One */ +/* NVAL element should contain f-1 and the other should contain l */ +/* , while C should contain a and b, resp. NAB(i,1) should be -1 */ +/* and NAB(i,2) should be N+1, to flag an error if the desired */ +/* interval does not lie in (a,b). DLAEBZ is then called with */ +/* IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals -- */ +/* j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while */ +/* if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r */ +/* >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and */ +/* N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and */ +/* w(l-r)=...=w(l+k) are handled similarly. */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + +/* Check for Errors */ + + /* Parameter adjustments */ + nab_dim1 = *mmax; + nab_offset = 1 + nab_dim1; + nab -= nab_offset; + ab_dim1 = *mmax; + ab_offset = 1 + ab_dim1; + ab -= ab_offset; + --d__; + --e; + --e2; + --nval; + --c__; + --work; + --iwork; + + /* Function Body */ + *info = 0; + if (*ijob < 1 || *ijob > 3) { + *info = -1; + return 0; + } + +/* Initialize NAB */ + + if (*ijob == 1) { + +/* Compute the number of eigenvalues in the initial intervals. */ + + *mout = 0; +/* DIR$ NOVECTOR */ + i__1 = *minp; + for (ji = 1; ji <= i__1; ++ji) { + for (jp = 1; jp <= 2; ++jp) { + tmp1 = d__[1] - ab[ji + jp * ab_dim1]; + if (abs(tmp1) < *pivmin) { + tmp1 = -(*pivmin); + } + nab[ji + jp * nab_dim1] = 0; + if (tmp1 <= 0.) { + nab[ji + jp * nab_dim1] = 1; + } + + i__2 = *n; + for (j = 2; j <= i__2; ++j) { + tmp1 = d__[j] - e2[j - 1] / tmp1 - ab[ji + jp * ab_dim1]; + if (abs(tmp1) < *pivmin) { + tmp1 = -(*pivmin); + } + if (tmp1 <= 0.) { + ++nab[ji + jp * nab_dim1]; + } +/* L10: */ + } +/* L20: */ + } + *mout = *mout + nab[ji + (nab_dim1 << 1)] - nab[ji + nab_dim1]; +/* L30: */ + } + return 0; + } + +/* Initialize for loop */ + +/* KF and KL have the following meaning: */ +/* Intervals 1,...,KF-1 have converged. */ +/* Intervals KF,...,KL still need to be refined. */ + + kf = 1; + kl = *minp; + +/* If IJOB=2, initialize C. */ +/* If IJOB=3, use the user-supplied starting point. */ + + if (*ijob == 2) { + i__1 = *minp; + for (ji = 1; ji <= i__1; ++ji) { + c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5; +/* L40: */ + } + } + +/* 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 */ + d__1 = work[ji], d__2 = -(*pivmin); + work[ji] = min(d__1,d__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 */ + d__1 = work[ji], d__2 = -(*pivmin); + work[ji] = min(d__1,d__2); + } +/* L50: */ + } +/* L60: */ + } + + 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[ji + nab_dim1], i__6 = iwork[ji]; + i__3 = nab[ji + (nab_dim1 << 1)], 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[ji + (nab_dim1 << 1)]) { + +/* No eigenvalue in the upper interval: */ +/* just use the lower interval. */ + + ab[ji + (ab_dim1 << 1)] = c__[ji]; + + } else if (iwork[ji] == nab[ji + nab_dim1]) { + +/* No eigenvalue in the lower interval: */ +/* just use the upper interval. */ + + ab[ji + ab_dim1] = c__[ji]; + } else { + ++klnew; + if (klnew <= *mmax) { + +/* Eigenvalue in both intervals -- add upper to */ +/* queue. */ + + ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << + 1)]; + nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 + << 1)]; + ab[klnew + ab_dim1] = c__[ji]; + nab[klnew + nab_dim1] = iwork[ji]; + ab[ji + (ab_dim1 << 1)] = c__[ji]; + nab[ji + (nab_dim1 << 1)] = 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[ji + ab_dim1] = c__[ji]; + nab[ji + nab_dim1] = iwork[ji]; + } + if (iwork[ji] >= nval[ji]) { + ab[ji + (ab_dim1 << 1)] = c__[ji]; + nab[ji + (nab_dim1 << 1)] = 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 */ + d__1 = tmp2, d__2 = -(*pivmin); + tmp2 = min(d__1,d__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 */ + d__1 = tmp2, d__2 = -(*pivmin); + tmp2 = min(d__1,d__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 + nab_dim1]; + i__3 = nab[ji + (nab_dim1 << 1)], 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[ji + (nab_dim1 << 1)]) { + +/* No eigenvalue in the upper interval: */ +/* just use the lower interval. */ + + ab[ji + (ab_dim1 << 1)] = tmp1; + + } else if (itmp1 == nab[ji + nab_dim1]) { + +/* No eigenvalue in the lower interval: */ +/* just use the upper interval. */ + + ab[ji + ab_dim1] = tmp1; + } else if (klnew < *mmax) { + +/* Eigenvalue in both intervals -- add upper to queue. */ + + ++klnew; + ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << 1)]; + nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 << + 1)]; + ab[klnew + ab_dim1] = tmp1; + nab[klnew + nab_dim1] = itmp1; + ab[ji + (ab_dim1 << 1)] = tmp1; + nab[ji + (nab_dim1 << 1)] = 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[ji + ab_dim1] = tmp1; + nab[ji + nab_dim1] = itmp1; + } + if (itmp1 >= nval[ji]) { + ab[ji + (ab_dim1 << 1)] = tmp1; + nab[ji + (nab_dim1 << 1)] = itmp1; + } + } +/* L100: */ + } + 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 = (d__1 = ab[ji + (ab_dim1 << 1)] - ab[ji + ab_dim1], abs( + d__1)); +/* Computing MAX */ + d__3 = (d__1 = ab[ji + (ab_dim1 << 1)], abs(d__1)), d__4 = (d__2 = + ab[ji + ab_dim1], abs(d__2)); + tmp2 = max(d__3,d__4); +/* Computing MAX */ + d__1 = max(*abstol,*pivmin), d__2 = *reltol * tmp2; + if (tmp1 < max(d__1,d__2) || nab[ji + nab_dim1] >= nab[ji + ( + nab_dim1 << 1)]) { + +/* Converged -- Swap with position KFNEW, */ +/* then increment KFNEW */ + + if (ji > kfnew) { + tmp1 = ab[ji + ab_dim1]; + tmp2 = ab[ji + (ab_dim1 << 1)]; + itmp1 = nab[ji + nab_dim1]; + itmp2 = nab[ji + (nab_dim1 << 1)]; + ab[ji + ab_dim1] = ab[kfnew + ab_dim1]; + ab[ji + (ab_dim1 << 1)] = ab[kfnew + (ab_dim1 << 1)]; + nab[ji + nab_dim1] = nab[kfnew + nab_dim1]; + nab[ji + (nab_dim1 << 1)] = nab[kfnew + (nab_dim1 << 1)]; + ab[kfnew + ab_dim1] = tmp1; + ab[kfnew + (ab_dim1 << 1)] = tmp2; + nab[kfnew + nab_dim1] = itmp1; + nab[kfnew + (nab_dim1 << 1)] = 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[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5; +/* 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 DLAEBZ */ + +} /* dlaebz_ */