git.maemo.org Git - opencv/blob - 3rdparty/lapack/dormbr.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
   9 /* Subroutine */ int dormbr_(char *vect, char *side, char *trans, integer *m,
  10         integer *n, integer *k, doublereal *a, integer *lda, doublereal *tau,
  11         doublereal *c__, integer *ldc, doublereal *work, integer *lwork,
  12         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];
  17     char ch__1[2];
  18
  19     /* Builtin functions */
  20     /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
  21
  22     /* Local variables */
  23     integer i1, i2, nb, mi, ni, nq, nw;
  24     logical left;
  25     extern logical lsame_(char *, char *);
  26     integer iinfo;
  27     extern /* Subroutine */ int xerbla_(char *, integer *);
  28     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  29             integer *, integer *);
  30     extern /* Subroutine */ int dormlq_(char *, char *, integer *, integer *,
  31             integer *, doublereal *, integer *, doublereal *, doublereal *,
  32             integer *, doublereal *, integer *, integer *);
  33     logical notran;
  34     extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *,
  35             integer *, doublereal *, integer *, doublereal *, doublereal *,
  36             integer *, doublereal *, integer *, integer *);
  37     logical applyq;
  38     char transt[1];
  39     integer lwkopt;
  40     logical lquery;
  41
  42
  43 /*  -- LAPACK routine (version 3.1) -- */
  44 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  45 /*     November 2006 */
  46
  47 /*     .. Scalar Arguments .. */
  48 /*     .. */
  49 /*     .. Array Arguments .. */
  50 /*     .. */
  51
  52 /*  Purpose */
  53 /*  ======= */
  54
  55 /*  If VECT = 'Q', DORMBR overwrites the general real M-by-N matrix C */
  56 /*  with */
  57 /*                  SIDE = 'L'     SIDE = 'R' */
  58 /*  TRANS = 'N':      Q * C          C * Q */
  59 /*  TRANS = 'T':      Q**T * C       C * Q**T */
  60
  61 /*  If VECT = 'P', DORMBR overwrites the general real M-by-N matrix C */
  62 /*  with */
  63 /*                  SIDE = 'L'     SIDE = 'R' */
  64 /*  TRANS = 'N':      P * C          C * P */
  65 /*  TRANS = 'T':      P**T * C       C * P**T */
  66
  67 /*  Here Q and P**T are the orthogonal matrices determined by DGEBRD when */
  68 /*  reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and */
  69 /*  P**T are defined as products of elementary reflectors H(i) and G(i) */
  70 /*  respectively. */
  71
  72 /*  Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the */
  73 /*  order of the orthogonal matrix Q or P**T that is applied. */
  74
  75 /*  If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: */
  76 /*  if nq >= k, Q = H(1) H(2) . . . H(k); */
  77 /*  if nq < k, Q = H(1) H(2) . . . H(nq-1). */
  78
  79 /*  If VECT = 'P', A is assumed to have been a K-by-NQ matrix: */
  80 /*  if k < nq, P = G(1) G(2) . . . G(k); */
  81 /*  if k >= nq, P = G(1) G(2) . . . G(nq-1). */
  82
  83 /*  Arguments */
  84 /*  ========= */
  85
  86 /*  VECT    (input) CHARACTER*1 */
  87 /*          = 'Q': apply Q or Q**T; */
  88 /*          = 'P': apply P or P**T. */
  89
  90 /*  SIDE    (input) CHARACTER*1 */
  91 /*          = 'L': apply Q, Q**T, P or P**T from the Left; */
  92 /*          = 'R': apply Q, Q**T, P or P**T from the Right. */
  93
  94 /*  TRANS   (input) CHARACTER*1 */
  95 /*          = 'N':  No transpose, apply Q  or P; */
  96 /*          = 'T':  Transpose, apply Q**T or P**T. */
  97
  98 /*  M       (input) INTEGER */
  99 /*          The number of rows of the matrix C. M >= 0. */
 100
 101 /*  N       (input) INTEGER */
 102 /*          The number of columns of the matrix C. N >= 0. */
 103
 104 /*  K       (input) INTEGER */
 105 /*          If VECT = 'Q', the number of columns in the original */
 106 /*          matrix reduced by DGEBRD. */
 107 /*          If VECT = 'P', the number of rows in the original */
 108 /*          matrix reduced by DGEBRD. */
 109 /*          K >= 0. */
 110
 111 /*  A       (input) DOUBLE PRECISION array, dimension */
 112 /*                                (LDA,min(nq,K)) if VECT = 'Q' */
 113 /*                                (LDA,nq)        if VECT = 'P' */
 114 /*          The vectors which define the elementary reflectors H(i) and */
 115 /*          G(i), whose products determine the matrices Q and P, as */
 116 /*          returned by DGEBRD. */
 117
 118 /*  LDA     (input) INTEGER */
 119 /*          The leading dimension of the array A. */
 120 /*          If VECT = 'Q', LDA >= max(1,nq); */
 121 /*          if VECT = 'P', LDA >= max(1,min(nq,K)). */
 122
 123 /*  TAU     (input) DOUBLE PRECISION array, dimension (min(nq,K)) */
 124 /*          TAU(i) must contain the scalar factor of the elementary */
 125 /*          reflector H(i) or G(i) which determines Q or P, as returned */
 126 /*          by DGEBRD in the array argument TAUQ or TAUP. */
 127
 128 /*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
 129 /*          On entry, the M-by-N matrix C. */
 130 /*          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q */
 131 /*          or P*C or P**T*C or C*P or C*P**T. */
 132
 133 /*  LDC     (input) INTEGER */
 134 /*          The leading dimension of the array C. LDC >= max(1,M). */
 135
 136 /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
 137 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 138
 139 /*  LWORK   (input) INTEGER */
 140 /*          The dimension of the array WORK. */
 141 /*          If SIDE = 'L', LWORK >= max(1,N); */
 142 /*          if SIDE = 'R', LWORK >= max(1,M). */
 143 /*          For optimum performance LWORK >= N*NB if SIDE = 'L', and */
 144 /*          LWORK >= M*NB if SIDE = 'R', where NB is the optimal */
 145 /*          blocksize. */
 146
 147 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 148 /*          only calculates the optimal size of the WORK array, returns */
 149 /*          this value as the first entry of the WORK array, and no error */
 150 /*          message related to LWORK is issued by XERBLA. */
 151
 152 /*  INFO    (output) INTEGER */
 153 /*          = 0:  successful exit */
 154 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 155
 156 /*  ===================================================================== */
 157
 158 /*     .. Local Scalars .. */
 159 /*     .. */
 160 /*     .. External Functions .. */
 161 /*     .. */
 162 /*     .. External Subroutines .. */
 163 /*     .. */
 164 /*     .. Intrinsic Functions .. */
 165 /*     .. */
 166 /*     .. Executable Statements .. */
 167
 168 /*     Test the input arguments */
 169
 170     /* Parameter adjustments */
 171     a_dim1 = *lda;
 172     a_offset = 1 + a_dim1;
 173     a -= a_offset;
 174     --tau;
 175     c_dim1 = *ldc;
 176     c_offset = 1 + c_dim1;
 177     c__ -= c_offset;
 178     --work;
 179
 180     /* Function Body */
 181     *info = 0;
 182     applyq = lsame_(vect, "Q");
 183     left = lsame_(side, "L");
 184     notran = lsame_(trans, "N");
 185     lquery = *lwork == -1;
 186
 187 /*     NQ is the order of Q or P and NW is the minimum dimension of WORK */
 188
 189     if (left) {
 190         nq = *m;
 191         nw = *n;
 192     } else {
 193         nq = *n;
 194         nw = *m;
 195     }
 196     if (! applyq && ! lsame_(vect, "P")) {
 197         *info = -1;
 198     } else if (! left && ! lsame_(side, "R")) {
 199         *info = -2;
 200     } else if (! notran && ! lsame_(trans, "T")) {
 201         *info = -3;
 202     } else if (*m < 0) {
 203         *info = -4;
 204     } else if (*n < 0) {
 205         *info = -5;
 206     } else if (*k < 0) {
 207         *info = -6;
 208     } else /* if(complicated condition) */ {
 209 /* Computing MAX */
 210         i__1 = 1, i__2 = min(nq,*k);
 211         if (applyq && *lda < max(1,nq) || ! applyq && *lda < max(i__1,i__2)) {
 212             *info = -8;
 213         } else if (*ldc < max(1,*m)) {
 214             *info = -11;
 215         } else if (*lwork < max(1,nw) && ! lquery) {
 216             *info = -13;
 217         }
 218     }
 219
 220     if (*info == 0) {
 221         if (applyq) {
 222             if (left) {
 223 /* Writing concatenation */
 224                 i__3[0] = 1, a__1[0] = side;
 225                 i__3[1] = 1, a__1[1] = trans;
 226                 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 227                 i__1 = *m - 1;
 228                 i__2 = *m - 1;
 229                 nb = ilaenv_(&c__1, "DORMQR", ch__1, &i__1, n, &i__2, &c_n1);
 230             } else {
 231 /* Writing concatenation */
 232                 i__3[0] = 1, a__1[0] = side;
 233                 i__3[1] = 1, a__1[1] = trans;
 234                 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 235                 i__1 = *n - 1;
 236                 i__2 = *n - 1;
 237                 nb = ilaenv_(&c__1, "DORMQR", ch__1, m, &i__1, &i__2, &c_n1);
 238             }
 239         } else {
 240             if (left) {
 241 /* Writing concatenation */
 242                 i__3[0] = 1, a__1[0] = side;
 243                 i__3[1] = 1, a__1[1] = trans;
 244                 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 245                 i__1 = *m - 1;
 246                 i__2 = *m - 1;
 247                 nb = ilaenv_(&c__1, "DORMLQ", ch__1, &i__1, n, &i__2, &c_n1);
 248             } else {
 249 /* Writing concatenation */
 250                 i__3[0] = 1, a__1[0] = side;
 251                 i__3[1] = 1, a__1[1] = trans;
 252                 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
 253                 i__1 = *n - 1;
 254                 i__2 = *n - 1;
 255                 nb = ilaenv_(&c__1, "DORMLQ", ch__1, m, &i__1, &i__2, &c_n1);
 256             }
 257         }
 258         lwkopt = max(1,nw) * nb;
 259         work[1] = (doublereal) lwkopt;
 260     }
 261
 262     if (*info != 0) {
 263         i__1 = -(*info);
 264         xerbla_("DORMBR", &i__1);
 265         return 0;
 266     } else if (lquery) {
 267         return 0;
 268     }
 269
 270 /*     Quick return if possible */
 271
 272     work[1] = 1.;
 273     if (*m == 0 || *n == 0) {
 274         return 0;
 275     }
 276
 277     if (applyq) {
 278
 279 /*        Apply Q */
 280
 281         if (nq >= *k) {
 282
 283 /*           Q was determined by a call to DGEBRD with nq >= k */
 284
 285             dormqr_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 286                     c_offset], ldc, &work[1], lwork, &iinfo);
 287         } else if (nq > 1) {
 288
 289 /*           Q was determined by a call to DGEBRD with nq < k */
 290
 291             if (left) {
 292                 mi = *m - 1;
 293                 ni = *n;
 294                 i1 = 2;
 295                 i2 = 1;
 296             } else {
 297                 mi = *m;
 298                 ni = *n - 1;
 299                 i1 = 1;
 300                 i2 = 2;
 301             }
 302             i__1 = nq - 1;
 303             dormqr_(side, trans, &mi, &ni, &i__1, &a[a_dim1 + 2], lda, &tau[1]
 304 , &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
 305         }
 306     } else {
 307
 308 /*        Apply P */
 309
 310         if (notran) {
 311             *(unsigned char *)transt = 'T';
 312         } else {
 313             *(unsigned char *)transt = 'N';
 314         }
 315         if (nq > *k) {
 316
 317 /*           P was determined by a call to DGEBRD with nq > k */
 318
 319             dormlq_(side, transt, m, n, k, &a[a_offset], lda, &tau[1], &c__[
 320                     c_offset], ldc, &work[1], lwork, &iinfo);
 321         } else if (nq > 1) {
 322
 323 /*           P was determined by a call to DGEBRD with nq <= k */
 324
 325             if (left) {
 326                 mi = *m - 1;
 327                 ni = *n;
 328                 i1 = 2;
 329                 i2 = 1;
 330             } else {
 331                 mi = *m;
 332                 ni = *n - 1;
 333                 i1 = 1;
 334                 i2 = 2;
 335             }
 336             i__1 = nq - 1;
 337             dormlq_(side, transt, &mi, &ni, &i__1, &a[(a_dim1 << 1) + 1], lda,
 338                      &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &
 339                     iinfo);
 340         }
 341     }
 342     work[1] = (doublereal) lwkopt;
 343     return 0;
 344
 345 /*     End of DORMBR */
 346
 347 } /* dormbr_ */