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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static doublereal c_b22 = 1.;
   7 static doublereal c_b23 = 0.;
   8
   9 /* Subroutine */ int dlaed3_(integer *k, integer *n, integer *n1, doublereal *
  10         d__, doublereal *q, integer *ldq, doublereal *rho, doublereal *dlamda,
  11          doublereal *q2, integer *indx, integer *ctot, doublereal *w,
  12         doublereal *s, integer *info)
  13 {
  14     /* System generated locals */
  15     integer q_dim1, q_offset, i__1, i__2;
  16     doublereal d__1;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
  20
  21     /* Local variables */
  22     integer i__, j, n2, n12, ii, n23, iq2;
  23     doublereal temp;
  24     extern doublereal dnrm2_(integer *, doublereal *, integer *);
  25     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
  26             integer *, doublereal *, doublereal *, integer *, doublereal *,
  27             integer *, doublereal *, doublereal *, integer *),
  28              dcopy_(integer *, doublereal *, integer *, doublereal *, integer
  29             *), dlaed4_(integer *, integer *, doublereal *, doublereal *,
  30             doublereal *, doublereal *, doublereal *, integer *);
  31     extern doublereal dlamc3_(doublereal *, doublereal *);
  32     extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
  33             doublereal *, integer *, doublereal *, integer *),
  34             dlaset_(char *, integer *, integer *, doublereal *, doublereal *,
  35             doublereal *, integer *), xerbla_(char *, integer *);
  36
  37
  38 /*  -- LAPACK routine (version 3.1) -- */
  39 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  40 /*     November 2006 */
  41
  42 /*     .. Scalar Arguments .. */
  43 /*     .. */
  44 /*     .. Array Arguments .. */
  45 /*     .. */
  46
  47 /*  Purpose */
  48 /*  ======= */
  49
  50 /*  DLAED3 finds the roots of the secular equation, as defined by the */
  51 /*  values in D, W, and RHO, between 1 and K.  It makes the */
  52 /*  appropriate calls to DLAED4 and then updates the eigenvectors by */
  53 /*  multiplying the matrix of eigenvectors of the pair of eigensystems */
  54 /*  being combined by the matrix of eigenvectors of the K-by-K system */
  55 /*  which is solved here. */
  56
  57 /*  This code makes very mild assumptions about floating point */
  58 /*  arithmetic. It will work on machines with a guard digit in */
  59 /*  add/subtract, or on those binary machines without guard digits */
  60 /*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. */
  61 /*  It could conceivably fail on hexadecimal or decimal machines */
  62 /*  without guard digits, but we know of none. */
  63
  64 /*  Arguments */
  65 /*  ========= */
  66
  67 /*  K       (input) INTEGER */
  68 /*          The number of terms in the rational function to be solved by */
  69 /*          DLAED4.  K >= 0. */
  70
  71 /*  N       (input) INTEGER */
  72 /*          The number of rows and columns in the Q matrix. */
  73 /*          N >= K (deflation may result in N>K). */
  74
  75 /*  N1      (input) INTEGER */
  76 /*          The location of the last eigenvalue in the leading submatrix. */
  77 /*          min(1,N) <= N1 <= N/2. */
  78
  79 /*  D       (output) DOUBLE PRECISION array, dimension (N) */
  80 /*          D(I) contains the updated eigenvalues for */
  81 /*          1 <= I <= K. */
  82
  83 /*  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N) */
  84 /*          Initially the first K columns are used as workspace. */
  85 /*          On output the columns 1 to K contain */
  86 /*          the updated eigenvectors. */
  87
  88 /*  LDQ     (input) INTEGER */
  89 /*          The leading dimension of the array Q.  LDQ >= max(1,N). */
  90
  91 /*  RHO     (input) DOUBLE PRECISION */
  92 /*          The value of the parameter in the rank one update equation. */
  93 /*          RHO >= 0 required. */
  94
  95 /*  DLAMDA  (input/output) DOUBLE PRECISION array, dimension (K) */
  96 /*          The first K elements of this array contain the old roots */
  97 /*          of the deflated updating problem.  These are the poles */
  98 /*          of the secular equation. May be changed on output by */
  99 /*          having lowest order bit set to zero on Cray X-MP, Cray Y-MP, */
 100 /*          Cray-2, or Cray C-90, as described above. */
 101
 102 /*  Q2      (input) DOUBLE PRECISION array, dimension (LDQ2, N) */
 103 /*          The first K columns of this matrix contain the non-deflated */
 104 /*          eigenvectors for the split problem. */
 105
 106 /*  INDX    (input) INTEGER array, dimension (N) */
 107 /*          The permutation used to arrange the columns of the deflated */
 108 /*          Q matrix into three groups (see DLAED2). */
 109 /*          The rows of the eigenvectors found by DLAED4 must be likewise */
 110 /*          permuted before the matrix multiply can take place. */
 111
 112 /*  CTOT    (input) INTEGER array, dimension (4) */
 113 /*          A count of the total number of the various types of columns */
 114 /*          in Q, as described in INDX.  The fourth column type is any */
 115 /*          column which has been deflated. */
 116
 117 /*  W       (input/output) DOUBLE PRECISION array, dimension (K) */
 118 /*          The first K elements of this array contain the components */
 119 /*          of the deflation-adjusted updating vector. Destroyed on */
 120 /*          output. */
 121
 122 /*  S       (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K */
 123 /*          Will contain the eigenvectors of the repaired matrix which */
 124 /*          will be multiplied by the previously accumulated eigenvectors */
 125 /*          to update the system. */
 126
 127 /*  LDS     (input) INTEGER */
 128 /*          The leading dimension of S.  LDS >= max(1,K). */
 129
 130 /*  INFO    (output) INTEGER */
 131 /*          = 0:  successful exit. */
 132 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 133 /*          > 0:  if INFO = 1, an eigenvalue did not converge */
 134
 135 /*  Further Details */
 136 /*  =============== */
 137
 138 /*  Based on contributions by */
 139 /*     Jeff Rutter, Computer Science Division, University of California */
 140 /*     at Berkeley, USA */
 141 /*  Modified by Francoise Tisseur, University of Tennessee. */
 142
 143 /*  ===================================================================== */
 144
 145 /*     .. Parameters .. */
 146 /*     .. */
 147 /*     .. Local Scalars .. */
 148 /*     .. */
 149 /*     .. External Functions .. */
 150 /*     .. */
 151 /*     .. External Subroutines .. */
 152 /*     .. */
 153 /*     .. Intrinsic Functions .. */
 154 /*     .. */
 155 /*     .. Executable Statements .. */
 156
 157 /*     Test the input parameters. */
 158
 159     /* Parameter adjustments */
 160     --d__;
 161     q_dim1 = *ldq;
 162     q_offset = 1 + q_dim1;
 163     q -= q_offset;
 164     --dlamda;
 165     --q2;
 166     --indx;
 167     --ctot;
 168     --w;
 169     --s;
 170
 171     /* Function Body */
 172     *info = 0;
 173
 174     if (*k < 0) {
 175         *info = -1;
 176     } else if (*n < *k) {
 177         *info = -2;
 178     } else if (*ldq < max(1,*n)) {
 179         *info = -6;
 180     }
 181     if (*info != 0) {
 182         i__1 = -(*info);
 183         xerbla_("DLAED3", &i__1);
 184         return 0;
 185     }
 186
 187 /*     Quick return if possible */
 188
 189     if (*k == 0) {
 190         return 0;
 191     }
 192
 193 /*     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can */
 194 /*     be computed with high relative accuracy (barring over/underflow). */
 195 /*     This is a problem on machines without a guard digit in */
 196 /*     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). */
 197 /*     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), */
 198 /*     which on any of these machines zeros out the bottommost */
 199 /*     bit of DLAMDA(I) if it is 1; this makes the subsequent */
 200 /*     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation */
 201 /*     occurs. On binary machines with a guard digit (almost all */
 202 /*     machines) it does not change DLAMDA(I) at all. On hexadecimal */
 203 /*     and decimal machines with a guard digit, it slightly */
 204 /*     changes the bottommost bits of DLAMDA(I). It does not account */
 205 /*     for hexadecimal or decimal machines without guard digits */
 206 /*     (we know of none). We use a subroutine call to compute */
 207 /*     2*DLAMBDA(I) to prevent optimizing compilers from eliminating */
 208 /*     this code. */
 209
 210     i__1 = *k;
 211     for (i__ = 1; i__ <= i__1; ++i__) {
 212         dlamda[i__] = dlamc3_(&dlamda[i__], &dlamda[i__]) - dlamda[i__];
 213 /* L10: */
 214     }
 215
 216     i__1 = *k;
 217     for (j = 1; j <= i__1; ++j) {
 218         dlaed4_(k, &j, &dlamda[1], &w[1], &q[j * q_dim1 + 1], rho, &d__[j],
 219                 info);
 220
 221 /*        If the zero finder fails, the computation is terminated. */
 222
 223         if (*info != 0) {
 224             goto L120;
 225         }
 226 /* L20: */
 227     }
 228
 229     if (*k == 1) {
 230         goto L110;
 231     }
 232     if (*k == 2) {
 233         i__1 = *k;
 234         for (j = 1; j <= i__1; ++j) {
 235             w[1] = q[j * q_dim1 + 1];
 236             w[2] = q[j * q_dim1 + 2];
 237             ii = indx[1];
 238             q[j * q_dim1 + 1] = w[ii];
 239             ii = indx[2];
 240             q[j * q_dim1 + 2] = w[ii];
 241 /* L30: */
 242         }
 243         goto L110;
 244     }
 245
 246 /*     Compute updated W. */
 247
 248     dcopy_(k, &w[1], &c__1, &s[1], &c__1);
 249
 250 /*     Initialize W(I) = Q(I,I) */
 251
 252     i__1 = *ldq + 1;
 253     dcopy_(k, &q[q_offset], &i__1, &w[1], &c__1);
 254     i__1 = *k;
 255     for (j = 1; j <= i__1; ++j) {
 256         i__2 = j - 1;
 257         for (i__ = 1; i__ <= i__2; ++i__) {
 258             w[i__] *= q[i__ + j * q_dim1] / (dlamda[i__] - dlamda[j]);
 259 /* L40: */
 260         }
 261         i__2 = *k;
 262         for (i__ = j + 1; i__ <= i__2; ++i__) {
 263             w[i__] *= q[i__ + j * q_dim1] / (dlamda[i__] - dlamda[j]);
 264 /* L50: */
 265         }
 266 /* L60: */
 267     }
 268     i__1 = *k;
 269     for (i__ = 1; i__ <= i__1; ++i__) {
 270         d__1 = sqrt(-w[i__]);
 271         w[i__] = d_sign(&d__1, &s[i__]);
 272 /* L70: */
 273     }
 274
 275 /*     Compute eigenvectors of the modified rank-1 modification. */
 276
 277     i__1 = *k;
 278     for (j = 1; j <= i__1; ++j) {
 279         i__2 = *k;
 280         for (i__ = 1; i__ <= i__2; ++i__) {
 281             s[i__] = w[i__] / q[i__ + j * q_dim1];
 282 /* L80: */
 283         }
 284         temp = dnrm2_(k, &s[1], &c__1);
 285         i__2 = *k;
 286         for (i__ = 1; i__ <= i__2; ++i__) {
 287             ii = indx[i__];
 288             q[i__ + j * q_dim1] = s[ii] / temp;
 289 /* L90: */
 290         }
 291 /* L100: */
 292     }
 293
 294 /*     Compute the updated eigenvectors. */
 295
 296 L110:
 297
 298     n2 = *n - *n1;
 299     n12 = ctot[1] + ctot[2];
 300     n23 = ctot[2] + ctot[3];
 301
 302     dlacpy_("A", &n23, k, &q[ctot[1] + 1 + q_dim1], ldq, &s[1], &n23);
 303     iq2 = *n1 * n12 + 1;
 304     if (n23 != 0) {
 305         dgemm_("N", "N", &n2, k, &n23, &c_b22, &q2[iq2], &n2, &s[1], &n23, &
 306                 c_b23, &q[*n1 + 1 + q_dim1], ldq);
 307     } else {
 308         dlaset_("A", &n2, k, &c_b23, &c_b23, &q[*n1 + 1 + q_dim1], ldq);
 309     }
 310
 311     dlacpy_("A", &n12, k, &q[q_offset], ldq, &s[1], &n12);
 312     if (n12 != 0) {
 313         dgemm_("N", "N", n1, k, &n12, &c_b22, &q2[1], n1, &s[1], &n12, &c_b23,
 314                  &q[q_offset], ldq);
 315     } else {
 316         dlaset_("A", n1, k, &c_b23, &c_b23, &q[q_dim1 + 1], ldq);
 317     }
 318
 319
 320 L120:
 321     return 0;
 322
 323 /*     End of DLAED3 */
 324
 325 } /* dlaed3_ */