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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6
   7 /* Subroutine */ int dgebd2_(integer *m, integer *n, doublereal *a, integer *
   8         lda, doublereal *d__, doublereal *e, doublereal *tauq, doublereal *
   9         taup, doublereal *work, integer *info)
  10 {
  11     /* System generated locals */
  12     integer a_dim1, a_offset, i__1, i__2, i__3;
  13
  14     /* Local variables */
  15     integer i__;
  16     extern /* Subroutine */ int dlarf_(char *, integer *, integer *,
  17             doublereal *, integer *, doublereal *, doublereal *, integer *,
  18             doublereal *), dlarfg_(integer *, doublereal *,
  19             doublereal *, integer *, doublereal *), xerbla_(char *, integer *);
  20
  21
  22 /*  -- LAPACK routine (version 3.1) -- */
  23 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  24 /*     November 2006 */
  25
  26 /*     .. Scalar Arguments .. */
  27 /*     .. */
  28 /*     .. Array Arguments .. */
  29 /*     .. */
  30
  31 /*  Purpose */
  32 /*  ======= */
  33
  34 /*  DGEBD2 reduces a real general m by n matrix A to upper or lower */
  35 /*  bidiagonal form B by an orthogonal transformation: Q' * A * P = B. */
  36
  37 /*  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. */
  38
  39 /*  Arguments */
  40 /*  ========= */
  41
  42 /*  M       (input) INTEGER */
  43 /*          The number of rows in the matrix A.  M >= 0. */
  44
  45 /*  N       (input) INTEGER */
  46 /*          The number of columns in the matrix A.  N >= 0. */
  47
  48 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
  49 /*          On entry, the m by n general matrix to be reduced. */
  50 /*          On exit, */
  51 /*          if m >= n, the diagonal and the first superdiagonal are */
  52 /*            overwritten with the upper bidiagonal matrix B; the */
  53 /*            elements below the diagonal, with the array TAUQ, represent */
  54 /*            the orthogonal matrix Q as a product of elementary */
  55 /*            reflectors, and the elements above the first superdiagonal, */
  56 /*            with the array TAUP, represent the orthogonal matrix P as */
  57 /*            a product of elementary reflectors; */
  58 /*          if m < n, the diagonal and the first subdiagonal are */
  59 /*            overwritten with the lower bidiagonal matrix B; the */
  60 /*            elements below the first subdiagonal, with the array TAUQ, */
  61 /*            represent the orthogonal matrix Q as a product of */
  62 /*            elementary reflectors, and the elements above the diagonal, */
  63 /*            with the array TAUP, represent the orthogonal matrix P as */
  64 /*            a product of elementary reflectors. */
  65 /*          See Further Details. */
  66
  67 /*  LDA     (input) INTEGER */
  68 /*          The leading dimension of the array A.  LDA >= max(1,M). */
  69
  70 /*  D       (output) DOUBLE PRECISION array, dimension (min(M,N)) */
  71 /*          The diagonal elements of the bidiagonal matrix B: */
  72 /*          D(i) = A(i,i). */
  73
  74 /*  E       (output) DOUBLE PRECISION array, dimension (min(M,N)-1) */
  75 /*          The off-diagonal elements of the bidiagonal matrix B: */
  76 /*          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; */
  77 /*          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. */
  78
  79 /*  TAUQ    (output) DOUBLE PRECISION array dimension (min(M,N)) */
  80 /*          The scalar factors of the elementary reflectors which */
  81 /*          represent the orthogonal matrix Q. See Further Details. */
  82
  83 /*  TAUP    (output) DOUBLE PRECISION array, dimension (min(M,N)) */
  84 /*          The scalar factors of the elementary reflectors which */
  85 /*          represent the orthogonal matrix P. See Further Details. */
  86
  87 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (max(M,N)) */
  88
  89 /*  INFO    (output) INTEGER */
  90 /*          = 0: successful exit. */
  91 /*          < 0: if INFO = -i, the i-th argument had an illegal value. */
  92
  93 /*  Further Details */
  94 /*  =============== */
  95
  96 /*  The matrices Q and P are represented as products of elementary */
  97 /*  reflectors: */
  98
  99 /*  If m >= n, */
 100
 101 /*     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1) */
 102
 103 /*  Each H(i) and G(i) has the form: */
 104
 105 /*     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u' */
 106
 107 /*  where tauq and taup are real scalars, and v and u are real vectors; */
 108 /*  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i); */
 109 /*  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n); */
 110 /*  tauq is stored in TAUQ(i) and taup in TAUP(i). */
 111
 112 /*  If m < n, */
 113
 114 /*     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m) */
 115
 116 /*  Each H(i) and G(i) has the form: */
 117
 118 /*     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u' */
 119
 120 /*  where tauq and taup are real scalars, and v and u are real vectors; */
 121 /*  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i); */
 122 /*  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n); */
 123 /*  tauq is stored in TAUQ(i) and taup in TAUP(i). */
 124
 125 /*  The contents of A on exit are illustrated by the following examples: */
 126
 127 /*  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n): */
 128
 129 /*    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 ) */
 130 /*    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 ) */
 131 /*    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 ) */
 132 /*    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 ) */
 133 /*    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 ) */
 134 /*    (  v1  v2  v3  v4  v5 ) */
 135
 136 /*  where d and e denote diagonal and off-diagonal elements of B, vi */
 137 /*  denotes an element of the vector defining H(i), and ui an element of */
 138 /*  the vector defining G(i). */
 139
 140 /*  ===================================================================== */
 141
 142 /*     .. Parameters .. */
 143 /*     .. */
 144 /*     .. Local Scalars .. */
 145 /*     .. */
 146 /*     .. External Subroutines .. */
 147 /*     .. */
 148 /*     .. Intrinsic Functions .. */
 149 /*     .. */
 150 /*     .. Executable Statements .. */
 151
 152 /*     Test the input parameters */
 153
 154     /* Parameter adjustments */
 155     a_dim1 = *lda;
 156     a_offset = 1 + a_dim1;
 157     a -= a_offset;
 158     --d__;
 159     --e;
 160     --tauq;
 161     --taup;
 162     --work;
 163
 164     /* Function Body */
 165     *info = 0;
 166     if (*m < 0) {
 167         *info = -1;
 168     } else if (*n < 0) {
 169         *info = -2;
 170     } else if (*lda < max(1,*m)) {
 171         *info = -4;
 172     }
 173     if (*info < 0) {
 174         i__1 = -(*info);
 175         xerbla_("DGEBD2", &i__1);
 176         return 0;
 177     }
 178
 179     if (*m >= *n) {
 180
 181 /*        Reduce to upper bidiagonal form */
 182
 183         i__1 = *n;
 184         for (i__ = 1; i__ <= i__1; ++i__) {
 185
 186 /*           Generate elementary reflector H(i) to annihilate A(i+1:m,i) */
 187
 188             i__2 = *m - i__ + 1;
 189 /* Computing MIN */
 190             i__3 = i__ + 1;
 191             dlarfg_(&i__2, &a[i__ + i__ * a_dim1], &a[min(i__3, *m)+ i__ *
 192                     a_dim1], &c__1, &tauq[i__]);
 193             d__[i__] = a[i__ + i__ * a_dim1];
 194             a[i__ + i__ * a_dim1] = 1.;
 195
 196 /*           Apply H(i) to A(i:m,i+1:n) from the left */
 197
 198             if (i__ < *n) {
 199                 i__2 = *m - i__ + 1;
 200                 i__3 = *n - i__;
 201                 dlarf_("Left", &i__2, &i__3, &a[i__ + i__ * a_dim1], &c__1, &
 202                         tauq[i__], &a[i__ + (i__ + 1) * a_dim1], lda, &work[1]
 203 );
 204             }
 205             a[i__ + i__ * a_dim1] = d__[i__];
 206
 207             if (i__ < *n) {
 208
 209 /*              Generate elementary reflector G(i) to annihilate */
 210 /*              A(i,i+2:n) */
 211
 212                 i__2 = *n - i__;
 213 /* Computing MIN */
 214                 i__3 = i__ + 2;
 215                 dlarfg_(&i__2, &a[i__ + (i__ + 1) * a_dim1], &a[i__ + min(
 216                         i__3, *n)* a_dim1], lda, &taup[i__]);
 217                 e[i__] = a[i__ + (i__ + 1) * a_dim1];
 218                 a[i__ + (i__ + 1) * a_dim1] = 1.;
 219
 220 /*              Apply G(i) to A(i+1:m,i+1:n) from the right */
 221
 222                 i__2 = *m - i__;
 223                 i__3 = *n - i__;
 224                 dlarf_("Right", &i__2, &i__3, &a[i__ + (i__ + 1) * a_dim1],
 225                         lda, &taup[i__], &a[i__ + 1 + (i__ + 1) * a_dim1],
 226                         lda, &work[1]);
 227                 a[i__ + (i__ + 1) * a_dim1] = e[i__];
 228             } else {
 229                 taup[i__] = 0.;
 230             }
 231 /* L10: */
 232         }
 233     } else {
 234
 235 /*        Reduce to lower bidiagonal form */
 236
 237         i__1 = *m;
 238         for (i__ = 1; i__ <= i__1; ++i__) {
 239
 240 /*           Generate elementary reflector G(i) to annihilate A(i,i+1:n) */
 241
 242             i__2 = *n - i__ + 1;
 243 /* Computing MIN */
 244             i__3 = i__ + 1;
 245             dlarfg_(&i__2, &a[i__ + i__ * a_dim1], &a[i__ + min(i__3, *n)*
 246                     a_dim1], lda, &taup[i__]);
 247             d__[i__] = a[i__ + i__ * a_dim1];
 248             a[i__ + i__ * a_dim1] = 1.;
 249
 250 /*           Apply G(i) to A(i+1:m,i:n) from the right */
 251
 252             if (i__ < *m) {
 253                 i__2 = *m - i__;
 254                 i__3 = *n - i__ + 1;
 255                 dlarf_("Right", &i__2, &i__3, &a[i__ + i__ * a_dim1], lda, &
 256                         taup[i__], &a[i__ + 1 + i__ * a_dim1], lda, &work[1]);
 257             }
 258             a[i__ + i__ * a_dim1] = d__[i__];
 259
 260             if (i__ < *m) {
 261
 262 /*              Generate elementary reflector H(i) to annihilate */
 263 /*              A(i+2:m,i) */
 264
 265                 i__2 = *m - i__;
 266 /* Computing MIN */
 267                 i__3 = i__ + 2;
 268                 dlarfg_(&i__2, &a[i__ + 1 + i__ * a_dim1], &a[min(i__3, *m)+
 269                         i__ * a_dim1], &c__1, &tauq[i__]);
 270                 e[i__] = a[i__ + 1 + i__ * a_dim1];
 271                 a[i__ + 1 + i__ * a_dim1] = 1.;
 272
 273 /*              Apply H(i) to A(i+1:m,i+1:n) from the left */
 274
 275                 i__2 = *m - i__;
 276                 i__3 = *n - i__;
 277                 dlarf_("Left", &i__2, &i__3, &a[i__ + 1 + i__ * a_dim1], &
 278                         c__1, &tauq[i__], &a[i__ + 1 + (i__ + 1) * a_dim1],
 279                         lda, &work[1]);
 280                 a[i__ + 1 + i__ * a_dim1] = e[i__];
 281             } else {
 282                 tauq[i__] = 0.;
 283             }
 284 /* L20: */
 285         }
 286     }
 287     return 0;
 288
 289 /*     End of DGEBD2 */
 290
 291 } /* dgebd2_ */