git.maemo.org Git - opencv/blob - 3rdparty/lapack/sormql.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__2 = 2;
   8 static integer c__65 = 65;
   9
  10 /* Subroutine */ int sormql_(char *side, char *trans, integer *m, integer *n,
  11         integer *k, real *a, integer *lda, real *tau, real *c__, integer *ldc,
  12          real *work, integer *lwork, integer *info)
  13 {
  14     /* System generated locals */
  15     address a__1[2];
  16     integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3[2], i__4,
  17             i__5;
  18     char ch__1[2];
  19
  20     /* Builtin functions */
  21     /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
  22
  23     /* Local variables */
  24     integer i__;
  25     real t[4160]        /* was [65][64] */;
  26     integer i1, i2, i3, ib, nb, mi, ni, nq, nw, iws;
  27     logical left;
  28     extern logical lsame_(char *, char *);
  29     integer nbmin, iinfo;
  30     extern /* Subroutine */ int sorm2l_(char *, char *, integer *, integer *,
  31             integer *, real *, integer *, real *, real *, integer *, real *,
  32             integer *), slarfb_(char *, char *, char *, char *
  33 , integer *, integer *, integer *, real *, integer *, real *,
  34             integer *, real *, integer *, real *, integer *), xerbla_(char *, integer *);
  35     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  36             integer *, integer *);
  37     extern /* Subroutine */ int slarft_(char *, char *, integer *, integer *,
  38             real *, integer *, real *, real *, integer *);
  39     logical notran;
  40     integer ldwork, lwkopt;
  41     logical lquery;
  42
  43
  44 /*  -- LAPACK routine (version 3.1) -- */
  45 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  46 /*     November 2006 */
  47
  48 /*     .. Scalar Arguments .. */
  49 /*     .. */
  50 /*     .. Array Arguments .. */
  51 /*     .. */
  52
  53 /*  Purpose */
  54 /*  ======= */
  55
  56 /*  SORMQL overwrites the general real M-by-N matrix C with */
  57
  58 /*                  SIDE = 'L'     SIDE = 'R' */
  59 /*  TRANS = 'N':      Q * C          C * Q */
  60 /*  TRANS = 'T':      Q**T * C       C * Q**T */
  61
  62 /*  where Q is a real orthogonal matrix defined as the product of k */
  63 /*  elementary reflectors */
  64
  65 /*        Q = H(k) . . . H(2) H(1) */
  66
  67 /*  as returned by SGEQLF. Q is of order M if SIDE = 'L' and of order N */
  68 /*  if SIDE = 'R'. */
  69
  70 /*  Arguments */
  71 /*  ========= */
  72
  73 /*  SIDE    (input) CHARACTER*1 */
  74 /*          = 'L': apply Q or Q**T from the Left; */
  75 /*          = 'R': apply Q or Q**T from the Right. */
  76
  77 /*  TRANS   (input) CHARACTER*1 */
  78 /*          = 'N':  No transpose, apply Q; */
  79 /*          = 'T':  Transpose, apply Q**T. */
  80
  81 /*  M       (input) INTEGER */
  82 /*          The number of rows of the matrix C. M >= 0. */
  83
  84 /*  N       (input) INTEGER */
  85 /*          The number of columns of the matrix C. N >= 0. */
  86
  87 /*  K       (input) INTEGER */
  88 /*          The number of elementary reflectors whose product defines */
  89 /*          the matrix Q. */
  90 /*          If SIDE = 'L', M >= K >= 0; */
  91 /*          if SIDE = 'R', N >= K >= 0. */
  92
  93 /*  A       (input) REAL array, dimension (LDA,K) */
  94 /*          The i-th column must contain the vector which defines the */
  95 /*          elementary reflector H(i), for i = 1,2,...,k, as returned by */
  96 /*          SGEQLF in the last k columns of its array argument A. */
  97 /*          A is modified by the routine but restored on exit. */
  98
  99 /*  LDA     (input) INTEGER */
 100 /*          The leading dimension of the array A. */
 101 /*          If SIDE = 'L', LDA >= max(1,M); */
 102 /*          if SIDE = 'R', LDA >= max(1,N). */
 103
 104 /*  TAU     (input) REAL array, dimension (K) */
 105 /*          TAU(i) must contain the scalar factor of the elementary */
 106 /*          reflector H(i), as returned by SGEQLF. */
 107
 108 /*  C       (input/output) REAL array, dimension (LDC,N) */
 109 /*          On entry, the M-by-N matrix C. */
 110 /*          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. */
 111
 112 /*  LDC     (input) INTEGER */
 113 /*          The leading dimension of the array C. LDC >= max(1,M). */
 114
 115 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
 116 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 117
 118 /*  LWORK   (input) INTEGER */
 119 /*          The dimension of the array WORK. */
 120 /*          If SIDE = 'L', LWORK >= max(1,N); */
 121 /*          if SIDE = 'R', LWORK >= max(1,M). */
 122 /*          For optimum performance LWORK >= N*NB if SIDE = 'L', and */
 123 /*          LWORK >= M*NB if SIDE = 'R', where NB is the optimal */
 124 /*          blocksize. */
 125
 126 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 127 /*          only calculates the optimal size of the WORK array, returns */
 128 /*          this value as the first entry of the WORK array, and no error */
 129 /*          message related to LWORK is issued by XERBLA. */
 130
 131 /*  INFO    (output) INTEGER */
 132 /*          = 0:  successful exit */
 133 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 134
 135 /*  ===================================================================== */
 136
 137 /*     .. Parameters .. */
 138 /*     .. */
 139 /*     .. Local Scalars .. */
 140 /*     .. */
 141 /*     .. Local Arrays .. */
 142 /*     .. */
 143 /*     .. External Functions .. */
 144 /*     .. */
 145 /*     .. External Subroutines .. */
 146 /*     .. */
 147 /*     .. Intrinsic Functions .. */
 148 /*     .. */
 149 /*     .. Executable Statements .. */
 150
 151 /*     Test the input arguments */
 152
 153     /* Parameter adjustments */
 154     a_dim1 = *lda;
 155     a_offset = 1 + a_dim1;
 156     a -= a_offset;
 157     --tau;
 158     c_dim1 = *ldc;
 159     c_offset = 1 + c_dim1;
 160     c__ -= c_offset;
 161     --work;
 162
 163     /* Function Body */
 164     *info = 0;
 165     left = lsame_(side, "L");
 166     notran = lsame_(trans, "N");
 167     lquery = *lwork == -1;
 168
 169 /*     NQ is the order of Q and NW is the minimum dimension of WORK */
 170
 171     if (left) {
 172         nq = *m;
 173         nw = max(1,*n);
 174     } else {
 175         nq = *n;
 176         nw = max(1,*m);
 177     }
 178     if (! left && ! lsame_(side, "R")) {
 179         *info = -1;
 180     } else if (! notran && ! lsame_(trans, "T")) {
 181         *info = -2;
 182     } else if (*m < 0) {
 183         *info = -3;
 184     } else if (*n < 0) {
 185         *info = -4;
 186     } else if (*k < 0 || *k > nq) {
 187         *info = -5;
 188     } else if (*lda < max(1,nq)) {
 189         *info = -7;
 190     } else if (*ldc < max(1,*m)) {
 191         *info = -10;
 192     }
 193
 194     if (*info == 0) {
 195         if (*m == 0 || *n == 0) {
 196             lwkopt = 1;
 197         } else {
 198
 199 /*           Determine the block size.  NB may be at most NBMAX, where */
 200 /*           NBMAX is used to define the local array T. */
 201
 202
 203 /* Computing MIN */
 204 /* Writing concatenation */
 205             i__3[0] = 1, a__1[0] = side;
 206             i__3[1] = 1, a__1[1] = trans;
 207             s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 208             i__1 = 64, i__2 = ilaenv_(&c__1, "SORMQL", ch__1, m, n, k, &c_n1);
 209             nb = min(i__1,i__2);
 210             lwkopt = nw * nb;
 211         }
 212         work[1] = (real) lwkopt;
 213
 214         if (*lwork < nw && ! lquery) {
 215             *info = -12;
 216         }
 217     }
 218
 219     if (*info != 0) {
 220         i__1 = -(*info);
 221         xerbla_("SORMQL", &i__1);
 222         return 0;
 223     } else if (lquery) {
 224         return 0;
 225     }
 226
 227 /*     Quick return if possible */
 228
 229     if (*m == 0 || *n == 0) {
 230         return 0;
 231     }
 232
 233     nbmin = 2;
 234     ldwork = nw;
 235     if (nb > 1 && nb < *k) {
 236         iws = nw * nb;
 237         if (*lwork < iws) {
 238             nb = *lwork / ldwork;
 239 /* Computing MAX */
 240 /* Writing concatenation */
 241             i__3[0] = 1, a__1[0] = side;
 242             i__3[1] = 1, a__1[1] = trans;
 243             s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 244             i__1 = 2, i__2 = ilaenv_(&c__2, "SORMQL", ch__1, m, n, k, &c_n1);
 245             nbmin = max(i__1,i__2);
 246         }
 247     } else {
 248         iws = nw;
 249     }
 250
 251     if (nb < nbmin || nb >= *k) {
 252
 253 /*        Use unblocked code */
 254
 255         sorm2l_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 256                 c_offset], ldc, &work[1], &iinfo);
 257     } else {
 258
 259 /*        Use blocked code */
 260
 261         if (left && notran || ! left && ! notran) {
 262             i1 = 1;
 263             i2 = *k;
 264             i3 = nb;
 265         } else {
 266             i1 = (*k - 1) / nb * nb + 1;
 267             i2 = 1;
 268             i3 = -nb;
 269         }
 270
 271         if (left) {
 272             ni = *n;
 273         } else {
 274             mi = *m;
 275         }
 276
 277         i__1 = i2;
 278         i__2 = i3;
 279         for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 280 /* Computing MIN */
 281             i__4 = nb, i__5 = *k - i__ + 1;
 282             ib = min(i__4,i__5);
 283
 284 /*           Form the triangular factor of the block reflector */
 285 /*           H = H(i+ib-1) . . . H(i+1) H(i) */
 286
 287             i__4 = nq - *k + i__ + ib - 1;
 288             slarft_("Backward", "Columnwise", &i__4, &ib, &a[i__ * a_dim1 + 1]
 289 , lda, &tau[i__], t, &c__65);
 290             if (left) {
 291
 292 /*              H or H' is applied to C(1:m-k+i+ib-1,1:n) */
 293
 294                 mi = *m - *k + i__ + ib - 1;
 295             } else {
 296
 297 /*              H or H' is applied to C(1:m,1:n-k+i+ib-1) */
 298
 299                 ni = *n - *k + i__ + ib - 1;
 300             }
 301
 302 /*           Apply H or H' */
 303
 304             slarfb_(side, trans, "Backward", "Columnwise", &mi, &ni, &ib, &a[
 305                     i__ * a_dim1 + 1], lda, t, &c__65, &c__[c_offset], ldc, &
 306                     work[1], &ldwork);
 307 /* L10: */
 308         }
 309     }
 310     work[1] = (real) lwkopt;
 311     return 0;
 312
 313 /*     End of SORMQL */
 314
 315 } /* sormql_ */