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