git.maemo.org Git - opencv/blob - 3rdparty/lapack/sgeqrf.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 static integer c__3 = 3;
   8 static integer c__2 = 2;
   9
  10 /* Subroutine */ int sgeqrf_(integer *m, integer *n, real *a, integer *lda,
  11         real *tau, real *work, integer *lwork, integer *info)
  12 {
  13     /* System generated locals */
  14     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
  15
  16     /* Local variables */
  17     integer i__, k, ib, nb, nx, iws, nbmin, iinfo;
  18     extern /* Subroutine */ int sgeqr2_(integer *, integer *, real *, integer
  19             *, real *, real *, integer *), slarfb_(char *, char *, char *,
  20             char *, integer *, integer *, integer *, real *, integer *, real *
  21 , integer *, real *, integer *, real *, integer *), xerbla_(char *, integer *);
  22     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  23             integer *, integer *);
  24     extern /* Subroutine */ int slarft_(char *, char *, integer *, integer *,
  25             real *, integer *, real *, real *, integer *);
  26     integer ldwork, lwkopt;
  27     logical lquery;
  28
  29
  30 /*  -- LAPACK routine (version 3.1) -- */
  31 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  32 /*     November 2006 */
  33
  34 /*     .. Scalar Arguments .. */
  35 /*     .. */
  36 /*     .. Array Arguments .. */
  37 /*     .. */
  38
  39 /*  Purpose */
  40 /*  ======= */
  41
  42 /*  SGEQRF computes a QR factorization of a real M-by-N matrix A: */
  43 /*  A = Q * R. */
  44
  45 /*  Arguments */
  46 /*  ========= */
  47
  48 /*  M       (input) INTEGER */
  49 /*          The number of rows of the matrix A.  M >= 0. */
  50
  51 /*  N       (input) INTEGER */
  52 /*          The number of columns of the matrix A.  N >= 0. */
  53
  54 /*  A       (input/output) REAL array, dimension (LDA,N) */
  55 /*          On entry, the M-by-N matrix A. */
  56 /*          On exit, the elements on and above the diagonal of the array */
  57 /*          contain the min(M,N)-by-N upper trapezoidal matrix R (R is */
  58 /*          upper triangular if m >= n); the elements below the diagonal, */
  59 /*          with the array TAU, represent the orthogonal matrix Q as a */
  60 /*          product of min(m,n) elementary reflectors (see Further */
  61 /*          Details). */
  62
  63 /*  LDA     (input) INTEGER */
  64 /*          The leading dimension of the array A.  LDA >= max(1,M). */
  65
  66 /*  TAU     (output) REAL array, dimension (min(M,N)) */
  67 /*          The scalar factors of the elementary reflectors (see Further */
  68 /*          Details). */
  69
  70 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
  71 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  72
  73 /*  LWORK   (input) INTEGER */
  74 /*          The dimension of the array WORK.  LWORK >= max(1,N). */
  75 /*          For optimum performance LWORK >= N*NB, where NB is */
  76 /*          the optimal blocksize. */
  77
  78 /*          If LWORK = -1, then a workspace query is assumed; the routine */
  79 /*          only calculates the optimal size of the WORK array, returns */
  80 /*          this value as the first entry of the WORK array, and no error */
  81 /*          message related to LWORK is issued by XERBLA. */
  82
  83 /*  INFO    (output) INTEGER */
  84 /*          = 0:  successful exit */
  85 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
  86
  87 /*  Further Details */
  88 /*  =============== */
  89
  90 /*  The matrix Q is represented as a product of elementary reflectors */
  91
  92 /*     Q = H(1) H(2) . . . H(k), where k = min(m,n). */
  93
  94 /*  Each H(i) has the form */
  95
  96 /*     H(i) = I - tau * v * v' */
  97
  98 /*  where tau is a real scalar, and v is a real vector with */
  99 /*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), */
 100 /*  and tau in TAU(i). */
 101
 102 /*  ===================================================================== */
 103
 104 /*     .. Local Scalars .. */
 105 /*     .. */
 106 /*     .. External Subroutines .. */
 107 /*     .. */
 108 /*     .. Intrinsic Functions .. */
 109 /*     .. */
 110 /*     .. External Functions .. */
 111 /*     .. */
 112 /*     .. Executable Statements .. */
 113
 114 /*     Test the input arguments */
 115
 116     /* Parameter adjustments */
 117     a_dim1 = *lda;
 118     a_offset = 1 + a_dim1;
 119     a -= a_offset;
 120     --tau;
 121     --work;
 122
 123     /* Function Body */
 124     *info = 0;
 125     nb = ilaenv_(&c__1, "SGEQRF", " ", m, n, &c_n1, &c_n1);
 126     lwkopt = *n * nb;
 127     work[1] = (real) lwkopt;
 128     lquery = *lwork == -1;
 129     if (*m < 0) {
 130         *info = -1;
 131     } else if (*n < 0) {
 132         *info = -2;
 133     } else if (*lda < max(1,*m)) {
 134         *info = -4;
 135     } else if (*lwork < max(1,*n) && ! lquery) {
 136         *info = -7;
 137     }
 138     if (*info != 0) {
 139         i__1 = -(*info);
 140         xerbla_("SGEQRF", &i__1);
 141         return 0;
 142     } else if (lquery) {
 143         return 0;
 144     }
 145
 146 /*     Quick return if possible */
 147
 148     k = min(*m,*n);
 149     if (k == 0) {
 150         work[1] = 1.f;
 151         return 0;
 152     }
 153
 154     nbmin = 2;
 155     nx = 0;
 156     iws = *n;
 157     if (nb > 1 && nb < k) {
 158
 159 /*        Determine when to cross over from blocked to unblocked code. */
 160
 161 /* Computing MAX */
 162         i__1 = 0, i__2 = ilaenv_(&c__3, "SGEQRF", " ", m, n, &c_n1, &c_n1);
 163         nx = max(i__1,i__2);
 164         if (nx < k) {
 165
 166 /*           Determine if workspace is large enough for blocked code. */
 167
 168             ldwork = *n;
 169             iws = ldwork * nb;
 170             if (*lwork < iws) {
 171
 172 /*              Not enough workspace to use optimal NB:  reduce NB and */
 173 /*              determine the minimum value of NB. */
 174
 175                 nb = *lwork / ldwork;
 176 /* Computing MAX */
 177                 i__1 = 2, i__2 = ilaenv_(&c__2, "SGEQRF", " ", m, n, &c_n1, &
 178                         c_n1);
 179                 nbmin = max(i__1,i__2);
 180             }
 181         }
 182     }
 183
 184     if (nb >= nbmin && nb < k && nx < k) {
 185
 186 /*        Use blocked code initially */
 187
 188         i__1 = k - nx;
 189         i__2 = nb;
 190         for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 191 /* Computing MIN */
 192             i__3 = k - i__ + 1;
 193             ib = min(i__3,nb);
 194
 195 /*           Compute the QR factorization of the current block */
 196 /*           A(i:m,i:i+ib-1) */
 197
 198             i__3 = *m - i__ + 1;
 199             sgeqr2_(&i__3, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
 200                     1], &iinfo);
 201             if (i__ + ib <= *n) {
 202
 203 /*              Form the triangular factor of the block reflector */
 204 /*              H = H(i) H(i+1) . . . H(i+ib-1) */
 205
 206                 i__3 = *m - i__ + 1;
 207                 slarft_("Forward", "Columnwise", &i__3, &ib, &a[i__ + i__ *
 208                         a_dim1], lda, &tau[i__], &work[1], &ldwork);
 209
 210 /*              Apply H' to A(i:m,i+ib:n) from the left */
 211
 212                 i__3 = *m - i__ + 1;
 213                 i__4 = *n - i__ - ib + 1;
 214                 slarfb_("Left", "Transpose", "Forward", "Columnwise", &i__3, &
 215                         i__4, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
 216                         ldwork, &a[i__ + (i__ + ib) * a_dim1], lda, &work[ib
 217                         + 1], &ldwork);
 218             }
 219 /* L10: */
 220         }
 221     } else {
 222         i__ = 1;
 223     }
 224
 225 /*     Use unblocked code to factor the last or only block. */
 226
 227     if (i__ <= k) {
 228         i__2 = *m - i__ + 1;
 229         i__1 = *n - i__ + 1;
 230         sgeqr2_(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
 231 , &iinfo);
 232     }
 233
 234     work[1] = (real) iws;
 235     return 0;
 236
 237 /*     End of SGEQRF */
 238
 239 } /* sgeqrf_ */