slaed8.c
Go to the documentation of this file.
00001 /* slaed8.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static real c_b3 = -1.f;
00019 static integer c__1 = 1;
00020 
00021 /* Subroutine */ int slaed8_(integer *icompq, integer *k, integer *n, integer 
00022         *qsiz, real *d__, real *q, integer *ldq, integer *indxq, real *rho, 
00023         integer *cutpnt, real *z__, real *dlamda, real *q2, integer *ldq2, 
00024         real *w, integer *perm, integer *givptr, integer *givcol, real *
00025         givnum, integer *indxp, integer *indx, integer *info)
00026 {
00027     /* System generated locals */
00028     integer q_dim1, q_offset, q2_dim1, q2_offset, i__1;
00029     real r__1;
00030 
00031     /* Builtin functions */
00032     double sqrt(doublereal);
00033 
00034     /* Local variables */
00035     real c__;
00036     integer i__, j;
00037     real s, t;
00038     integer k2, n1, n2, jp, n1p1;
00039     real eps, tau, tol;
00040     integer jlam, imax, jmax;
00041     extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 
00042             integer *, real *, real *), sscal_(integer *, real *, real *, 
00043             integer *), scopy_(integer *, real *, integer *, real *, integer *
00044 );
00045     extern doublereal slapy2_(real *, real *), slamch_(char *);
00046     extern /* Subroutine */ int xerbla_(char *, integer *);
00047     extern integer isamax_(integer *, real *, integer *);
00048     extern /* Subroutine */ int slamrg_(integer *, integer *, real *, integer 
00049             *, integer *, integer *), slacpy_(char *, integer *, integer *, 
00050             real *, integer *, real *, integer *);
00051 
00052 
00053 /*  -- LAPACK routine (version 3.2) -- */
00054 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00055 /*     November 2006 */
00056 
00057 /*     .. Scalar Arguments .. */
00058 /*     .. */
00059 /*     .. Array Arguments .. */
00060 /*     .. */
00061 
00062 /*  Purpose */
00063 /*  ======= */
00064 
00065 /*  SLAED8 merges the two sets of eigenvalues together into a single */
00066 /*  sorted set.  Then it tries to deflate the size of the problem. */
00067 /*  There are two ways in which deflation can occur:  when two or more */
00068 /*  eigenvalues are close together or if there is a tiny element in the */
00069 /*  Z vector.  For each such occurrence the order of the related secular */
00070 /*  equation problem is reduced by one. */
00071 
00072 /*  Arguments */
00073 /*  ========= */
00074 
00075 /*  ICOMPQ  (input) INTEGER */
00076 /*          = 0:  Compute eigenvalues only. */
00077 /*          = 1:  Compute eigenvectors of original dense symmetric matrix */
00078 /*                also.  On entry, Q contains the orthogonal matrix used */
00079 /*                to reduce the original matrix to tridiagonal form. */
00080 
00081 /*  K      (output) INTEGER */
00082 /*         The number of non-deflated eigenvalues, and the order of the */
00083 /*         related secular equation. */
00084 
00085 /*  N      (input) INTEGER */
00086 /*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */
00087 
00088 /*  QSIZ   (input) INTEGER */
00089 /*         The dimension of the orthogonal matrix used to reduce */
00090 /*         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1. */
00091 
00092 /*  D      (input/output) REAL array, dimension (N) */
00093 /*         On entry, the eigenvalues of the two submatrices to be */
00094 /*         combined.  On exit, the trailing (N-K) updated eigenvalues */
00095 /*         (those which were deflated) sorted into increasing order. */
00096 
00097 /*  Q      (input/output) REAL array, dimension (LDQ,N) */
00098 /*         If ICOMPQ = 0, Q is not referenced.  Otherwise, */
00099 /*         on entry, Q contains the eigenvectors of the partially solved */
00100 /*         system which has been previously updated in matrix */
00101 /*         multiplies with other partially solved eigensystems. */
00102 /*         On exit, Q contains the trailing (N-K) updated eigenvectors */
00103 /*         (those which were deflated) in its last N-K columns. */
00104 
00105 /*  LDQ    (input) INTEGER */
00106 /*         The leading dimension of the array Q.  LDQ >= max(1,N). */
00107 
00108 /*  INDXQ  (input) INTEGER array, dimension (N) */
00109 /*         The permutation which separately sorts the two sub-problems */
00110 /*         in D into ascending order.  Note that elements in the second */
00111 /*         half of this permutation must first have CUTPNT added to */
00112 /*         their values in order to be accurate. */
00113 
00114 /*  RHO    (input/output) REAL */
00115 /*         On entry, the off-diagonal element associated with the rank-1 */
00116 /*         cut which originally split the two submatrices which are now */
00117 /*         being recombined. */
00118 /*         On exit, RHO has been modified to the value required by */
00119 /*         SLAED3. */
00120 
00121 /*  CUTPNT (input) INTEGER */
00122 /*         The location of the last eigenvalue in the leading */
00123 /*         sub-matrix.  min(1,N) <= CUTPNT <= N. */
00124 
00125 /*  Z      (input) REAL array, dimension (N) */
00126 /*         On entry, Z contains the updating vector (the last row of */
00127 /*         the first sub-eigenvector matrix and the first row of the */
00128 /*         second sub-eigenvector matrix). */
00129 /*         On exit, the contents of Z are destroyed by the updating */
00130 /*         process. */
00131 
00132 /*  DLAMDA (output) REAL array, dimension (N) */
00133 /*         A copy of the first K eigenvalues which will be used by */
00134 /*         SLAED3 to form the secular equation. */
00135 
00136 /*  Q2     (output) REAL array, dimension (LDQ2,N) */
00137 /*         If ICOMPQ = 0, Q2 is not referenced.  Otherwise, */
00138 /*         a copy of the first K eigenvectors which will be used by */
00139 /*         SLAED7 in a matrix multiply (SGEMM) to update the new */
00140 /*         eigenvectors. */
00141 
00142 /*  LDQ2   (input) INTEGER */
00143 /*         The leading dimension of the array Q2.  LDQ2 >= max(1,N). */
00144 
00145 /*  W      (output) REAL array, dimension (N) */
00146 /*         The first k values of the final deflation-altered z-vector and */
00147 /*         will be passed to SLAED3. */
00148 
00149 /*  PERM   (output) INTEGER array, dimension (N) */
00150 /*         The permutations (from deflation and sorting) to be applied */
00151 /*         to each eigenblock. */
00152 
00153 /*  GIVPTR (output) INTEGER */
00154 /*         The number of Givens rotations which took place in this */
00155 /*         subproblem. */
00156 
00157 /*  GIVCOL (output) INTEGER array, dimension (2, N) */
00158 /*         Each pair of numbers indicates a pair of columns to take place */
00159 /*         in a Givens rotation. */
00160 
00161 /*  GIVNUM (output) REAL array, dimension (2, N) */
00162 /*         Each number indicates the S value to be used in the */
00163 /*         corresponding Givens rotation. */
00164 
00165 /*  INDXP  (workspace) INTEGER array, dimension (N) */
00166 /*         The permutation used to place deflated values of D at the end */
00167 /*         of the array.  INDXP(1:K) points to the nondeflated D-values */
00168 /*         and INDXP(K+1:N) points to the deflated eigenvalues. */
00169 
00170 /*  INDX   (workspace) INTEGER array, dimension (N) */
00171 /*         The permutation used to sort the contents of D into ascending */
00172 /*         order. */
00173 
00174 /*  INFO   (output) INTEGER */
00175 /*          = 0:  successful exit. */
00176 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00177 
00178 /*  Further Details */
00179 /*  =============== */
00180 
00181 /*  Based on contributions by */
00182 /*     Jeff Rutter, Computer Science Division, University of California */
00183 /*     at Berkeley, USA */
00184 
00185 /*  ===================================================================== */
00186 
00187 /*     .. Parameters .. */
00188 /*     .. */
00189 /*     .. Local Scalars .. */
00190 
00191 /*     .. */
00192 /*     .. External Functions .. */
00193 /*     .. */
00194 /*     .. External Subroutines .. */
00195 /*     .. */
00196 /*     .. Intrinsic Functions .. */
00197 /*     .. */
00198 /*     .. Executable Statements .. */
00199 
00200 /*     Test the input parameters. */
00201 
00202     /* Parameter adjustments */
00203     --d__;
00204     q_dim1 = *ldq;
00205     q_offset = 1 + q_dim1;
00206     q -= q_offset;
00207     --indxq;
00208     --z__;
00209     --dlamda;
00210     q2_dim1 = *ldq2;
00211     q2_offset = 1 + q2_dim1;
00212     q2 -= q2_offset;
00213     --w;
00214     --perm;
00215     givcol -= 3;
00216     givnum -= 3;
00217     --indxp;
00218     --indx;
00219 
00220     /* Function Body */
00221     *info = 0;
00222 
00223     if (*icompq < 0 || *icompq > 1) {
00224         *info = -1;
00225     } else if (*n < 0) {
00226         *info = -3;
00227     } else if (*icompq == 1 && *qsiz < *n) {
00228         *info = -4;
00229     } else if (*ldq < max(1,*n)) {
00230         *info = -7;
00231     } else if (*cutpnt < min(1,*n) || *cutpnt > *n) {
00232         *info = -10;
00233     } else if (*ldq2 < max(1,*n)) {
00234         *info = -14;
00235     }
00236     if (*info != 0) {
00237         i__1 = -(*info);
00238         xerbla_("SLAED8", &i__1);
00239         return 0;
00240     }
00241 
00242 /*     Quick return if possible */
00243 
00244     if (*n == 0) {
00245         return 0;
00246     }
00247 
00248     n1 = *cutpnt;
00249     n2 = *n - n1;
00250     n1p1 = n1 + 1;
00251 
00252     if (*rho < 0.f) {
00253         sscal_(&n2, &c_b3, &z__[n1p1], &c__1);
00254     }
00255 
00256 /*     Normalize z so that norm(z) = 1 */
00257 
00258     t = 1.f / sqrt(2.f);
00259     i__1 = *n;
00260     for (j = 1; j <= i__1; ++j) {
00261         indx[j] = j;
00262 /* L10: */
00263     }
00264     sscal_(n, &t, &z__[1], &c__1);
00265     *rho = (r__1 = *rho * 2.f, dabs(r__1));
00266 
00267 /*     Sort the eigenvalues into increasing order */
00268 
00269     i__1 = *n;
00270     for (i__ = *cutpnt + 1; i__ <= i__1; ++i__) {
00271         indxq[i__] += *cutpnt;
00272 /* L20: */
00273     }
00274     i__1 = *n;
00275     for (i__ = 1; i__ <= i__1; ++i__) {
00276         dlamda[i__] = d__[indxq[i__]];
00277         w[i__] = z__[indxq[i__]];
00278 /* L30: */
00279     }
00280     i__ = 1;
00281     j = *cutpnt + 1;
00282     slamrg_(&n1, &n2, &dlamda[1], &c__1, &c__1, &indx[1]);
00283     i__1 = *n;
00284     for (i__ = 1; i__ <= i__1; ++i__) {
00285         d__[i__] = dlamda[indx[i__]];
00286         z__[i__] = w[indx[i__]];
00287 /* L40: */
00288     }
00289 
00290 /*     Calculate the allowable deflation tolerence */
00291 
00292     imax = isamax_(n, &z__[1], &c__1);
00293     jmax = isamax_(n, &d__[1], &c__1);
00294     eps = slamch_("Epsilon");
00295     tol = eps * 8.f * (r__1 = d__[jmax], dabs(r__1));
00296 
00297 /*     If the rank-1 modifier is small enough, no more needs to be done */
00298 /*     except to reorganize Q so that its columns correspond with the */
00299 /*     elements in D. */
00300 
00301     if (*rho * (r__1 = z__[imax], dabs(r__1)) <= tol) {
00302         *k = 0;
00303         if (*icompq == 0) {
00304             i__1 = *n;
00305             for (j = 1; j <= i__1; ++j) {
00306                 perm[j] = indxq[indx[j]];
00307 /* L50: */
00308             }
00309         } else {
00310             i__1 = *n;
00311             for (j = 1; j <= i__1; ++j) {
00312                 perm[j] = indxq[indx[j]];
00313                 scopy_(qsiz, &q[perm[j] * q_dim1 + 1], &c__1, &q2[j * q2_dim1 
00314                         + 1], &c__1);
00315 /* L60: */
00316             }
00317             slacpy_("A", qsiz, n, &q2[q2_dim1 + 1], ldq2, &q[q_dim1 + 1], ldq);
00318         }
00319         return 0;
00320     }
00321 
00322 /*     If there are multiple eigenvalues then the problem deflates.  Here */
00323 /*     the number of equal eigenvalues are found.  As each equal */
00324 /*     eigenvalue is found, an elementary reflector is computed to rotate */
00325 /*     the corresponding eigensubspace so that the corresponding */
00326 /*     components of Z are zero in this new basis. */
00327 
00328     *k = 0;
00329     *givptr = 0;
00330     k2 = *n + 1;
00331     i__1 = *n;
00332     for (j = 1; j <= i__1; ++j) {
00333         if (*rho * (r__1 = z__[j], dabs(r__1)) <= tol) {
00334 
00335 /*           Deflate due to small z component. */
00336 
00337             --k2;
00338             indxp[k2] = j;
00339             if (j == *n) {
00340                 goto L110;
00341             }
00342         } else {
00343             jlam = j;
00344             goto L80;
00345         }
00346 /* L70: */
00347     }
00348 L80:
00349     ++j;
00350     if (j > *n) {
00351         goto L100;
00352     }
00353     if (*rho * (r__1 = z__[j], dabs(r__1)) <= tol) {
00354 
00355 /*        Deflate due to small z component. */
00356 
00357         --k2;
00358         indxp[k2] = j;
00359     } else {
00360 
00361 /*        Check if eigenvalues are close enough to allow deflation. */
00362 
00363         s = z__[jlam];
00364         c__ = z__[j];
00365 
00366 /*        Find sqrt(a**2+b**2) without overflow or */
00367 /*        destructive underflow. */
00368 
00369         tau = slapy2_(&c__, &s);
00370         t = d__[j] - d__[jlam];
00371         c__ /= tau;
00372         s = -s / tau;
00373         if ((r__1 = t * c__ * s, dabs(r__1)) <= tol) {
00374 
00375 /*           Deflation is possible. */
00376 
00377             z__[j] = tau;
00378             z__[jlam] = 0.f;
00379 
00380 /*           Record the appropriate Givens rotation */
00381 
00382             ++(*givptr);
00383             givcol[(*givptr << 1) + 1] = indxq[indx[jlam]];
00384             givcol[(*givptr << 1) + 2] = indxq[indx[j]];
00385             givnum[(*givptr << 1) + 1] = c__;
00386             givnum[(*givptr << 1) + 2] = s;
00387             if (*icompq == 1) {
00388                 srot_(qsiz, &q[indxq[indx[jlam]] * q_dim1 + 1], &c__1, &q[
00389                         indxq[indx[j]] * q_dim1 + 1], &c__1, &c__, &s);
00390             }
00391             t = d__[jlam] * c__ * c__ + d__[j] * s * s;
00392             d__[j] = d__[jlam] * s * s + d__[j] * c__ * c__;
00393             d__[jlam] = t;
00394             --k2;
00395             i__ = 1;
00396 L90:
00397             if (k2 + i__ <= *n) {
00398                 if (d__[jlam] < d__[indxp[k2 + i__]]) {
00399                     indxp[k2 + i__ - 1] = indxp[k2 + i__];
00400                     indxp[k2 + i__] = jlam;
00401                     ++i__;
00402                     goto L90;
00403                 } else {
00404                     indxp[k2 + i__ - 1] = jlam;
00405                 }
00406             } else {
00407                 indxp[k2 + i__ - 1] = jlam;
00408             }
00409             jlam = j;
00410         } else {
00411             ++(*k);
00412             w[*k] = z__[jlam];
00413             dlamda[*k] = d__[jlam];
00414             indxp[*k] = jlam;
00415             jlam = j;
00416         }
00417     }
00418     goto L80;
00419 L100:
00420 
00421 /*     Record the last eigenvalue. */
00422 
00423     ++(*k);
00424     w[*k] = z__[jlam];
00425     dlamda[*k] = d__[jlam];
00426     indxp[*k] = jlam;
00427 
00428 L110:
00429 
00430 /*     Sort the eigenvalues and corresponding eigenvectors into DLAMDA */
00431 /*     and Q2 respectively.  The eigenvalues/vectors which were not */
00432 /*     deflated go into the first K slots of DLAMDA and Q2 respectively, */
00433 /*     while those which were deflated go into the last N - K slots. */
00434 
00435     if (*icompq == 0) {
00436         i__1 = *n;
00437         for (j = 1; j <= i__1; ++j) {
00438             jp = indxp[j];
00439             dlamda[j] = d__[jp];
00440             perm[j] = indxq[indx[jp]];
00441 /* L120: */
00442         }
00443     } else {
00444         i__1 = *n;
00445         for (j = 1; j <= i__1; ++j) {
00446             jp = indxp[j];
00447             dlamda[j] = d__[jp];
00448             perm[j] = indxq[indx[jp]];
00449             scopy_(qsiz, &q[perm[j] * q_dim1 + 1], &c__1, &q2[j * q2_dim1 + 1]
00450 , &c__1);
00451 /* L130: */
00452         }
00453     }
00454 
00455 /*     The deflated eigenvalues and their corresponding vectors go back */
00456 /*     into the last N - K slots of D and Q respectively. */
00457 
00458     if (*k < *n) {
00459         if (*icompq == 0) {
00460             i__1 = *n - *k;
00461             scopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
00462         } else {
00463             i__1 = *n - *k;
00464             scopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
00465             i__1 = *n - *k;
00466             slacpy_("A", qsiz, &i__1, &q2[(*k + 1) * q2_dim1 + 1], ldq2, &q[(*
00467                     k + 1) * q_dim1 + 1], ldq);
00468         }
00469     }
00470 
00471     return 0;
00472 
00473 /*     End of SLAED8 */
00474 
00475 } /* slaed8_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:10