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