git.maemo.org Git - opencv/blob - 3rdparty/lapack/sorgbr.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 sorgbr_(char *vect, integer *m, integer *n, integer *k,
   9         real *a, integer *lda, real *tau, real *work, integer *lwork, integer
  10         *info)
  11 {
  12     /* System generated locals */
  13     integer a_dim1, a_offset, i__1, i__2, i__3;
  14
  15     /* Local variables */
  16     integer i__, j, nb, mn;
  17     extern logical lsame_(char *, char *);
  18     integer iinfo;
  19     logical wantq;
  20     extern /* Subroutine */ int xerbla_(char *, integer *);
  21     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  22             integer *, integer *);
  23     extern /* Subroutine */ int sorglq_(integer *, integer *, integer *, real
  24             *, integer *, real *, real *, integer *, integer *), sorgqr_(
  25             integer *, integer *, integer *, real *, integer *, real *, real *
  26 , integer *, integer *);
  27     integer lwkopt;
  28     logical lquery;
  29
  30
  31 /*  -- LAPACK routine (version 3.1) -- */
  32 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  33 /*     November 2006 */
  34
  35 /*     .. Scalar Arguments .. */
  36 /*     .. */
  37 /*     .. Array Arguments .. */
  38 /*     .. */
  39
  40 /*  Purpose */
  41 /*  ======= */
  42
  43 /*  SORGBR generates one of the real orthogonal matrices Q or P**T */
  44 /*  determined by SGEBRD when reducing a real matrix A to bidiagonal */
  45 /*  form: A = Q * B * P**T.  Q and P**T are defined as products of */
  46 /*  elementary reflectors H(i) or G(i) respectively. */
  47
  48 /*  If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q */
  49 /*  is of order M: */
  50 /*  if m >= k, Q = H(1) H(2) . . . H(k) and SORGBR returns the first n */
  51 /*  columns of Q, where m >= n >= k; */
  52 /*  if m < k, Q = H(1) H(2) . . . H(m-1) and SORGBR returns Q as an */
  53 /*  M-by-M matrix. */
  54
  55 /*  If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T */
  56 /*  is of order N: */
  57 /*  if k < n, P**T = G(k) . . . G(2) G(1) and SORGBR returns the first m */
  58 /*  rows of P**T, where n >= m >= k; */
  59 /*  if k >= n, P**T = G(n-1) . . . G(2) G(1) and SORGBR returns P**T as */
  60 /*  an N-by-N matrix. */
  61
  62 /*  Arguments */
  63 /*  ========= */
  64
  65 /*  VECT    (input) CHARACTER*1 */
  66 /*          Specifies whether the matrix Q or the matrix P**T is */
  67 /*          required, as defined in the transformation applied by SGEBRD: */
  68 /*          = 'Q':  generate Q; */
  69 /*          = 'P':  generate P**T. */
  70
  71 /*  M       (input) INTEGER */
  72 /*          The number of rows of the matrix Q or P**T to be returned. */
  73 /*          M >= 0. */
  74
  75 /*  N       (input) INTEGER */
  76 /*          The number of columns of the matrix Q or P**T to be returned. */
  77 /*          N >= 0. */
  78 /*          If VECT = 'Q', M >= N >= min(M,K); */
  79 /*          if VECT = 'P', N >= M >= min(N,K). */
  80
  81 /*  K       (input) INTEGER */
  82 /*          If VECT = 'Q', the number of columns in the original M-by-K */
  83 /*          matrix reduced by SGEBRD. */
  84 /*          If VECT = 'P', the number of rows in the original K-by-N */
  85 /*          matrix reduced by SGEBRD. */
  86 /*          K >= 0. */
  87
  88 /*  A       (input/output) REAL array, dimension (LDA,N) */
  89 /*          On entry, the vectors which define the elementary reflectors, */
  90 /*          as returned by SGEBRD. */
  91 /*          On exit, the M-by-N matrix Q or P**T. */
  92
  93 /*  LDA     (input) INTEGER */
  94 /*          The leading dimension of the array A. LDA >= max(1,M). */
  95
  96 /*  TAU     (input) REAL array, dimension */
  97 /*                                (min(M,K)) if VECT = 'Q' */
  98 /*                                (min(N,K)) if VECT = 'P' */
  99 /*          TAU(i) must contain the scalar factor of the elementary */
 100 /*          reflector H(i) or G(i), which determines Q or P**T, as */
 101 /*          returned by SGEBRD in its array argument TAUQ or TAUP. */
 102
 103 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
 104 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 105
 106 /*  LWORK   (input) INTEGER */
 107 /*          The dimension of the array WORK. LWORK >= max(1,min(M,N)). */
 108 /*          For optimum performance LWORK >= min(M,N)*NB, where NB */
 109 /*          is the optimal blocksize. */
 110
 111 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 112 /*          only calculates the optimal size of the WORK array, returns */
 113 /*          this value as the first entry of the WORK array, and no error */
 114 /*          message related to LWORK is issued by XERBLA. */
 115
 116 /*  INFO    (output) INTEGER */
 117 /*          = 0:  successful exit */
 118 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 119
 120 /*  ===================================================================== */
 121
 122 /*     .. Parameters .. */
 123 /*     .. */
 124 /*     .. Local Scalars .. */
 125 /*     .. */
 126 /*     .. External Functions .. */
 127 /*     .. */
 128 /*     .. External Subroutines .. */
 129 /*     .. */
 130 /*     .. Intrinsic Functions .. */
 131 /*     .. */
 132 /*     .. Executable Statements .. */
 133
 134 /*     Test the input arguments */
 135
 136     /* Parameter adjustments */
 137     a_dim1 = *lda;
 138     a_offset = 1 + a_dim1;
 139     a -= a_offset;
 140     --tau;
 141     --work;
 142
 143     /* Function Body */
 144     *info = 0;
 145     wantq = lsame_(vect, "Q");
 146     mn = min(*m,*n);
 147     lquery = *lwork == -1;
 148     if (! wantq && ! lsame_(vect, "P")) {
 149         *info = -1;
 150     } else if (*m < 0) {
 151         *info = -2;
 152     } else if (*n < 0 || wantq && (*n > *m || *n < min(*m,*k)) || ! wantq && (
 153             *m > *n || *m < min(*n,*k))) {
 154         *info = -3;
 155     } else if (*k < 0) {
 156         *info = -4;
 157     } else if (*lda < max(1,*m)) {
 158         *info = -6;
 159     } else if (*lwork < max(1,mn) && ! lquery) {
 160         *info = -9;
 161     }
 162
 163     if (*info == 0) {
 164         if (wantq) {
 165             nb = ilaenv_(&c__1, "SORGQR", " ", m, n, k, &c_n1);
 166         } else {
 167             nb = ilaenv_(&c__1, "SORGLQ", " ", m, n, k, &c_n1);
 168         }
 169         lwkopt = max(1,mn) * nb;
 170         work[1] = (real) lwkopt;
 171     }
 172
 173     if (*info != 0) {
 174         i__1 = -(*info);
 175         xerbla_("SORGBR", &i__1);
 176         return 0;
 177     } else if (lquery) {
 178         return 0;
 179     }
 180
 181 /*     Quick return if possible */
 182
 183     if (*m == 0 || *n == 0) {
 184         work[1] = 1.f;
 185         return 0;
 186     }
 187
 188     if (wantq) {
 189
 190 /*        Form Q, determined by a call to SGEBRD to reduce an m-by-k */
 191 /*        matrix */
 192
 193         if (*m >= *k) {
 194
 195 /*           If m >= k, assume m >= n >= k */
 196
 197             sorgqr_(m, n, k, &a[a_offset], lda, &tau[1], &work[1], lwork, &
 198                     iinfo);
 199
 200         } else {
 201
 202 /*           If m < k, assume m = n */
 203
 204 /*           Shift the vectors which define the elementary reflectors one */
 205 /*           column to the right, and set the first row and column of Q */
 206 /*           to those of the unit matrix */
 207
 208             for (j = *m; j >= 2; --j) {
 209                 a[j * a_dim1 + 1] = 0.f;
 210                 i__1 = *m;
 211                 for (i__ = j + 1; i__ <= i__1; ++i__) {
 212                     a[i__ + j * a_dim1] = a[i__ + (j - 1) * a_dim1];
 213 /* L10: */
 214                 }
 215 /* L20: */
 216             }
 217             a[a_dim1 + 1] = 1.f;
 218             i__1 = *m;
 219             for (i__ = 2; i__ <= i__1; ++i__) {
 220                 a[i__ + a_dim1] = 0.f;
 221 /* L30: */
 222             }
 223             if (*m > 1) {
 224
 225 /*              Form Q(2:m,2:m) */
 226
 227                 i__1 = *m - 1;
 228                 i__2 = *m - 1;
 229                 i__3 = *m - 1;
 230                 sorgqr_(&i__1, &i__2, &i__3, &a[(a_dim1 << 1) + 2], lda, &tau[
 231                         1], &work[1], lwork, &iinfo);
 232             }
 233         }
 234     } else {
 235
 236 /*        Form P', determined by a call to SGEBRD to reduce a k-by-n */
 237 /*        matrix */
 238
 239         if (*k < *n) {
 240
 241 /*           If k < n, assume k <= m <= n */
 242
 243             sorglq_(m, n, k, &a[a_offset], lda, &tau[1], &work[1], lwork, &
 244                     iinfo);
 245
 246         } else {
 247
 248 /*           If k >= n, assume m = n */
 249
 250 /*           Shift the vectors which define the elementary reflectors one */
 251 /*           row downward, and set the first row and column of P' to */
 252 /*           those of the unit matrix */
 253
 254             a[a_dim1 + 1] = 1.f;
 255             i__1 = *n;
 256             for (i__ = 2; i__ <= i__1; ++i__) {
 257                 a[i__ + a_dim1] = 0.f;
 258 /* L40: */
 259             }
 260             i__1 = *n;
 261             for (j = 2; j <= i__1; ++j) {
 262                 for (i__ = j - 1; i__ >= 2; --i__) {
 263                     a[i__ + j * a_dim1] = a[i__ - 1 + j * a_dim1];
 264 /* L50: */
 265                 }
 266                 a[j * a_dim1 + 1] = 0.f;
 267 /* L60: */
 268             }
 269             if (*n > 1) {
 270
 271 /*              Form P'(2:n,2:n) */
 272
 273                 i__1 = *n - 1;
 274                 i__2 = *n - 1;
 275                 i__3 = *n - 1;
 276                 sorglq_(&i__1, &i__2, &i__3, &a[(a_dim1 << 1) + 2], lda, &tau[
 277                         1], &work[1], lwork, &iinfo);
 278             }
 279         }
 280     }
 281     work[1] = (real) lwkopt;
 282     return 0;
 283
 284 /*     End of SORGBR */
 285
 286 } /* sorgbr_ */