git.maemo.org Git - opencv/blob - 3rdparty/lapack/dormqr.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 dormqr_(char *side, char *trans, integer *m, integer *n,
  11         integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *
  12         c__, integer *ldc, doublereal *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     doublereal t[4160]  /* was [65][64] */;
  26     integer i1, i2, i3, ib, ic, jc, nb, mi, ni, nq, nw, iws;
  27     logical left;
  28     extern logical lsame_(char *, char *);
  29     integer nbmin, iinfo;
  30     extern /* Subroutine */ int dorm2r_(char *, char *, integer *, integer *,
  31             integer *, doublereal *, integer *, doublereal *, doublereal *,
  32             integer *, doublereal *, integer *), dlarfb_(char
  33             *, char *, char *, char *, integer *, integer *, integer *,
  34             doublereal *, integer *, doublereal *, integer *, doublereal *,
  35             integer *, doublereal *, integer *), dlarft_(char *, char *, integer *, integer *, doublereal
  36             *, integer *, doublereal *, doublereal *, integer *), xerbla_(char *, integer *);
  37     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  38             integer *, 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 /*  DORMQR 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(1) H(2) . . . H(k) */
  66
  67 /*  as returned by DGEQRF. 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) DOUBLE PRECISION 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 /*          DGEQRF in the first 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) DOUBLE PRECISION array, dimension (K) */
 105 /*          TAU(i) must contain the scalar factor of the elementary */
 106 /*          reflector H(i), as returned by DGEQRF. */
 107
 108 /*  C       (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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 = *n;
 174     } else {
 175         nq = *n;
 176         nw = *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     } else if (*lwork < max(1,nw) && ! lquery) {
 193         *info = -12;
 194     }
 195
 196     if (*info == 0) {
 197
 198 /*        Determine the block size.  NB may be at most NBMAX, where NBMAX */
 199 /*        is used to define the local array T. */
 200
 201 /* Computing MIN */
 202 /* Writing concatenation */
 203         i__3[0] = 1, a__1[0] = side;
 204         i__3[1] = 1, a__1[1] = trans;
 205         s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 206         i__1 = 64, i__2 = ilaenv_(&c__1, "DORMQR", ch__1, m, n, k, &c_n1);
 207         nb = min(i__1,i__2);
 208         lwkopt = max(1,nw) * nb;
 209         work[1] = (doublereal) lwkopt;
 210     }
 211
 212     if (*info != 0) {
 213         i__1 = -(*info);
 214         xerbla_("DORMQR", &i__1);
 215         return 0;
 216     } else if (lquery) {
 217         return 0;
 218     }
 219
 220 /*     Quick return if possible */
 221
 222     if (*m == 0 || *n == 0 || *k == 0) {
 223         work[1] = 1.;
 224         return 0;
 225     }
 226
 227     nbmin = 2;
 228     ldwork = nw;
 229     if (nb > 1 && nb < *k) {
 230         iws = nw * nb;
 231         if (*lwork < iws) {
 232             nb = *lwork / ldwork;
 233 /* Computing MAX */
 234 /* Writing concatenation */
 235             i__3[0] = 1, a__1[0] = side;
 236             i__3[1] = 1, a__1[1] = trans;
 237             s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 238             i__1 = 2, i__2 = ilaenv_(&c__2, "DORMQR", ch__1, m, n, k, &c_n1);
 239             nbmin = max(i__1,i__2);
 240         }
 241     } else {
 242         iws = nw;
 243     }
 244
 245     if (nb < nbmin || nb >= *k) {
 246
 247 /*        Use unblocked code */
 248
 249         dorm2r_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 250                 c_offset], ldc, &work[1], &iinfo);
 251     } else {
 252
 253 /*        Use blocked code */
 254
 255         if (left && ! notran || ! left && notran) {
 256             i1 = 1;
 257             i2 = *k;
 258             i3 = nb;
 259         } else {
 260             i1 = (*k - 1) / nb * nb + 1;
 261             i2 = 1;
 262             i3 = -nb;
 263         }
 264
 265         if (left) {
 266             ni = *n;
 267             jc = 1;
 268         } else {
 269             mi = *m;
 270             ic = 1;
 271         }
 272
 273         i__1 = i2;
 274         i__2 = i3;
 275         for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 276 /* Computing MIN */
 277             i__4 = nb, i__5 = *k - i__ + 1;
 278             ib = min(i__4,i__5);
 279
 280 /*           Form the triangular factor of the block reflector */
 281 /*           H = H(i) H(i+1) . . . H(i+ib-1) */
 282
 283             i__4 = nq - i__ + 1;
 284             dlarft_("Forward", "Columnwise", &i__4, &ib, &a[i__ + i__ *
 285                     a_dim1], lda, &tau[i__], t, &c__65)
 286                     ;
 287             if (left) {
 288
 289 /*              H or H' is applied to C(i:m,1:n) */
 290
 291                 mi = *m - i__ + 1;
 292                 ic = i__;
 293             } else {
 294
 295 /*              H or H' is applied to C(1:m,i:n) */
 296
 297                 ni = *n - i__ + 1;
 298                 jc = i__;
 299             }
 300
 301 /*           Apply H or H' */
 302
 303             dlarfb_(side, trans, "Forward", "Columnwise", &mi, &ni, &ib, &a[
 304                     i__ + i__ * a_dim1], lda, t, &c__65, &c__[ic + jc *
 305                     c_dim1], ldc, &work[1], &ldwork);
 306 /* L10: */
 307         }
 308     }
 309     work[1] = (doublereal) lwkopt;
 310     return 0;
 311
 312 /*     End of DORMQR */
 313
 314 } /* dormqr_ */