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