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