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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__2 = 2;
   6 static integer c__1 = 1;
   7 static real c_b24 = 1.f;
   8 static real c_b26 = 0.f;
   9
  10 /* Subroutine */ int slaeda_(integer *n, integer *tlvls, integer *curlvl,
  11         integer *curpbm, integer *prmptr, integer *perm, integer *givptr,
  12         integer *givcol, real *givnum, real *q, integer *qptr, real *z__,
  13         real *ztemp, integer *info)
  14 {
  15     /* System generated locals */
  16     integer i__1, i__2, i__3;
  17
  18     /* Builtin functions */
  19     integer pow_ii(integer *, integer *);
  20     double sqrt(doublereal);
  21
  22     /* Local variables */
  23     integer i__, k, mid, ptr, curr;
  24     extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
  25             integer *, real *, real *);
  26     integer bsiz1, bsiz2, psiz1, psiz2, zptr1;
  27     extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *,
  28             real *, integer *, real *, integer *, real *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *),
  29             xerbla_(char *, integer *);
  30
  31
  32 /*  -- LAPACK routine (version 3.1) -- */
  33 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  34 /*     November 2006 */
  35
  36 /*     .. Scalar Arguments .. */
  37 /*     .. */
  38 /*     .. Array Arguments .. */
  39 /*     .. */
  40
  41 /*  Purpose */
  42 /*  ======= */
  43
  44 /*  SLAEDA computes the Z vector corresponding to the merge step in the */
  45 /*  CURLVLth step of the merge process with TLVLS steps for the CURPBMth */
  46 /*  problem. */
  47
  48 /*  Arguments */
  49 /*  ========= */
  50
  51 /*  N      (input) INTEGER */
  52 /*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */
  53
  54 /*  TLVLS  (input) INTEGER */
  55 /*         The total number of merging levels in the overall divide and */
  56 /*         conquer tree. */
  57
  58 /*  CURLVL (input) INTEGER */
  59 /*         The current level in the overall merge routine, */
  60 /*         0 <= curlvl <= tlvls. */
  61
  62 /*  CURPBM (input) INTEGER */
  63 /*         The current problem in the current level in the overall */
  64 /*         merge routine (counting from upper left to lower right). */
  65
  66 /*  PRMPTR (input) INTEGER array, dimension (N lg N) */
  67 /*         Contains a list of pointers which indicate where in PERM a */
  68 /*         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i) */
  69 /*         indicates the size of the permutation and incidentally the */
  70 /*         size of the full, non-deflated problem. */
  71
  72 /*  PERM   (input) INTEGER array, dimension (N lg N) */
  73 /*         Contains the permutations (from deflation and sorting) to be */
  74 /*         applied to each eigenblock. */
  75
  76 /*  GIVPTR (input) INTEGER array, dimension (N lg N) */
  77 /*         Contains a list of pointers which indicate where in GIVCOL a */
  78 /*         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i) */
  79 /*         indicates the number of Givens rotations. */
  80
  81 /*  GIVCOL (input) INTEGER array, dimension (2, N lg N) */
  82 /*         Each pair of numbers indicates a pair of columns to take place */
  83 /*         in a Givens rotation. */
  84
  85 /*  GIVNUM (input) REAL array, dimension (2, N lg N) */
  86 /*         Each number indicates the S value to be used in the */
  87 /*         corresponding Givens rotation. */
  88
  89 /*  Q      (input) REAL array, dimension (N**2) */
  90 /*         Contains the square eigenblocks from previous levels, the */
  91 /*         starting positions for blocks are given by QPTR. */
  92
  93 /*  QPTR   (input) INTEGER array, dimension (N+2) */
  94 /*         Contains a list of pointers which indicate where in Q an */
  95 /*         eigenblock is stored.  SQRT( QPTR(i+1) - QPTR(i) ) indicates */
  96 /*         the size of the block. */
  97
  98 /*  Z      (output) REAL array, dimension (N) */
  99 /*         On output this vector contains the updating vector (the last */
 100 /*         row of the first sub-eigenvector matrix and the first row of */
 101 /*         the second sub-eigenvector matrix). */
 102
 103 /*  ZTEMP  (workspace) REAL array, dimension (N) */
 104
 105 /*  INFO   (output) INTEGER */
 106 /*          = 0:  successful exit. */
 107 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
 108
 109 /*  Further Details */
 110 /*  =============== */
 111
 112 /*  Based on contributions by */
 113 /*     Jeff Rutter, Computer Science Division, University of California */
 114 /*     at Berkeley, USA */
 115
 116 /*  ===================================================================== */
 117
 118 /*     .. Parameters .. */
 119 /*     .. */
 120 /*     .. Local Scalars .. */
 121 /*     .. */
 122 /*     .. External Subroutines .. */
 123 /*     .. */
 124 /*     .. Intrinsic Functions .. */
 125 /*     .. */
 126 /*     .. Executable Statements .. */
 127
 128 /*     Test the input parameters. */
 129
 130     /* Parameter adjustments */
 131     --ztemp;
 132     --z__;
 133     --qptr;
 134     --q;
 135     givnum -= 3;
 136     givcol -= 3;
 137     --givptr;
 138     --perm;
 139     --prmptr;
 140
 141     /* Function Body */
 142     *info = 0;
 143
 144     if (*n < 0) {
 145         *info = -1;
 146     }
 147     if (*info != 0) {
 148         i__1 = -(*info);
 149         xerbla_("SLAEDA", &i__1);
 150         return 0;
 151     }
 152
 153 /*     Quick return if possible */
 154
 155     if (*n == 0) {
 156         return 0;
 157     }
 158
 159 /*     Determine location of first number in second half. */
 160
 161     mid = *n / 2 + 1;
 162
 163 /*     Gather last/first rows of appropriate eigenblocks into center of Z */
 164
 165     ptr = 1;
 166
 167 /*     Determine location of lowest level subproblem in the full storage */
 168 /*     scheme */
 169
 170     i__1 = *curlvl - 1;
 171     curr = ptr + *curpbm * pow_ii(&c__2, curlvl) + pow_ii(&c__2, &i__1) - 1;
 172
 173 /*     Determine size of these matrices.  We add HALF to the value of */
 174 /*     the SQRT in case the machine underestimates one of these square */
 175 /*     roots. */
 176
 177     bsiz1 = (integer) (sqrt((real) (qptr[curr + 1] - qptr[curr])) + .5f);
 178     bsiz2 = (integer) (sqrt((real) (qptr[curr + 2] - qptr[curr + 1])) + .5f);
 179     i__1 = mid - bsiz1 - 1;
 180     for (k = 1; k <= i__1; ++k) {
 181         z__[k] = 0.f;
 182 /* L10: */
 183     }
 184     scopy_(&bsiz1, &q[qptr[curr] + bsiz1 - 1], &bsiz1, &z__[mid - bsiz1], &
 185             c__1);
 186     scopy_(&bsiz2, &q[qptr[curr + 1]], &bsiz2, &z__[mid], &c__1);
 187     i__1 = *n;
 188     for (k = mid + bsiz2; k <= i__1; ++k) {
 189         z__[k] = 0.f;
 190 /* L20: */
 191     }
 192
 193 /*     Loop thru remaining levels 1 -> CURLVL applying the Givens */
 194 /*     rotations and permutation and then multiplying the center matrices */
 195 /*     against the current Z. */
 196
 197     ptr = pow_ii(&c__2, tlvls) + 1;
 198     i__1 = *curlvl - 1;
 199     for (k = 1; k <= i__1; ++k) {
 200         i__2 = *curlvl - k;
 201         i__3 = *curlvl - k - 1;
 202         curr = ptr + *curpbm * pow_ii(&c__2, &i__2) + pow_ii(&c__2, &i__3) -
 203                 1;
 204         psiz1 = prmptr[curr + 1] - prmptr[curr];
 205         psiz2 = prmptr[curr + 2] - prmptr[curr + 1];
 206         zptr1 = mid - psiz1;
 207
 208 /*       Apply Givens at CURR and CURR+1 */
 209
 210         i__2 = givptr[curr + 1] - 1;
 211         for (i__ = givptr[curr]; i__ <= i__2; ++i__) {
 212             srot_(&c__1, &z__[zptr1 + givcol[(i__ << 1) + 1] - 1], &c__1, &
 213                     z__[zptr1 + givcol[(i__ << 1) + 2] - 1], &c__1, &givnum[(
 214                     i__ << 1) + 1], &givnum[(i__ << 1) + 2]);
 215 /* L30: */
 216         }
 217         i__2 = givptr[curr + 2] - 1;
 218         for (i__ = givptr[curr + 1]; i__ <= i__2; ++i__) {
 219             srot_(&c__1, &z__[mid - 1 + givcol[(i__ << 1) + 1]], &c__1, &z__[
 220                     mid - 1 + givcol[(i__ << 1) + 2]], &c__1, &givnum[(i__ <<
 221                     1) + 1], &givnum[(i__ << 1) + 2]);
 222 /* L40: */
 223         }
 224         psiz1 = prmptr[curr + 1] - prmptr[curr];
 225         psiz2 = prmptr[curr + 2] - prmptr[curr + 1];
 226         i__2 = psiz1 - 1;
 227         for (i__ = 0; i__ <= i__2; ++i__) {
 228             ztemp[i__ + 1] = z__[zptr1 + perm[prmptr[curr] + i__] - 1];
 229 /* L50: */
 230         }
 231         i__2 = psiz2 - 1;
 232         for (i__ = 0; i__ <= i__2; ++i__) {
 233             ztemp[psiz1 + i__ + 1] = z__[mid + perm[prmptr[curr + 1] + i__] -
 234                     1];
 235 /* L60: */
 236         }
 237
 238 /*        Multiply Blocks at CURR and CURR+1 */
 239
 240 /*        Determine size of these matrices.  We add HALF to the value of */
 241 /*        the SQRT in case the machine underestimates one of these */
 242 /*        square roots. */
 243
 244         bsiz1 = (integer) (sqrt((real) (qptr[curr + 1] - qptr[curr])) + .5f);
 245         bsiz2 = (integer) (sqrt((real) (qptr[curr + 2] - qptr[curr + 1])) +
 246                 .5f);
 247         if (bsiz1 > 0) {
 248             sgemv_("T", &bsiz1, &bsiz1, &c_b24, &q[qptr[curr]], &bsiz1, &
 249                     ztemp[1], &c__1, &c_b26, &z__[zptr1], &c__1);
 250         }
 251         i__2 = psiz1 - bsiz1;
 252         scopy_(&i__2, &ztemp[bsiz1 + 1], &c__1, &z__[zptr1 + bsiz1], &c__1);
 253         if (bsiz2 > 0) {
 254             sgemv_("T", &bsiz2, &bsiz2, &c_b24, &q[qptr[curr + 1]], &bsiz2, &
 255                     ztemp[psiz1 + 1], &c__1, &c_b26, &z__[mid], &c__1);
 256         }
 257         i__2 = psiz2 - bsiz2;
 258         scopy_(&i__2, &ztemp[psiz1 + bsiz2 + 1], &c__1, &z__[mid + bsiz2], &
 259                 c__1);
 260
 261         i__2 = *tlvls - k;
 262         ptr += pow_ii(&c__2, &i__2);
 263 /* L70: */
 264     }
 265
 266     return 0;
 267
 268 /*     End of SLAEDA */
 269
 270 } /* slaeda_ */