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


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