git.maemo.org Git - opencv/blob - 3rdparty/lapack/sorglq.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 sorglq_(integer *m, integer *n, integer *k, real *a,
  11         integer *lda, 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;
  15
  16     /* Local variables */
  17     integer i__, j, l, ib, nb, ki, kk, nx, iws, nbmin, iinfo;
  18     extern /* Subroutine */ int sorgl2_(integer *, integer *, integer *, real
  19             *, integer *, real *, real *, integer *), slarfb_(char *, char *,
  20             char *, char *, integer *, integer *, integer *, real *, integer *
  21 , real *, 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 /*  SORGLQ generates an M-by-N real matrix Q with orthonormal rows, */
  43 /*  which is defined as the first M rows of a product of K elementary */
  44 /*  reflectors of order N */
  45
  46 /*        Q  =  H(k) . . . H(2) H(1) */
  47
  48 /*  as returned by SGELQF. */
  49
  50 /*  Arguments */
  51 /*  ========= */
  52
  53 /*  M       (input) INTEGER */
  54 /*          The number of rows of the matrix Q. M >= 0. */
  55
  56 /*  N       (input) INTEGER */
  57 /*          The number of columns of the matrix Q. N >= M. */
  58
  59 /*  K       (input) INTEGER */
  60 /*          The number of elementary reflectors whose product defines the */
  61 /*          matrix Q. M >= K >= 0. */
  62
  63 /*  A       (input/output) REAL array, dimension (LDA,N) */
  64 /*          On entry, the i-th row must contain the vector which defines */
  65 /*          the elementary reflector H(i), for i = 1,2,...,k, as returned */
  66 /*          by SGELQF in the first k rows of its array argument A. */
  67 /*          On exit, the M-by-N matrix Q. */
  68
  69 /*  LDA     (input) INTEGER */
  70 /*          The first dimension of the array A. LDA >= max(1,M). */
  71
  72 /*  TAU     (input) REAL array, dimension (K) */
  73 /*          TAU(i) must contain the scalar factor of the elementary */
  74 /*          reflector H(i), as returned by SGELQF. */
  75
  76 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
  77 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  78
  79 /*  LWORK   (input) INTEGER */
  80 /*          The dimension of the array WORK. LWORK >= max(1,M). */
  81 /*          For optimum performance LWORK >= M*NB, where NB is */
  82 /*          the optimal blocksize. */
  83
  84 /*          If LWORK = -1, then a workspace query is assumed; the routine */
  85 /*          only calculates the optimal size of the WORK array, returns */
  86 /*          this value as the first entry of the WORK array, and no error */
  87 /*          message related to LWORK is issued by XERBLA. */
  88
  89 /*  INFO    (output) INTEGER */
  90 /*          = 0:  successful exit */
  91 /*          < 0:  if INFO = -i, the i-th argument has an illegal value */
  92
  93 /*  ===================================================================== */
  94
  95 /*     .. Parameters .. */
  96 /*     .. */
  97 /*     .. Local Scalars .. */
  98 /*     .. */
  99 /*     .. External Subroutines .. */
 100 /*     .. */
 101 /*     .. Intrinsic Functions .. */
 102 /*     .. */
 103 /*     .. External Functions .. */
 104 /*     .. */
 105 /*     .. Executable Statements .. */
 106
 107 /*     Test the input arguments */
 108
 109     /* Parameter adjustments */
 110     a_dim1 = *lda;
 111     a_offset = 1 + a_dim1;
 112     a -= a_offset;
 113     --tau;
 114     --work;
 115
 116     /* Function Body */
 117     *info = 0;
 118     nb = ilaenv_(&c__1, "SORGLQ", " ", m, n, k, &c_n1);
 119     lwkopt = max(1,*m) * nb;
 120     work[1] = (real) lwkopt;
 121     lquery = *lwork == -1;
 122     if (*m < 0) {
 123         *info = -1;
 124     } else if (*n < *m) {
 125         *info = -2;
 126     } else if (*k < 0 || *k > *m) {
 127         *info = -3;
 128     } else if (*lda < max(1,*m)) {
 129         *info = -5;
 130     } else if (*lwork < max(1,*m) && ! lquery) {
 131         *info = -8;
 132     }
 133     if (*info != 0) {
 134         i__1 = -(*info);
 135         xerbla_("SORGLQ", &i__1);
 136         return 0;
 137     } else if (lquery) {
 138         return 0;
 139     }
 140
 141 /*     Quick return if possible */
 142
 143     if (*m <= 0) {
 144         work[1] = 1.f;
 145         return 0;
 146     }
 147
 148     nbmin = 2;
 149     nx = 0;
 150     iws = *m;
 151     if (nb > 1 && nb < *k) {
 152
 153 /*        Determine when to cross over from blocked to unblocked code. */
 154
 155 /* Computing MAX */
 156         i__1 = 0, i__2 = ilaenv_(&c__3, "SORGLQ", " ", m, n, k, &c_n1);
 157         nx = max(i__1,i__2);
 158         if (nx < *k) {
 159
 160 /*           Determine if workspace is large enough for blocked code. */
 161
 162             ldwork = *m;
 163             iws = ldwork * nb;
 164             if (*lwork < iws) {
 165
 166 /*              Not enough workspace to use optimal NB:  reduce NB and */
 167 /*              determine the minimum value of NB. */
 168
 169                 nb = *lwork / ldwork;
 170 /* Computing MAX */
 171                 i__1 = 2, i__2 = ilaenv_(&c__2, "SORGLQ", " ", m, n, k, &c_n1);
 172                 nbmin = max(i__1,i__2);
 173             }
 174         }
 175     }
 176
 177     if (nb >= nbmin && nb < *k && nx < *k) {
 178
 179 /*        Use blocked code after the last block. */
 180 /*        The first kk rows are handled by the block method. */
 181
 182         ki = (*k - nx - 1) / nb * nb;
 183 /* Computing MIN */
 184         i__1 = *k, i__2 = ki + nb;
 185         kk = min(i__1,i__2);
 186
 187 /*        Set A(kk+1:m,1:kk) to zero. */
 188
 189         i__1 = kk;
 190         for (j = 1; j <= i__1; ++j) {
 191             i__2 = *m;
 192             for (i__ = kk + 1; i__ <= i__2; ++i__) {
 193                 a[i__ + j * a_dim1] = 0.f;
 194 /* L10: */
 195             }
 196 /* L20: */
 197         }
 198     } else {
 199         kk = 0;
 200     }
 201
 202 /*     Use unblocked code for the last or only block. */
 203
 204     if (kk < *m) {
 205         i__1 = *m - kk;
 206         i__2 = *n - kk;
 207         i__3 = *k - kk;
 208         sorgl2_(&i__1, &i__2, &i__3, &a[kk + 1 + (kk + 1) * a_dim1], lda, &
 209                 tau[kk + 1], &work[1], &iinfo);
 210     }
 211
 212     if (kk > 0) {
 213
 214 /*        Use blocked code */
 215
 216         i__1 = -nb;
 217         for (i__ = ki + 1; i__1 < 0 ? i__ >= 1 : i__ <= 1; i__ += i__1) {
 218 /* Computing MIN */
 219             i__2 = nb, i__3 = *k - i__ + 1;
 220             ib = min(i__2,i__3);
 221             if (i__ + ib <= *m) {
 222
 223 /*              Form the triangular factor of the block reflector */
 224 /*              H = H(i) H(i+1) . . . H(i+ib-1) */
 225
 226                 i__2 = *n - i__ + 1;
 227                 slarft_("Forward", "Rowwise", &i__2, &ib, &a[i__ + i__ *
 228                         a_dim1], lda, &tau[i__], &work[1], &ldwork);
 229
 230 /*              Apply H' to A(i+ib:m,i:n) from the right */
 231
 232                 i__2 = *m - i__ - ib + 1;
 233                 i__3 = *n - i__ + 1;
 234                 slarfb_("Right", "Transpose", "Forward", "Rowwise", &i__2, &
 235                         i__3, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
 236                         ldwork, &a[i__ + ib + i__ * a_dim1], lda, &work[ib +
 237                         1], &ldwork);
 238             }
 239
 240 /*           Apply H' to columns i:n of current block */
 241
 242             i__2 = *n - i__ + 1;
 243             sorgl2_(&ib, &i__2, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &
 244                     work[1], &iinfo);
 245
 246 /*           Set columns 1:i-1 of current block to zero */
 247
 248             i__2 = i__ - 1;
 249             for (j = 1; j <= i__2; ++j) {
 250                 i__3 = i__ + ib - 1;
 251                 for (l = i__; l <= i__3; ++l) {
 252                     a[l + j * a_dim1] = 0.f;
 253 /* L30: */
 254                 }
 255 /* L40: */
 256             }
 257 /* L50: */
 258         }
 259     }
 260
 261     work[1] = (real) iws;
 262     return 0;
 263
 264 /*     End of SORGLQ */
 265
 266 } /* sorglq_ */