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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c_n1 = -1;
   7
   8 /* Subroutine */ int dlaed1_(integer *n, doublereal *d__, doublereal *q,
   9         integer *ldq, integer *indxq, doublereal *rho, integer *cutpnt,
  10         doublereal *work, integer *iwork, integer *info)
  11 {
  12     /* System generated locals */
  13     integer q_dim1, q_offset, i__1, i__2;
  14
  15     /* Local variables */
  16     integer i__, k, n1, n2, is, iw, iz, iq2, zpp1, indx, indxc;
  17     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  18             doublereal *, integer *);
  19     integer indxp;
  20     extern /* Subroutine */ int dlaed2_(integer *, integer *, integer *,
  21             doublereal *, doublereal *, integer *, integer *, doublereal *,
  22             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
  23              integer *, integer *, integer *, integer *), dlaed3_(integer *,
  24             integer *, integer *, doublereal *, doublereal *, integer *,
  25             doublereal *, doublereal *, doublereal *, integer *, integer *,
  26             doublereal *, doublereal *, integer *);
  27     integer idlmda;
  28     extern /* Subroutine */ int dlamrg_(integer *, integer *, doublereal *,
  29             integer *, integer *, integer *), xerbla_(char *, integer *);
  30     integer coltyp;
  31
  32
  33 /*  -- LAPACK routine (version 3.1) -- */
  34 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  35 /*     November 2006 */
  36
  37 /*     .. Scalar Arguments .. */
  38 /*     .. */
  39 /*     .. Array Arguments .. */
  40 /*     .. */
  41
  42 /*  Purpose */
  43 /*  ======= */
  44
  45 /*  DLAED1 computes the updated eigensystem of a diagonal */
  46 /*  matrix after modification by a rank-one symmetric matrix.  This */
  47 /*  routine is used only for the eigenproblem which requires all */
  48 /*  eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles */
  49 /*  the case in which eigenvalues only or eigenvalues and eigenvectors */
  50 /*  of a full symmetric matrix (which was reduced to tridiagonal form) */
  51 /*  are desired. */
  52
  53 /*    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) */
  54
  55 /*     where Z = Q'u, u is a vector of length N with ones in the */
  56 /*     CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. */
  57
  58 /*     The eigenvectors of the original matrix are stored in Q, and the */
  59 /*     eigenvalues are in D.  The algorithm consists of three stages: */
  60
  61 /*        The first stage consists of deflating the size of the problem */
  62 /*        when there are multiple eigenvalues or if there is a zero in */
  63 /*        the Z vector.  For each such occurence the dimension of the */
  64 /*        secular equation problem is reduced by one.  This stage is */
  65 /*        performed by the routine DLAED2. */
  66
  67 /*        The second stage consists of calculating the updated */
  68 /*        eigenvalues. This is done by finding the roots of the secular */
  69 /*        equation via the routine DLAED4 (as called by DLAED3). */
  70 /*        This routine also calculates the eigenvectors of the current */
  71 /*        problem. */
  72
  73 /*        The final stage consists of computing the updated eigenvectors */
  74 /*        directly using the updated eigenvalues.  The eigenvectors for */
  75 /*        the current problem are multiplied with the eigenvectors from */
  76 /*        the overall problem. */
  77
  78 /*  Arguments */
  79 /*  ========= */
  80
  81 /*  N      (input) INTEGER */
  82 /*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */
  83
  84 /*  D      (input/output) DOUBLE PRECISION array, dimension (N) */
  85 /*         On entry, the eigenvalues of the rank-1-perturbed matrix. */
  86 /*         On exit, the eigenvalues of the repaired matrix. */
  87
  88 /*  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ,N) */
  89 /*         On entry, the eigenvectors of the rank-1-perturbed matrix. */
  90 /*         On exit, the eigenvectors of the repaired tridiagonal matrix. */
  91
  92 /*  LDQ    (input) INTEGER */
  93 /*         The leading dimension of the array Q.  LDQ >= max(1,N). */
  94
  95 /*  INDXQ  (input/output) INTEGER array, dimension (N) */
  96 /*         On entry, the permutation which separately sorts the two */
  97 /*         subproblems in D into ascending order. */
  98 /*         On exit, the permutation which will reintegrate the */
  99 /*         subproblems back into sorted order, */
 100 /*         i.e. D( INDXQ( I = 1, N ) ) will be in ascending order. */
 101
 102 /*  RHO    (input) DOUBLE PRECISION */
 103 /*         The subdiagonal entry used to create the rank-1 modification. */
 104
 105 /*  CUTPNT (input) INTEGER */
 106 /*         The location of the last eigenvalue in the leading sub-matrix. */
 107 /*         min(1,N) <= CUTPNT <= N/2. */
 108
 109 /*  WORK   (workspace) DOUBLE PRECISION array, dimension (4*N + N**2) */
 110
 111 /*  IWORK  (workspace) INTEGER array, dimension (4*N) */
 112
 113 /*  INFO   (output) INTEGER */
 114 /*          = 0:  successful exit. */
 115 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 116 /*          > 0:  if INFO = 1, an eigenvalue did not converge */
 117
 118 /*  Further Details */
 119 /*  =============== */
 120
 121 /*  Based on contributions by */
 122 /*     Jeff Rutter, Computer Science Division, University of California */
 123 /*     at Berkeley, USA */
 124 /*  Modified by Francoise Tisseur, University of Tennessee. */
 125
 126 /*  ===================================================================== */
 127
 128 /*     .. Local Scalars .. */
 129 /*     .. */
 130 /*     .. External Subroutines .. */
 131 /*     .. */
 132 /*     .. Intrinsic Functions .. */
 133 /*     .. */
 134 /*     .. Executable Statements .. */
 135
 136 /*     Test the input parameters. */
 137
 138     /* Parameter adjustments */
 139     --d__;
 140     q_dim1 = *ldq;
 141     q_offset = 1 + q_dim1;
 142     q -= q_offset;
 143     --indxq;
 144     --work;
 145     --iwork;
 146
 147     /* Function Body */
 148     *info = 0;
 149
 150     if (*n < 0) {
 151         *info = -1;
 152     } else if (*ldq < max(1,*n)) {
 153         *info = -4;
 154     } else /* if(complicated condition) */ {
 155 /* Computing MIN */
 156         i__1 = 1, i__2 = *n / 2;
 157         if (min(i__1,i__2) > *cutpnt || *n / 2 < *cutpnt) {
 158             *info = -7;
 159         }
 160     }
 161     if (*info != 0) {
 162         i__1 = -(*info);
 163         xerbla_("DLAED1", &i__1);
 164         return 0;
 165     }
 166
 167 /*     Quick return if possible */
 168
 169     if (*n == 0) {
 170         return 0;
 171     }
 172
 173 /*     The following values are integer pointers which indicate */
 174 /*     the portion of the workspace */
 175 /*     used by a particular array in DLAED2 and DLAED3. */
 176
 177     iz = 1;
 178     idlmda = iz + *n;
 179     iw = idlmda + *n;
 180     iq2 = iw + *n;
 181
 182     indx = 1;
 183     indxc = indx + *n;
 184     coltyp = indxc + *n;
 185     indxp = coltyp + *n;
 186
 187
 188 /*     Form the z-vector which consists of the last row of Q_1 and the */
 189 /*     first row of Q_2. */
 190
 191     dcopy_(cutpnt, &q[*cutpnt + q_dim1], ldq, &work[iz], &c__1);
 192     zpp1 = *cutpnt + 1;
 193     i__1 = *n - *cutpnt;
 194     dcopy_(&i__1, &q[zpp1 + zpp1 * q_dim1], ldq, &work[iz + *cutpnt], &c__1);
 195
 196 /*     Deflate eigenvalues. */
 197
 198     dlaed2_(&k, n, cutpnt, &d__[1], &q[q_offset], ldq, &indxq[1], rho, &work[
 199             iz], &work[idlmda], &work[iw], &work[iq2], &iwork[indx], &iwork[
 200             indxc], &iwork[indxp], &iwork[coltyp], info);
 201
 202     if (*info != 0) {
 203         goto L20;
 204     }
 205
 206 /*     Solve Secular Equation. */
 207
 208     if (k != 0) {
 209         is = (iwork[coltyp] + iwork[coltyp + 1]) * *cutpnt + (iwork[coltyp +
 210                 1] + iwork[coltyp + 2]) * (*n - *cutpnt) + iq2;
 211         dlaed3_(&k, n, cutpnt, &d__[1], &q[q_offset], ldq, rho, &work[idlmda],
 212                  &work[iq2], &iwork[indxc], &iwork[coltyp], &work[iw], &work[
 213                 is], info);
 214         if (*info != 0) {
 215             goto L20;
 216         }
 217
 218 /*     Prepare the INDXQ sorting permutation. */
 219
 220         n1 = k;
 221         n2 = *n - k;
 222         dlamrg_(&n1, &n2, &d__[1], &c__1, &c_n1, &indxq[1]);
 223     } else {
 224         i__1 = *n;
 225         for (i__ = 1; i__ <= i__1; ++i__) {
 226             indxq[i__] = i__;
 227 /* L10: */
 228         }
 229     }
 230
 231 L20:
 232     return 0;
 233
 234 /*     End of DLAED1 */
 235
 236 } /* dlaed1_ */