git.maemo.org Git - opencv/blob - 3rdparty/lapack/dstein.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__2 = 2;
   6 static integer c__1 = 1;
   7 static integer c_n1 = -1;
   8
   9 /* Subroutine */ int dstein_(integer *n, doublereal *d__, doublereal *e,
  10         integer *m, doublereal *w, integer *iblock, integer *isplit,
  11         doublereal *z__, integer *ldz, doublereal *work, integer *iwork,
  12         integer *ifail, integer *info)
  13 {
  14     /* System generated locals */
  15     integer z_dim1, z_offset, i__1, i__2, i__3;
  16     doublereal d__1, d__2, d__3, d__4, d__5;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal);
  20
  21     /* Local variables */
  22     integer i__, j, b1, j1, bn;
  23     doublereal xj, scl, eps, sep, nrm, tol;
  24     integer its;
  25     doublereal xjm, ztr, eps1;
  26     integer jblk, nblk;
  27     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  28             integer *);
  29     integer jmax;
  30     extern doublereal dnrm2_(integer *, doublereal *, integer *);
  31     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  32             integer *);
  33     integer iseed[4], gpind, iinfo;
  34     extern doublereal dasum_(integer *, doublereal *, integer *);
  35     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  36             doublereal *, integer *), daxpy_(integer *, doublereal *,
  37             doublereal *, integer *, doublereal *, integer *);
  38     doublereal ortol;
  39     integer indrv1, indrv2, indrv3, indrv4, indrv5;
  40     extern doublereal dlamch_(char *);
  41     extern /* Subroutine */ int dlagtf_(integer *, doublereal *, doublereal *,
  42              doublereal *, doublereal *, doublereal *, doublereal *, integer *
  43 , integer *);
  44     extern integer idamax_(integer *, doublereal *, integer *);
  45     extern /* Subroutine */ int xerbla_(char *, integer *), dlagts_(
  46             integer *, integer *, doublereal *, doublereal *, doublereal *,
  47             doublereal *, integer *, doublereal *, doublereal *, integer *);
  48     integer nrmchk;
  49     extern /* Subroutine */ int dlarnv_(integer *, integer *, integer *,
  50             doublereal *);
  51     integer blksiz;
  52     doublereal onenrm, dtpcrt, pertol;
  53
  54
  55 /*  -- LAPACK routine (version 3.1) -- */
  56 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  57 /*     November 2006 */
  58
  59 /*     .. Scalar Arguments .. */
  60 /*     .. */
  61 /*     .. Array Arguments .. */
  62 /*     .. */
  63
  64 /*  Purpose */
  65 /*  ======= */
  66
  67 /*  DSTEIN computes the eigenvectors of a real symmetric tridiagonal */
  68 /*  matrix T corresponding to specified eigenvalues, using inverse */
  69 /*  iteration. */
  70
  71 /*  The maximum number of iterations allowed for each eigenvector is */
  72 /*  specified by an internal parameter MAXITS (currently set to 5). */
  73
  74 /*  Arguments */
  75 /*  ========= */
  76
  77 /*  N       (input) INTEGER */
  78 /*          The order of the matrix.  N >= 0. */
  79
  80 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
  81 /*          The n diagonal elements of the tridiagonal matrix T. */
  82
  83 /*  E       (input) DOUBLE PRECISION array, dimension (N-1) */
  84 /*          The (n-1) subdiagonal elements of the tridiagonal matrix */
  85 /*          T, in elements 1 to N-1. */
  86
  87 /*  M       (input) INTEGER */
  88 /*          The number of eigenvectors to be found.  0 <= M <= N. */
  89
  90 /*  W       (input) DOUBLE PRECISION array, dimension (N) */
  91 /*          The first M elements of W contain the eigenvalues for */
  92 /*          which eigenvectors are to be computed.  The eigenvalues */
  93 /*          should be grouped by split-off block and ordered from */
  94 /*          smallest to largest within the block.  ( The output array */
  95 /*          W from DSTEBZ with ORDER = 'B' is expected here. ) */
  96
  97 /*  IBLOCK  (input) INTEGER array, dimension (N) */
  98 /*          The submatrix indices associated with the corresponding */
  99 /*          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to */
 100 /*          the first submatrix from the top, =2 if W(i) belongs to */
 101 /*          the second submatrix, etc.  ( The output array IBLOCK */
 102 /*          from DSTEBZ is expected here. ) */
 103
 104 /*  ISPLIT  (input) INTEGER array, dimension (N) */
 105 /*          The splitting points, at which T breaks up into submatrices. */
 106 /*          The first submatrix consists of rows/columns 1 to */
 107 /*          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
 108 /*          through ISPLIT( 2 ), etc. */
 109 /*          ( The output array ISPLIT from DSTEBZ is expected here. ) */
 110
 111 /*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, M) */
 112 /*          The computed eigenvectors.  The eigenvector associated */
 113 /*          with the eigenvalue W(i) is stored in the i-th column of */
 114 /*          Z.  Any vector which fails to converge is set to its current */
 115 /*          iterate after MAXITS iterations. */
 116
 117 /*  LDZ     (input) INTEGER */
 118 /*          The leading dimension of the array Z.  LDZ >= max(1,N). */
 119
 120 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (5*N) */
 121
 122 /*  IWORK   (workspace) INTEGER array, dimension (N) */
 123
 124 /*  IFAIL   (output) INTEGER array, dimension (M) */
 125 /*          On normal exit, all elements of IFAIL are zero. */
 126 /*          If one or more eigenvectors fail to converge after */
 127 /*          MAXITS iterations, then their indices are stored in */
 128 /*          array IFAIL. */
 129
 130 /*  INFO    (output) INTEGER */
 131 /*          = 0: successful exit. */
 132 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
 133 /*          > 0: if INFO = i, then i eigenvectors failed to converge */
 134 /*               in MAXITS iterations.  Their indices are stored in */
 135 /*               array IFAIL. */
 136
 137 /*  Internal Parameters */
 138 /*  =================== */
 139
 140 /*  MAXITS  INTEGER, default = 5 */
 141 /*          The maximum number of iterations performed. */
 142
 143 /*  EXTRA   INTEGER, default = 2 */
 144 /*          The number of iterations performed after norm growth */
 145 /*          criterion is satisfied, should be at least 1. */
 146
 147 /*  ===================================================================== */
 148
 149 /*     .. Parameters .. */
 150 /*     .. */
 151 /*     .. Local Scalars .. */
 152 /*     .. */
 153 /*     .. Local Arrays .. */
 154 /*     .. */
 155 /*     .. External Functions .. */
 156 /*     .. */
 157 /*     .. External Subroutines .. */
 158 /*     .. */
 159 /*     .. Intrinsic Functions .. */
 160 /*     .. */
 161 /*     .. Executable Statements .. */
 162
 163 /*     Test the input parameters. */
 164
 165     /* Parameter adjustments */
 166     --d__;
 167     --e;
 168     --w;
 169     --iblock;
 170     --isplit;
 171     z_dim1 = *ldz;
 172     z_offset = 1 + z_dim1;
 173     z__ -= z_offset;
 174     --work;
 175     --iwork;
 176     --ifail;
 177
 178     /* Function Body */
 179     *info = 0;
 180     i__1 = *m;
 181     for (i__ = 1; i__ <= i__1; ++i__) {
 182         ifail[i__] = 0;
 183 /* L10: */
 184     }
 185
 186     if (*n < 0) {
 187         *info = -1;
 188     } else if (*m < 0 || *m > *n) {
 189         *info = -4;
 190     } else if (*ldz < max(1,*n)) {
 191         *info = -9;
 192     } else {
 193         i__1 = *m;
 194         for (j = 2; j <= i__1; ++j) {
 195             if (iblock[j] < iblock[j - 1]) {
 196                 *info = -6;
 197                 goto L30;
 198             }
 199             if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
 200                 *info = -5;
 201                 goto L30;
 202             }
 203 /* L20: */
 204         }
 205 L30:
 206         ;
 207     }
 208
 209     if (*info != 0) {
 210         i__1 = -(*info);
 211         xerbla_("DSTEIN", &i__1);
 212         return 0;
 213     }
 214
 215 /*     Quick return if possible */
 216
 217     if (*n == 0 || *m == 0) {
 218         return 0;
 219     } else if (*n == 1) {
 220         z__[z_dim1 + 1] = 1.;
 221         return 0;
 222     }
 223
 224 /*     Get machine constants. */
 225
 226     eps = dlamch_("Precision");
 227
 228 /*     Initialize seed for random number generator DLARNV. */
 229
 230     for (i__ = 1; i__ <= 4; ++i__) {
 231         iseed[i__ - 1] = 1;
 232 /* L40: */
 233     }
 234
 235 /*     Initialize pointers. */
 236
 237     indrv1 = 0;
 238     indrv2 = indrv1 + *n;
 239     indrv3 = indrv2 + *n;
 240     indrv4 = indrv3 + *n;
 241     indrv5 = indrv4 + *n;
 242
 243 /*     Compute eigenvectors of matrix blocks. */
 244
 245     j1 = 1;
 246     i__1 = iblock[*m];
 247     for (nblk = 1; nblk <= i__1; ++nblk) {
 248
 249 /*        Find starting and ending indices of block nblk. */
 250
 251         if (nblk == 1) {
 252             b1 = 1;
 253         } else {
 254             b1 = isplit[nblk - 1] + 1;
 255         }
 256         bn = isplit[nblk];
 257         blksiz = bn - b1 + 1;
 258         if (blksiz == 1) {
 259             goto L60;
 260         }
 261         gpind = b1;
 262
 263 /*        Compute reorthogonalization criterion and stopping criterion. */
 264
 265         onenrm = (d__1 = d__[b1], abs(d__1)) + (d__2 = e[b1], abs(d__2));
 266 /* Computing MAX */
 267         d__3 = onenrm, d__4 = (d__1 = d__[bn], abs(d__1)) + (d__2 = e[bn - 1],
 268                  abs(d__2));
 269         onenrm = max(d__3,d__4);
 270         i__2 = bn - 1;
 271         for (i__ = b1 + 1; i__ <= i__2; ++i__) {
 272 /* Computing MAX */
 273             d__4 = onenrm, d__5 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[
 274                     i__ - 1], abs(d__2)) + (d__3 = e[i__], abs(d__3));
 275             onenrm = max(d__4,d__5);
 276 /* L50: */
 277         }
 278         ortol = onenrm * .001;
 279
 280         dtpcrt = sqrt(.1 / blksiz);
 281
 282 /*        Loop through eigenvalues of block nblk. */
 283
 284 L60:
 285         jblk = 0;
 286         i__2 = *m;
 287         for (j = j1; j <= i__2; ++j) {
 288             if (iblock[j] != nblk) {
 289                 j1 = j;
 290                 goto L160;
 291             }
 292             ++jblk;
 293             xj = w[j];
 294
 295 /*           Skip all the work if the block size is one. */
 296
 297             if (blksiz == 1) {
 298                 work[indrv1 + 1] = 1.;
 299                 goto L120;
 300             }
 301
 302 /*           If eigenvalues j and j-1 are too close, add a relatively */
 303 /*           small perturbation. */
 304
 305             if (jblk > 1) {
 306                 eps1 = (d__1 = eps * xj, abs(d__1));
 307                 pertol = eps1 * 10.;
 308                 sep = xj - xjm;
 309                 if (sep < pertol) {
 310                     xj = xjm + pertol;
 311                 }
 312             }
 313
 314             its = 0;
 315             nrmchk = 0;
 316
 317 /*           Get random starting vector. */
 318
 319             dlarnv_(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
 320
 321 /*           Copy the matrix T so it won't be destroyed in factorization. */
 322
 323             dcopy_(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
 324             i__3 = blksiz - 1;
 325             dcopy_(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
 326             i__3 = blksiz - 1;
 327             dcopy_(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
 328
 329 /*           Compute LU factors with partial pivoting  ( PT = LU ) */
 330
 331             tol = 0.;
 332             dlagtf_(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
 333                     indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
 334
 335 /*           Update iteration count. */
 336
 337 L70:
 338             ++its;
 339             if (its > 5) {
 340                 goto L100;
 341             }
 342
 343 /*           Normalize and scale the righthand side vector Pb. */
 344
 345 /* Computing MAX */
 346             d__2 = eps, d__3 = (d__1 = work[indrv4 + blksiz], abs(d__1));
 347             scl = blksiz * onenrm * max(d__2,d__3) / dasum_(&blksiz, &work[
 348                     indrv1 + 1], &c__1);
 349             dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
 350
 351 /*           Solve the system LU = Pb. */
 352
 353             dlagts_(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
 354                     work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
 355                     indrv1 + 1], &tol, &iinfo);
 356
 357 /*           Reorthogonalize by modified Gram-Schmidt if eigenvalues are */
 358 /*           close enough. */
 359
 360             if (jblk == 1) {
 361                 goto L90;
 362             }
 363             if ((d__1 = xj - xjm, abs(d__1)) > ortol) {
 364                 gpind = j;
 365             }
 366             if (gpind != j) {
 367                 i__3 = j - 1;
 368                 for (i__ = gpind; i__ <= i__3; ++i__) {
 369                     ztr = -ddot_(&blksiz, &work[indrv1 + 1], &c__1, &z__[b1 +
 370                             i__ * z_dim1], &c__1);
 371                     daxpy_(&blksiz, &ztr, &z__[b1 + i__ * z_dim1], &c__1, &
 372                             work[indrv1 + 1], &c__1);
 373 /* L80: */
 374                 }
 375             }
 376
 377 /*           Check the infinity norm of the iterate. */
 378
 379 L90:
 380             jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
 381             nrm = (d__1 = work[indrv1 + jmax], abs(d__1));
 382
 383 /*           Continue for additional iterations after norm reaches */
 384 /*           stopping criterion. */
 385
 386             if (nrm < dtpcrt) {
 387                 goto L70;
 388             }
 389             ++nrmchk;
 390             if (nrmchk < 3) {
 391                 goto L70;
 392             }
 393
 394             goto L110;
 395
 396 /*           If stopping criterion was not satisfied, update info and */
 397 /*           store eigenvector number in array ifail. */
 398
 399 L100:
 400             ++(*info);
 401             ifail[*info] = j;
 402
 403 /*           Accept iterate as jth eigenvector. */
 404
 405 L110:
 406             scl = 1. / dnrm2_(&blksiz, &work[indrv1 + 1], &c__1);
 407             jmax = idamax_(&blksiz, &work[indrv1 + 1], &c__1);
 408             if (work[indrv1 + jmax] < 0.) {
 409                 scl = -scl;
 410             }
 411             dscal_(&blksiz, &scl, &work[indrv1 + 1], &c__1);
 412 L120:
 413             i__3 = *n;
 414             for (i__ = 1; i__ <= i__3; ++i__) {
 415                 z__[i__ + j * z_dim1] = 0.;
 416 /* L130: */
 417             }
 418             i__3 = blksiz;
 419             for (i__ = 1; i__ <= i__3; ++i__) {
 420                 z__[b1 + i__ - 1 + j * z_dim1] = work[indrv1 + i__];
 421 /* L140: */
 422             }
 423
 424 /*           Save the shift to check eigenvalue spacing at next */
 425 /*           iteration. */
 426
 427             xjm = xj;
 428
 429 /* L150: */
 430         }
 431 L160:
 432         ;
 433     }
 434
 435     return 0;
 436
 437 /*     End of DSTEIN */
 438
 439 } /* dstein_ */