--- /dev/null
+#include "clapack.h"
+
+/* Table of constant values */
+
+static integer c__2 = 2;
+static integer c__1 = 1;
+static integer c_n1 = -1;
+
+/* Subroutine */ int dstein_(integer *n, doublereal *d__, doublereal *e,
+ integer *m, doublereal *w, integer *iblock, integer *isplit,
+ doublereal *z__, integer *ldz, doublereal *work, integer *iwork,
+ integer *ifail, integer *info)
+{
+ /* System generated locals */
+ integer z_dim1, z_offset, i__1, i__2, i__3;
+ doublereal d__1, d__2, d__3, d__4, d__5;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ integer i__, j, b1, j1, bn;
+ doublereal xj, scl, eps, sep, nrm, tol;
+ integer its;
+ doublereal xjm, ztr, eps1;
+ integer jblk, nblk;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ integer jmax;
+ extern doublereal dnrm2_(integer *, doublereal *, integer *);
+ extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
+ integer *);
+ integer iseed[4], gpind, iinfo;
+ extern doublereal dasum_(integer *, doublereal *, integer *);
+ extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
+ doublereal *, integer *), daxpy_(integer *, doublereal *,
+ doublereal *, integer *, doublereal *, integer *);
+ doublereal ortol;
+ integer indrv1, indrv2, indrv3, indrv4, indrv5;
+ extern doublereal dlamch_(char *);
+ extern /* Subroutine */ int dlagtf_(integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *, integer *
+, integer *);
+ extern integer idamax_(integer *, doublereal *, integer *);
+ extern /* Subroutine */ int xerbla_(char *, integer *), dlagts_(
+ integer *, integer *, doublereal *, doublereal *, doublereal *,
+ doublereal *, integer *, doublereal *, doublereal *, integer *);
+ integer nrmchk;
+ extern /* Subroutine */ int dlarnv_(integer *, integer *, integer *,
+ doublereal *);
+ integer blksiz;
+ doublereal onenrm, dtpcrt, pertol;
+
+
+/* -- LAPACK routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+/* .. Array Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* DSTEIN computes the eigenvectors of a real symmetric tridiagonal */
+/* matrix T corresponding to specified eigenvalues, using inverse */
+/* iteration. */
+
+/* The maximum number of iterations allowed for each eigenvector is */
+/* specified by an internal parameter MAXITS (currently set to 5). */
+
+/* Arguments */
+/* ========= */
+
+/* N (input) INTEGER */
+/* The order of the matrix. N >= 0. */
+
+/* D (input) DOUBLE PRECISION array, dimension (N) */
+/* The n diagonal elements of the tridiagonal matrix T. */
+
+/* E (input) DOUBLE PRECISION array, dimension (N-1) */
+/* The (n-1) subdiagonal elements of the tridiagonal matrix */
+/* T, in elements 1 to N-1. */
+
+/* M (input) INTEGER */
+/* The number of eigenvectors to be found. 0 <= M <= N. */
+
+/* W (input) DOUBLE PRECISION array, dimension (N) */
+/* The first M elements of W contain the eigenvalues for */
+/* which eigenvectors are to be computed. The eigenvalues */
+/* should be grouped by split-off block and ordered from */
+/* smallest to largest within the block. ( The output array */
+/* W from DSTEBZ with ORDER = 'B' is expected here. ) */
+
+/* IBLOCK (input) INTEGER array, dimension (N) */
+/* The submatrix indices associated with the corresponding */
+/* eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to */
+/* the first submatrix from the top, =2 if W(i) belongs to */
+/* the second submatrix, etc. ( The output array IBLOCK */
+/* from DSTEBZ is expected here. ) */
+
+/* ISPLIT (input) INTEGER array, dimension (N) */
+/* The splitting points, at which T breaks up into submatrices. */
+/* The first submatrix consists of rows/columns 1 to */
+/* ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
+/* through ISPLIT( 2 ), etc. */
+/* ( The output array ISPLIT from DSTEBZ is expected here. ) */
+
+/* Z (output) DOUBLE PRECISION array, dimension (LDZ, M) */
+/* The computed eigenvectors. The eigenvector associated */
+/* with the eigenvalue W(i) is stored in the i-th column of */
+/* Z. Any vector which fails to converge is set to its current */
+/* iterate after MAXITS iterations. */
+
+/* LDZ (input) INTEGER */
+/* The leading dimension of the array Z. LDZ >= max(1,N). */
+
+/* WORK (workspace) DOUBLE PRECISION array, dimension (5*N) */
+
+/* IWORK (workspace) INTEGER array, dimension (N) */
+
+/* IFAIL (output) INTEGER array, dimension (M) */
+/* On normal exit, all elements of IFAIL are zero. */
+/* If one or more eigenvectors fail to converge after */
+/* MAXITS iterations, then their indices are stored in */
+/* array IFAIL. */
+
+/* INFO (output) INTEGER */
+/* = 0: successful exit. */
+/* < 0: if INFO = -i, the i-th argument had an illegal value */
+/* > 0: if INFO = i, then i eigenvectors failed to converge */
+/* in MAXITS iterations. Their indices are stored in */
+/* array IFAIL. */
+
+/* Internal Parameters */
+/* =================== */
+
+/* MAXITS INTEGER, default = 5 */
+/* The maximum number of iterations performed. */
+
+/* EXTRA INTEGER, default = 2 */
+/* The number of iterations performed after norm growth */
+/* criterion is satisfied, should be at least 1. */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. Local Arrays .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. External Subroutines .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+/* Test the input parameters. */
+
+ /* Parameter adjustments */
+ --d__;
+ --e;
+ --w;
+ --iblock;
+ --isplit;
+ z_dim1 = *ldz;
+ z_offset = 1 + z_dim1;
+ z__ -= z_offset;
+ --work;
+ --iwork;
+ --ifail;
+
+ /* Function Body */
+ *info = 0;
+ i__1 = *m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ ifail[i__] = 0;
+/* L10: */
+ }
+
+ if (*n < 0) {
+ *info = -1;
+ } else if (*m < 0 || *m > *n) {
+ *info = -4;
+ } else if (*ldz < max(1,*n)) {
+ *info = -9;
+ } else {
+ i__1 = *m;
+ for (j = 2; j <= i__1; ++j) {
+ if (iblock[j] < iblock[j - 1]) {
+ *info = -6;
+ goto L30;
+ }
+ if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
+ *info = -5;
+ goto L30;
+ }
+/* L20: */
+ }
+L30:
+ ;
+ }
+
+ if (*info != 0) {
+ i__1 = -(*info);
+ xerbla_("DSTEIN", &i__1);
+ return 0;
+ }
+
+/* Quick return if possible */
+
+ if (*n == 0 || *m == 0) {
+ return 0;
+ } else if (*n == 1) {
+ z__[z_dim1 + 1] = 1.;
+ return 0;
+ }
+
+/* Get machine constants. */
+
+ eps = dlamch_("Precision");
+
+/* Initialize seed for random number generator DLARNV. */
+
+ for (i__ = 1; i__ <= 4; ++i__) {
+ iseed[i__ - 1] = 1;
+/* L40: */
+ }
+
+/* Initialize pointers. */
+
+ indrv1 = 0;
+ indrv2 = indrv1 + *n;
+ indrv3 = indrv2 + *n;
+ indrv4 = indrv3 + *n;
+ indrv5 = indrv4 + *n;
+
+/* Compute eigenvectors of matrix blocks. */
+
+ j1 = 1;
+ i__1 = iblock[*m];
+ for (nblk = 1; nblk <= i__1; ++nblk) {
+
+/* Find starting and ending indices of block nblk. */
+
+ if (nblk == 1) {
+ b1 = 1;
+ } else {
+ b1 = isplit[nblk - 1] + 1;
+ }
+ bn = isplit[nblk];
+ blksiz = bn - b1 + 1;
+ if (blksiz == 1) {
+ goto L60;
+ }
+ gpind = b1;
+
+/* Compute reorthogonalization criterion and stopping criterion. */
+
+ onenrm = (d__1 = d__[b1], abs(d__1)) + (d__2 = e[b1], abs(d__2));
+/* Computing MAX */
+ d__3 = onenrm, d__4 = (d__1 = d__[bn], abs(d__1)) + (d__2 = e[bn - 1],
+ abs(d__2));
+ onenrm = max(d__3,d__4);
+ i__2 = bn - 1;
+ for (i__ = b1 + 1; i__ <= i__2; ++i__) {
+/* Computing MAX */
+ d__4 = onenrm, d__5 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[
+ i__ - 1], abs(d__2)) + (d__3 = e[i__], abs(d__3));
+ onenrm = max(d__4,d__5);
+/* L50: */
+ }
+ ortol = onenrm * .001;
+
+ dtpcrt = sqrt(.1 / blksiz);
+
+/* Loop through eigenvalues of block nblk. */
+
+L60:
+ jblk = 0;
+ i__2 = *m;
+ for (j = j1; j <= i__2; ++j) {
+ if (iblock[j] != nblk) {
+ j1 = j;
+ goto L160;
+ }
+ ++jblk;
+ xj = w[j];
+
+/* Skip all the work if the block size is one. */
+
+ if (blksiz == 1) {
+ work[indrv1 + 1] = 1.;
+ goto L120;
+ }
+
+/* If eigenvalues j and j-1 are too close, add a relatively */
+/* small perturbation. */
+
+ if (jblk > 1) {
+ eps1 = (d__1 = eps * xj, abs(d__1));
+ pertol = eps1 * 10.;
+ sep = xj - xjm;
+ if (sep < pertol) {
+ xj = xjm + pertol;
+ }
+ }
+
+ its = 0;
+ nrmchk = 0;
+
+/* Get random starting vector. */
+
+ dlarnv_(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
+
+/* Copy the matrix T so it won't be destroyed in factorization. */
+
+ dcopy_(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
+ i__3 = blksiz - 1;
+ dcopy_(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
+ i__3 = blksiz - 1;
+ dcopy_(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
+
+/* Compute LU factors with partial pivoting ( PT = LU ) */
+
+ tol = 0.;
+ dlagtf_(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
+ indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
+
+/* Update iteration count. */
+
+L70:
+ ++its;
+ if (its > 5) {
+ goto L100;
+ }
+
+/* Normalize and scale the righthand side vector Pb. */
+
+/* Computing MAX */
+ d__2 = eps, d__3 = (d__1 = work[indrv4 + blksiz], abs(d__1));
+ scl = blksiz * onenrm * max(d__2,d__3) / dasum_(&blksiz, &work[
+ indrv1 + 1], &c__1);
+ dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
+
+/* Solve the system LU = Pb. */
+
+ dlagts_(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
+ work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
+ indrv1 + 1], &tol, &iinfo);
+
+/* Reorthogonalize by modified Gram-Schmidt if eigenvalues are */
+/* close enough. */
+
+ if (jblk == 1) {
+ goto L90;
+ }
+ if ((d__1 = xj - xjm, abs(d__1)) > ortol) {
+ gpind = j;
+ }
+ if (gpind != j) {
+ i__3 = j - 1;
+ for (i__ = gpind; i__ <= i__3; ++i__) {
+ ztr = -ddot_(&blksiz, &work[indrv1 + 1], &c__1, &z__[b1 +
+ i__ * z_dim1], &c__1);
+ daxpy_(&blksiz, &ztr, &z__[b1 + i__ * z_dim1], &c__1, &
+ work[indrv1 + 1], &c__1);
+/* L80: */
+ }
+ }
+
+/* Check the infinity norm of the iterate. */
+
+L90:
+ jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
+ nrm = (d__1 = work[indrv1 + jmax], abs(d__1));
+
+/* Continue for additional iterations after norm reaches */
+/* stopping criterion. */
+
+ if (nrm < dtpcrt) {
+ goto L70;
+ }
+ ++nrmchk;
+ if (nrmchk < 3) {
+ goto L70;
+ }
+
+ goto L110;
+
+/* If stopping criterion was not satisfied, update info and */
+/* store eigenvector number in array ifail. */
+
+L100:
+ ++(*info);
+ ifail[*info] = j;
+
+/* Accept iterate as jth eigenvector. */
+
+L110:
+ scl = 1. / dnrm2_(&blksiz, &work[indrv1 + 1], &c__1);
+ jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
+ if (work[indrv1 + jmax] < 0.) {
+ scl = -scl;
+ }
+ dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
+L120:
+ i__3 = *n;
+ for (i__ = 1; i__ <= i__3; ++i__) {
+ z__[i__ + j * z_dim1] = 0.;
+/* L130: */
+ }
+ i__3 = blksiz;
+ for (i__ = 1; i__ <= i__3; ++i__) {
+ z__[b1 + i__ - 1 + j * z_dim1] = work[indrv1 + i__];
+/* L140: */
+ }
+
+/* Save the shift to check eigenvalue spacing at next */
+/* iteration. */
+
+ xjm = xj;
+
+/* L150: */
+ }
+L160:
+ ;
+ }
+
+ return 0;
+
+/* End of DSTEIN */
+
+} /* dstein_ */
projects / opencv / blobdiff