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