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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__9 = 9;
   6 static integer c__0 = 0;
   7 static integer c__2 = 2;
   8 static real c_b23 = 1.f;
   9 static real c_b24 = 0.f;
  10 static integer c__1 = 1;
  11
  12 /* Subroutine */ int slaed0_(integer *icompq, integer *qsiz, integer *n, real
  13         *d__, real *e, real *q, integer *ldq, real *qstore, integer *ldqs,
  14         real *work, integer *iwork, integer *info)
  15 {
  16     /* System generated locals */
  17     integer q_dim1, q_offset, qstore_dim1, qstore_offset, i__1, i__2;
  18     real r__1;
  19
  20     /* Builtin functions */
  21     double log(doublereal);
  22     integer pow_ii(integer *, integer *);
  23
  24     /* Local variables */
  25     integer i__, j, k, iq, lgn, msd2, smm1, spm1, spm2;
  26     real temp;
  27     integer curr;
  28     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
  29             integer *, real *, real *, integer *, real *, integer *, real *,
  30             real *, integer *);
  31     integer iperm, indxq, iwrem;
  32     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  33             integer *);
  34     integer iqptr, tlvls;
  35     extern /* Subroutine */ int slaed1_(integer *, real *, real *, integer *,
  36             integer *, real *, integer *, real *, integer *, integer *),
  37             slaed7_(integer *, integer *, integer *, integer *, integer *,
  38             integer *, real *, real *, integer *, integer *, real *, integer *
  39 , real *, integer *, integer *, integer *, integer *, integer *,
  40             real *, real *, integer *, integer *);
  41     integer igivcl;
  42     extern /* Subroutine */ int xerbla_(char *, integer *);
  43     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  44             integer *, integer *);
  45     integer igivnm, submat;
  46     extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *,
  47             integer *, real *, integer *);
  48     integer curprb, subpbs, igivpt, curlvl, matsiz, iprmpt, smlsiz;
  49     extern /* Subroutine */ int ssteqr_(char *, integer *, real *, real *,
  50             real *, integer *, real *, integer *);
  51
  52
  53 /*  -- LAPACK routine (version 3.1) -- */
  54 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  55 /*     November 2006 */
  56
  57 /*     .. Scalar Arguments .. */
  58 /*     .. */
  59 /*     .. Array Arguments .. */
  60 /*     .. */
  61
  62 /*  Purpose */
  63 /*  ======= */
  64
  65 /*  SLAED0 computes all eigenvalues and corresponding eigenvectors of a */
  66 /*  symmetric tridiagonal matrix using the divide and conquer method. */
  67
  68 /*  Arguments */
  69 /*  ========= */
  70
  71 /*  ICOMPQ  (input) INTEGER */
  72 /*          = 0:  Compute eigenvalues only. */
  73 /*          = 1:  Compute eigenvectors of original dense symmetric matrix */
  74 /*                also.  On entry, Q contains the orthogonal matrix used */
  75 /*                to reduce the original matrix to tridiagonal form. */
  76 /*          = 2:  Compute eigenvalues and eigenvectors of tridiagonal */
  77 /*                matrix. */
  78
  79 /*  QSIZ   (input) INTEGER */
  80 /*         The dimension of the orthogonal matrix used to reduce */
  81 /*         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1. */
  82
  83 /*  N      (input) INTEGER */
  84 /*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */
  85
  86 /*  D      (input/output) REAL array, dimension (N) */
  87 /*         On entry, the main diagonal of the tridiagonal matrix. */
  88 /*         On exit, its eigenvalues. */
  89
  90 /*  E      (input) REAL array, dimension (N-1) */
  91 /*         The off-diagonal elements of the tridiagonal matrix. */
  92 /*         On exit, E has been destroyed. */
  93
  94 /*  Q      (input/output) REAL array, dimension (LDQ, N) */
  95 /*         On entry, Q must contain an N-by-N orthogonal matrix. */
  96 /*         If ICOMPQ = 0    Q is not referenced. */
  97 /*         If ICOMPQ = 1    On entry, Q is a subset of the columns of the */
  98 /*                          orthogonal matrix used to reduce the full */
  99 /*                          matrix to tridiagonal form corresponding to */
 100 /*                          the subset of the full matrix which is being */
 101 /*                          decomposed at this time. */
 102 /*         If ICOMPQ = 2    On entry, Q will be the identity matrix. */
 103 /*                          On exit, Q contains the eigenvectors of the */
 104 /*                          tridiagonal matrix. */
 105
 106 /*  LDQ    (input) INTEGER */
 107 /*         The leading dimension of the array Q.  If eigenvectors are */
 108 /*         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1. */
 109
 110 /*  QSTORE (workspace) REAL array, dimension (LDQS, N) */
 111 /*         Referenced only when ICOMPQ = 1.  Used to store parts of */
 112 /*         the eigenvector matrix when the updating matrix multiplies */
 113 /*         take place. */
 114
 115 /*  LDQS   (input) INTEGER */
 116 /*         The leading dimension of the array QSTORE.  If ICOMPQ = 1, */
 117 /*         then  LDQS >= max(1,N).  In any case,  LDQS >= 1. */
 118
 119 /*  WORK   (workspace) REAL array, */
 120 /*         If ICOMPQ = 0 or 1, the dimension of WORK must be at least */
 121 /*                     1 + 3*N + 2*N*lg N + 2*N**2 */
 122 /*                     ( lg( N ) = smallest integer k */
 123 /*                                 such that 2^k >= N ) */
 124 /*         If ICOMPQ = 2, the dimension of WORK must be at least */
 125 /*                     4*N + N**2. */
 126
 127 /*  IWORK  (workspace) INTEGER array, */
 128 /*         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least */
 129 /*                        6 + 6*N + 5*N*lg N. */
 130 /*                        ( lg( N ) = smallest integer k */
 131 /*                                    such that 2^k >= N ) */
 132 /*         If ICOMPQ = 2, the dimension of IWORK must be at least */
 133 /*                        3 + 5*N. */
 134
 135 /*  INFO   (output) INTEGER */
 136 /*          = 0:  successful exit. */
 137 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 138 /*          > 0:  The algorithm failed to compute an eigenvalue while */
 139 /*                working on the submatrix lying in rows and columns */
 140 /*                INFO/(N+1) through mod(INFO,N+1). */
 141
 142 /*  Further Details */
 143 /*  =============== */
 144
 145 /*  Based on contributions by */
 146 /*     Jeff Rutter, Computer Science Division, University of California */
 147 /*     at Berkeley, USA */
 148
 149 /*  ===================================================================== */
 150
 151 /*     .. Parameters .. */
 152 /*     .. */
 153 /*     .. Local Scalars .. */
 154 /*     .. */
 155 /*     .. External Subroutines .. */
 156 /*     .. */
 157 /*     .. External Functions .. */
 158 /*     .. */
 159 /*     .. Intrinsic Functions .. */
 160 /*     .. */
 161 /*     .. Executable Statements .. */
 162
 163 /*     Test the input parameters. */
 164
 165     /* Parameter adjustments */
 166     --d__;
 167     --e;
 168     q_dim1 = *ldq;
 169     q_offset = 1 + q_dim1;
 170     q -= q_offset;
 171     qstore_dim1 = *ldqs;
 172     qstore_offset = 1 + qstore_dim1;
 173     qstore -= qstore_offset;
 174     --work;
 175     --iwork;
 176
 177     /* Function Body */
 178     *info = 0;
 179
 180     if (*icompq < 0 || *icompq > 2) {
 181         *info = -1;
 182     } else if (*icompq == 1 && *qsiz < max(0,*n)) {
 183         *info = -2;
 184     } else if (*n < 0) {
 185         *info = -3;
 186     } else if (*ldq < max(1,*n)) {
 187         *info = -7;
 188     } else if (*ldqs < max(1,*n)) {
 189         *info = -9;
 190     }
 191     if (*info != 0) {
 192         i__1 = -(*info);
 193         xerbla_("SLAED0", &i__1);
 194         return 0;
 195     }
 196
 197 /*     Quick return if possible */
 198
 199     if (*n == 0) {
 200         return 0;
 201     }
 202
 203     smlsiz = ilaenv_(&c__9, "SLAED0", " ", &c__0, &c__0, &c__0, &c__0);
 204
 205 /*     Determine the size and placement of the submatrices, and save in */
 206 /*     the leading elements of IWORK. */
 207
 208     iwork[1] = *n;
 209     subpbs = 1;
 210     tlvls = 0;
 211 L10:
 212     if (iwork[subpbs] > smlsiz) {
 213         for (j = subpbs; j >= 1; --j) {
 214             iwork[j * 2] = (iwork[j] + 1) / 2;
 215             iwork[(j << 1) - 1] = iwork[j] / 2;
 216 /* L20: */
 217         }
 218         ++tlvls;
 219         subpbs <<= 1;
 220         goto L10;
 221     }
 222     i__1 = subpbs;
 223     for (j = 2; j <= i__1; ++j) {
 224         iwork[j] += iwork[j - 1];
 225 /* L30: */
 226     }
 227
 228 /*     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 */
 229 /*     using rank-1 modifications (cuts). */
 230
 231     spm1 = subpbs - 1;
 232     i__1 = spm1;
 233     for (i__ = 1; i__ <= i__1; ++i__) {
 234         submat = iwork[i__] + 1;
 235         smm1 = submat - 1;
 236         d__[smm1] -= (r__1 = e[smm1], dabs(r__1));
 237         d__[submat] -= (r__1 = e[smm1], dabs(r__1));
 238 /* L40: */
 239     }
 240
 241     indxq = (*n << 2) + 3;
 242     if (*icompq != 2) {
 243
 244 /*        Set up workspaces for eigenvalues only/accumulate new vectors */
 245 /*        routine */
 246
 247         temp = log((real) (*n)) / log(2.f);
 248         lgn = (integer) temp;
 249         if (pow_ii(&c__2, &lgn) < *n) {
 250             ++lgn;
 251         }
 252         if (pow_ii(&c__2, &lgn) < *n) {
 253             ++lgn;
 254         }
 255         iprmpt = indxq + *n + 1;
 256         iperm = iprmpt + *n * lgn;
 257         iqptr = iperm + *n * lgn;
 258         igivpt = iqptr + *n + 2;
 259         igivcl = igivpt + *n * lgn;
 260
 261         igivnm = 1;
 262         iq = igivnm + (*n << 1) * lgn;
 263 /* Computing 2nd power */
 264         i__1 = *n;
 265         iwrem = iq + i__1 * i__1 + 1;
 266
 267 /*        Initialize pointers */
 268
 269         i__1 = subpbs;
 270         for (i__ = 0; i__ <= i__1; ++i__) {
 271             iwork[iprmpt + i__] = 1;
 272             iwork[igivpt + i__] = 1;
 273 /* L50: */
 274         }
 275         iwork[iqptr] = 1;
 276     }
 277
 278 /*     Solve each submatrix eigenproblem at the bottom of the divide and */
 279 /*     conquer tree. */
 280
 281     curr = 0;
 282     i__1 = spm1;
 283     for (i__ = 0; i__ <= i__1; ++i__) {
 284         if (i__ == 0) {
 285             submat = 1;
 286             matsiz = iwork[1];
 287         } else {
 288             submat = iwork[i__] + 1;
 289             matsiz = iwork[i__ + 1] - iwork[i__];
 290         }
 291         if (*icompq == 2) {
 292             ssteqr_("I", &matsiz, &d__[submat], &e[submat], &q[submat +
 293                     submat * q_dim1], ldq, &work[1], info);
 294             if (*info != 0) {
 295                 goto L130;
 296             }
 297         } else {
 298             ssteqr_("I", &matsiz, &d__[submat], &e[submat], &work[iq - 1 +
 299                     iwork[iqptr + curr]], &matsiz, &work[1], info);
 300             if (*info != 0) {
 301                 goto L130;
 302             }
 303             if (*icompq == 1) {
 304                 sgemm_("N", "N", qsiz, &matsiz, &matsiz, &c_b23, &q[submat *
 305                         q_dim1 + 1], ldq, &work[iq - 1 + iwork[iqptr + curr]],
 306                          &matsiz, &c_b24, &qstore[submat * qstore_dim1 + 1],
 307                         ldqs);
 308             }
 309 /* Computing 2nd power */
 310             i__2 = matsiz;
 311             iwork[iqptr + curr + 1] = iwork[iqptr + curr] + i__2 * i__2;
 312             ++curr;
 313         }
 314         k = 1;
 315         i__2 = iwork[i__ + 1];
 316         for (j = submat; j <= i__2; ++j) {
 317             iwork[indxq + j] = k;
 318             ++k;
 319 /* L60: */
 320         }
 321 /* L70: */
 322     }
 323
 324 /*     Successively merge eigensystems of adjacent submatrices */
 325 /*     into eigensystem for the corresponding larger matrix. */
 326
 327 /*     while ( SUBPBS > 1 ) */
 328
 329     curlvl = 1;
 330 L80:
 331     if (subpbs > 1) {
 332         spm2 = subpbs - 2;
 333         i__1 = spm2;
 334         for (i__ = 0; i__ <= i__1; i__ += 2) {
 335             if (i__ == 0) {
 336                 submat = 1;
 337                 matsiz = iwork[2];
 338                 msd2 = iwork[1];
 339                 curprb = 0;
 340             } else {
 341                 submat = iwork[i__] + 1;
 342                 matsiz = iwork[i__ + 2] - iwork[i__];
 343                 msd2 = matsiz / 2;
 344                 ++curprb;
 345             }
 346
 347 /*     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) */
 348 /*     into an eigensystem of size MATSIZ. */
 349 /*     SLAED1 is used only for the full eigensystem of a tridiagonal */
 350 /*     matrix. */
 351 /*     SLAED7 handles the cases in which eigenvalues only or eigenvalues */
 352 /*     and eigenvectors of a full symmetric matrix (which was reduced to */
 353 /*     tridiagonal form) are desired. */
 354
 355             if (*icompq == 2) {
 356                 slaed1_(&matsiz, &d__[submat], &q[submat + submat * q_dim1],
 357                         ldq, &iwork[indxq + submat], &e[submat + msd2 - 1], &
 358                         msd2, &work[1], &iwork[subpbs + 1], info);
 359             } else {
 360                 slaed7_(icompq, &matsiz, qsiz, &tlvls, &curlvl, &curprb, &d__[
 361                         submat], &qstore[submat * qstore_dim1 + 1], ldqs, &
 362                         iwork[indxq + submat], &e[submat + msd2 - 1], &msd2, &
 363                         work[iq], &iwork[iqptr], &iwork[iprmpt], &iwork[iperm]
 364 , &iwork[igivpt], &iwork[igivcl], &work[igivnm], &
 365                         work[iwrem], &iwork[subpbs + 1], info);
 366             }
 367             if (*info != 0) {
 368                 goto L130;
 369             }
 370             iwork[i__ / 2 + 1] = iwork[i__ + 2];
 371 /* L90: */
 372         }
 373         subpbs /= 2;
 374         ++curlvl;
 375         goto L80;
 376     }
 377
 378 /*     end while */
 379
 380 /*     Re-merge the eigenvalues/vectors which were deflated at the final */
 381 /*     merge step. */
 382
 383     if (*icompq == 1) {
 384         i__1 = *n;
 385         for (i__ = 1; i__ <= i__1; ++i__) {
 386             j = iwork[indxq + i__];
 387             work[i__] = d__[j];
 388             scopy_(qsiz, &qstore[j * qstore_dim1 + 1], &c__1, &q[i__ * q_dim1
 389                     + 1], &c__1);
 390 /* L100: */
 391         }
 392         scopy_(n, &work[1], &c__1, &d__[1], &c__1);
 393     } else if (*icompq == 2) {
 394         i__1 = *n;
 395         for (i__ = 1; i__ <= i__1; ++i__) {
 396             j = iwork[indxq + i__];
 397             work[i__] = d__[j];
 398             scopy_(n, &q[j * q_dim1 + 1], &c__1, &work[*n * i__ + 1], &c__1);
 399 /* L110: */
 400         }
 401         scopy_(n, &work[1], &c__1, &d__[1], &c__1);
 402         slacpy_("A", n, n, &work[*n + 1], n, &q[q_offset], ldq);
 403     } else {
 404         i__1 = *n;
 405         for (i__ = 1; i__ <= i__1; ++i__) {
 406             j = iwork[indxq + i__];
 407             work[i__] = d__[j];
 408 /* L120: */
 409         }
 410         scopy_(n, &work[1], &c__1, &d__[1], &c__1);
 411     }
 412     goto L140;
 413
 414 L130:
 415     *info = submat * (*n + 1) + submat + matsiz - 1;
 416
 417 L140:
 418     return 0;
 419
 420 /*     End of SLAED0 */
 421
 422 } /* slaed0_ */