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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static doublereal c_b4 = -1.;
   6 static doublereal c_b5 = 1.;
   7 static integer c__1 = 1;
   8 static doublereal c_b16 = 0.;
   9
  10 /* Subroutine */ int dlabrd_(integer *m, integer *n, integer *nb, doublereal *
  11         a, integer *lda, doublereal *d__, doublereal *e, doublereal *tauq,
  12         doublereal *taup, doublereal *x, integer *ldx, doublereal *y, integer
  13         *ldy)
  14 {
  15     /* System generated locals */
  16     integer a_dim1, a_offset, x_dim1, x_offset, y_dim1, y_offset, i__1, i__2,
  17             i__3;
  18
  19     /* Local variables */
  20     integer i__;
  21     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  22             integer *), dgemv_(char *, integer *, integer *, doublereal *,
  23             doublereal *, integer *, doublereal *, integer *, doublereal *,
  24             doublereal *, integer *), dlarfg_(integer *, doublereal *,
  25              doublereal *, integer *, doublereal *);
  26
  27
  28 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  29 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  30 /*     November 2006 */
  31
  32 /*     .. Scalar Arguments .. */
  33 /*     .. */
  34 /*     .. Array Arguments .. */
  35 /*     .. */
  36
  37 /*  Purpose */
  38 /*  ======= */
  39
  40 /*  DLABRD reduces the first NB rows and columns of a real general */
  41 /*  m by n matrix A to upper or lower bidiagonal form by an orthogonal */
  42 /*  transformation Q' * A * P, and returns the matrices X and Y which */
  43 /*  are needed to apply the transformation to the unreduced part of A. */
  44
  45 /*  If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower */
  46 /*  bidiagonal form. */
  47
  48 /*  This is an auxiliary routine called by DGEBRD */
  49
  50 /*  Arguments */
  51 /*  ========= */
  52
  53 /*  M       (input) INTEGER */
  54 /*          The number of rows in the matrix A. */
  55
  56 /*  N       (input) INTEGER */
  57 /*          The number of columns in the matrix A. */
  58
  59 /*  NB      (input) INTEGER */
  60 /*          The number of leading rows and columns of A to be reduced. */
  61
  62 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
  63 /*          On entry, the m by n general matrix to be reduced. */
  64 /*          On exit, the first NB rows and columns of the matrix are */
  65 /*          overwritten; the rest of the array is unchanged. */
  66 /*          If m >= n, elements on and below the diagonal in the first NB */
  67 /*            columns, with the array TAUQ, represent the orthogonal */
  68 /*            matrix Q as a product of elementary reflectors; and */
  69 /*            elements above the diagonal in the first NB rows, with the */
  70 /*            array TAUP, represent the orthogonal matrix P as a product */
  71 /*            of elementary reflectors. */
  72 /*          If m < n, elements below the diagonal in the first NB */
  73 /*            columns, with the array TAUQ, represent the orthogonal */
  74 /*            matrix Q as a product of elementary reflectors, and */
  75 /*            elements on and above the diagonal in the first NB rows, */
  76 /*            with the array TAUP, represent the orthogonal matrix P as */
  77 /*            a product of elementary reflectors. */
  78 /*          See Further Details. */
  79
  80 /*  LDA     (input) INTEGER */
  81 /*          The leading dimension of the array A.  LDA >= max(1,M). */
  82
  83 /*  D       (output) DOUBLE PRECISION array, dimension (NB) */
  84 /*          The diagonal elements of the first NB rows and columns of */
  85 /*          the reduced matrix.  D(i) = A(i,i). */
  86
  87 /*  E       (output) DOUBLE PRECISION array, dimension (NB) */
  88 /*          The off-diagonal elements of the first NB rows and columns of */
  89 /*          the reduced matrix. */
  90
  91 /*  TAUQ    (output) DOUBLE PRECISION array dimension (NB) */
  92 /*          The scalar factors of the elementary reflectors which */
  93 /*          represent the orthogonal matrix Q. See Further Details. */
  94
  95 /*  TAUP    (output) DOUBLE PRECISION array, dimension (NB) */
  96 /*          The scalar factors of the elementary reflectors which */
  97 /*          represent the orthogonal matrix P. See Further Details. */
  98
  99 /*  X       (output) DOUBLE PRECISION array, dimension (LDX,NB) */
 100 /*          The m-by-nb matrix X required to update the unreduced part */
 101 /*          of A. */
 102
 103 /*  LDX     (input) INTEGER */
 104 /*          The leading dimension of the array X. LDX >= M. */
 105
 106 /*  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB) */
 107 /*          The n-by-nb matrix Y required to update the unreduced part */
 108 /*          of A. */
 109
 110 /*  LDY     (input) INTEGER */
 111 /*          The leading dimension of the array Y. LDY >= N. */
 112
 113 /*  Further Details */
 114 /*  =============== */
 115
 116 /*  The matrices Q and P are represented as products of elementary */
 117 /*  reflectors: */
 118
 119 /*     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb) */
 120
 121 /*  Each H(i) and G(i) has the form: */
 122
 123 /*     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u' */
 124
 125 /*  where tauq and taup are real scalars, and v and u are real vectors. */
 126
 127 /*  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in */
 128 /*  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in */
 129 /*  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). */
 130
 131 /*  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in */
 132 /*  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in */
 133 /*  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). */
 134
 135 /*  The elements of the vectors v and u together form the m-by-nb matrix */
 136 /*  V and the nb-by-n matrix U' which are needed, with X and Y, to apply */
 137 /*  the transformation to the unreduced part of the matrix, using a block */
 138 /*  update of the form:  A := A - V*Y' - X*U'. */
 139
 140 /*  The contents of A on exit are illustrated by the following examples */
 141 /*  with nb = 2: */
 142
 143 /*  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n): */
 144
 145 /*    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 ) */
 146 /*    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 ) */
 147 /*    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  ) */
 148 /*    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  ) */
 149 /*    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  ) */
 150 /*    (  v1  v2  a   a   a  ) */
 151
 152 /*  where a denotes an element of the original matrix which is unchanged, */
 153 /*  vi denotes an element of the vector defining H(i), and ui an element */
 154 /*  of the vector defining G(i). */
 155
 156 /*  ===================================================================== */
 157
 158 /*     .. Parameters .. */
 159 /*     .. */
 160 /*     .. Local Scalars .. */
 161 /*     .. */
 162 /*     .. External Subroutines .. */
 163 /*     .. */
 164 /*     .. Intrinsic Functions .. */
 165 /*     .. */
 166 /*     .. Executable Statements .. */
 167
 168 /*     Quick return if possible */
 169
 170     /* Parameter adjustments */
 171     a_dim1 = *lda;
 172     a_offset = 1 + a_dim1;
 173     a -= a_offset;
 174     --d__;
 175     --e;
 176     --tauq;
 177     --taup;
 178     x_dim1 = *ldx;
 179     x_offset = 1 + x_dim1;
 180     x -= x_offset;
 181     y_dim1 = *ldy;
 182     y_offset = 1 + y_dim1;
 183     y -= y_offset;
 184
 185     /* Function Body */
 186     if (*m <= 0 || *n <= 0) {
 187         return 0;
 188     }
 189
 190     if (*m >= *n) {
 191
 192 /*        Reduce to upper bidiagonal form */
 193
 194         i__1 = *nb;
 195         for (i__ = 1; i__ <= i__1; ++i__) {
 196
 197 /*           Update A(i:m,i) */
 198
 199             i__2 = *m - i__ + 1;
 200             i__3 = i__ - 1;
 201             dgemv_("No transpose", &i__2, &i__3, &c_b4, &a[i__ + a_dim1], lda,
 202                      &y[i__ + y_dim1], ldy, &c_b5, &a[i__ + i__ * a_dim1], &
 203                     c__1);
 204             i__2 = *m - i__ + 1;
 205             i__3 = i__ - 1;
 206             dgemv_("No transpose", &i__2, &i__3, &c_b4, &x[i__ + x_dim1], ldx,
 207                      &a[i__ * a_dim1 + 1], &c__1, &c_b5, &a[i__ + i__ *
 208                     a_dim1], &c__1);
 209
 210 /*           Generate reflection Q(i) to annihilate A(i+1:m,i) */
 211
 212             i__2 = *m - i__ + 1;
 213 /* Computing MIN */
 214             i__3 = i__ + 1;
 215             dlarfg_(&i__2, &a[i__ + i__ * a_dim1], &a[min(i__3, *m)+ i__ *
 216                     a_dim1], &c__1, &tauq[i__]);
 217             d__[i__] = a[i__ + i__ * a_dim1];
 218             if (i__ < *n) {
 219                 a[i__ + i__ * a_dim1] = 1.;
 220
 221 /*              Compute Y(i+1:n,i) */
 222
 223                 i__2 = *m - i__ + 1;
 224                 i__3 = *n - i__;
 225                 dgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + (i__ + 1) *
 226                         a_dim1], lda, &a[i__ + i__ * a_dim1], &c__1, &c_b16, &
 227                         y[i__ + 1 + i__ * y_dim1], &c__1);
 228                 i__2 = *m - i__ + 1;
 229                 i__3 = i__ - 1;
 230                 dgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + a_dim1],
 231                         lda, &a[i__ + i__ * a_dim1], &c__1, &c_b16, &y[i__ *
 232                         y_dim1 + 1], &c__1);
 233                 i__2 = *n - i__;
 234                 i__3 = i__ - 1;
 235                 dgemv_("No transpose", &i__2, &i__3, &c_b4, &y[i__ + 1 +
 236                         y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &c_b5, &y[
 237                         i__ + 1 + i__ * y_dim1], &c__1);
 238                 i__2 = *m - i__ + 1;
 239                 i__3 = i__ - 1;
 240                 dgemv_("Transpose", &i__2, &i__3, &c_b5, &x[i__ + x_dim1],
 241                         ldx, &a[i__ + i__ * a_dim1], &c__1, &c_b16, &y[i__ *
 242                         y_dim1 + 1], &c__1);
 243                 i__2 = i__ - 1;
 244                 i__3 = *n - i__;
 245                 dgemv_("Transpose", &i__2, &i__3, &c_b4, &a[(i__ + 1) *
 246                         a_dim1 + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &c_b5,
 247                         &y[i__ + 1 + i__ * y_dim1], &c__1);
 248                 i__2 = *n - i__;
 249                 dscal_(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
 250
 251 /*              Update A(i,i+1:n) */
 252
 253                 i__2 = *n - i__;
 254                 dgemv_("No transpose", &i__2, &i__, &c_b4, &y[i__ + 1 +
 255                         y_dim1], ldy, &a[i__ + a_dim1], lda, &c_b5, &a[i__ + (
 256                         i__ + 1) * a_dim1], lda);
 257                 i__2 = i__ - 1;
 258                 i__3 = *n - i__;
 259                 dgemv_("Transpose", &i__2, &i__3, &c_b4, &a[(i__ + 1) *
 260                         a_dim1 + 1], lda, &x[i__ + x_dim1], ldx, &c_b5, &a[
 261                         i__ + (i__ + 1) * a_dim1], lda);
 262
 263 /*              Generate reflection P(i) to annihilate A(i,i+2:n) */
 264
 265                 i__2 = *n - i__;
 266 /* Computing MIN */
 267                 i__3 = i__ + 2;
 268                 dlarfg_(&i__2, &a[i__ + (i__ + 1) * a_dim1], &a[i__ + min(
 269                         i__3, *n)* a_dim1], lda, &taup[i__]);
 270                 e[i__] = a[i__ + (i__ + 1) * a_dim1];
 271                 a[i__ + (i__ + 1) * a_dim1] = 1.;
 272
 273 /*              Compute X(i+1:m,i) */
 274
 275                 i__2 = *m - i__;
 276                 i__3 = *n - i__;
 277                 dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + (i__
 278                         + 1) * a_dim1], lda, &a[i__ + (i__ + 1) * a_dim1],
 279                         lda, &c_b16, &x[i__ + 1 + i__ * x_dim1], &c__1);
 280                 i__2 = *n - i__;
 281                 dgemv_("Transpose", &i__2, &i__, &c_b5, &y[i__ + 1 + y_dim1],
 282                         ldy, &a[i__ + (i__ + 1) * a_dim1], lda, &c_b16, &x[
 283                         i__ * x_dim1 + 1], &c__1);
 284                 i__2 = *m - i__;
 285                 dgemv_("No transpose", &i__2, &i__, &c_b4, &a[i__ + 1 +
 286                         a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
 287                         i__ + 1 + i__ * x_dim1], &c__1);
 288                 i__2 = i__ - 1;
 289                 i__3 = *n - i__;
 290                 dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[(i__ + 1) *
 291                         a_dim1 + 1], lda, &a[i__ + (i__ + 1) * a_dim1], lda, &
 292                         c_b16, &x[i__ * x_dim1 + 1], &c__1);
 293                 i__2 = *m - i__;
 294                 i__3 = i__ - 1;
 295                 dgemv_("No transpose", &i__2, &i__3, &c_b4, &x[i__ + 1 +
 296                         x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
 297                         i__ + 1 + i__ * x_dim1], &c__1);
 298                 i__2 = *m - i__;
 299                 dscal_(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
 300             }
 301 /* L10: */
 302         }
 303     } else {
 304
 305 /*        Reduce to lower bidiagonal form */
 306
 307         i__1 = *nb;
 308         for (i__ = 1; i__ <= i__1; ++i__) {
 309
 310 /*           Update A(i,i:n) */
 311
 312             i__2 = *n - i__ + 1;
 313             i__3 = i__ - 1;
 314             dgemv_("No transpose", &i__2, &i__3, &c_b4, &y[i__ + y_dim1], ldy,
 315                      &a[i__ + a_dim1], lda, &c_b5, &a[i__ + i__ * a_dim1],
 316                     lda);
 317             i__2 = i__ - 1;
 318             i__3 = *n - i__ + 1;
 319             dgemv_("Transpose", &i__2, &i__3, &c_b4, &a[i__ * a_dim1 + 1],
 320                     lda, &x[i__ + x_dim1], ldx, &c_b5, &a[i__ + i__ * a_dim1],
 321                      lda);
 322
 323 /*           Generate reflection P(i) to annihilate A(i,i+1:n) */
 324
 325             i__2 = *n - i__ + 1;
 326 /* Computing MIN */
 327             i__3 = i__ + 1;
 328             dlarfg_(&i__2, &a[i__ + i__ * a_dim1], &a[i__ + min(i__3, *n)*
 329                     a_dim1], lda, &taup[i__]);
 330             d__[i__] = a[i__ + i__ * a_dim1];
 331             if (i__ < *m) {
 332                 a[i__ + i__ * a_dim1] = 1.;
 333
 334 /*              Compute X(i+1:m,i) */
 335
 336                 i__2 = *m - i__;
 337                 i__3 = *n - i__ + 1;
 338                 dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + i__ *
 339                          a_dim1], lda, &a[i__ + i__ * a_dim1], lda, &c_b16, &
 340                         x[i__ + 1 + i__ * x_dim1], &c__1);
 341                 i__2 = *n - i__ + 1;
 342                 i__3 = i__ - 1;
 343                 dgemv_("Transpose", &i__2, &i__3, &c_b5, &y[i__ + y_dim1],
 344                         ldy, &a[i__ + i__ * a_dim1], lda, &c_b16, &x[i__ *
 345                         x_dim1 + 1], &c__1);
 346                 i__2 = *m - i__;
 347                 i__3 = i__ - 1;
 348                 dgemv_("No transpose", &i__2, &i__3, &c_b4, &a[i__ + 1 +
 349                         a_dim1], lda, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
 350                         i__ + 1 + i__ * x_dim1], &c__1);
 351                 i__2 = i__ - 1;
 352                 i__3 = *n - i__ + 1;
 353                 dgemv_("No transpose", &i__2, &i__3, &c_b5, &a[i__ * a_dim1 +
 354                         1], lda, &a[i__ + i__ * a_dim1], lda, &c_b16, &x[i__ *
 355                          x_dim1 + 1], &c__1);
 356                 i__2 = *m - i__;
 357                 i__3 = i__ - 1;
 358                 dgemv_("No transpose", &i__2, &i__3, &c_b4, &x[i__ + 1 +
 359                         x_dim1], ldx, &x[i__ * x_dim1 + 1], &c__1, &c_b5, &x[
 360                         i__ + 1 + i__ * x_dim1], &c__1);
 361                 i__2 = *m - i__;
 362                 dscal_(&i__2, &taup[i__], &x[i__ + 1 + i__ * x_dim1], &c__1);
 363
 364 /*              Update A(i+1:m,i) */
 365
 366                 i__2 = *m - i__;
 367                 i__3 = i__ - 1;
 368                 dgemv_("No transpose", &i__2, &i__3, &c_b4, &a[i__ + 1 +
 369                         a_dim1], lda, &y[i__ + y_dim1], ldy, &c_b5, &a[i__ +
 370                         1 + i__ * a_dim1], &c__1);
 371                 i__2 = *m - i__;
 372                 dgemv_("No transpose", &i__2, &i__, &c_b4, &x[i__ + 1 +
 373                         x_dim1], ldx, &a[i__ * a_dim1 + 1], &c__1, &c_b5, &a[
 374                         i__ + 1 + i__ * a_dim1], &c__1);
 375
 376 /*              Generate reflection Q(i) to annihilate A(i+2:m,i) */
 377
 378                 i__2 = *m - i__;
 379 /* Computing MIN */
 380                 i__3 = i__ + 2;
 381                 dlarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *m)+
 382                         i__ * a_dim1], &c__1, &tauq[i__]);
 383                 e[i__] = a[i__ + 1 + i__ * a_dim1];
 384                 a[i__ + 1 + i__ * a_dim1] = 1.;
 385
 386 /*              Compute Y(i+1:n,i) */
 387
 388                 i__2 = *m - i__;
 389                 i__3 = *n - i__;
 390                 dgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + (i__ +
 391                         1) * a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1,
 392                         &c_b16, &y[i__ + 1 + i__ * y_dim1], &c__1);
 393                 i__2 = *m - i__;
 394                 i__3 = i__ - 1;
 395                 dgemv_("Transpose", &i__2, &i__3, &c_b5, &a[i__ + 1 + a_dim1],
 396                          lda, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &y[
 397                         i__ * y_dim1 + 1], &c__1);
 398                 i__2 = *n - i__;
 399                 i__3 = i__ - 1;
 400                 dgemv_("No transpose", &i__2, &i__3, &c_b4, &y[i__ + 1 +
 401                         y_dim1], ldy, &y[i__ * y_dim1 + 1], &c__1, &c_b5, &y[
 402                         i__ + 1 + i__ * y_dim1], &c__1);
 403                 i__2 = *m - i__;
 404                 dgemv_("Transpose", &i__2, &i__, &c_b5, &x[i__ + 1 + x_dim1],
 405                         ldx, &a[i__ + 1 + i__ * a_dim1], &c__1, &c_b16, &y[
 406                         i__ * y_dim1 + 1], &c__1);
 407                 i__2 = *n - i__;
 408                 dgemv_("Transpose", &i__, &i__2, &c_b4, &a[(i__ + 1) * a_dim1
 409                         + 1], lda, &y[i__ * y_dim1 + 1], &c__1, &c_b5, &y[i__
 410                         + 1 + i__ * y_dim1], &c__1);
 411                 i__2 = *n - i__;
 412                 dscal_(&i__2, &tauq[i__], &y[i__ + 1 + i__ * y_dim1], &c__1);
 413             }
 414 /* L20: */
 415         }
 416     }
 417     return 0;
 418
 419 /*     End of DLABRD */
 420
 421 } /* dlabrd_ */