git.maemo.org Git - opencv/blob - 3rdparty/lapack/dorm2r.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6
   7 /* Subroutine */ int dorm2r_(char *side, char *trans, integer *m, integer *n,
   8         integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *
   9         c__, integer *ldc, doublereal *work, integer *info)
  10 {
  11     /* System generated locals */
  12     integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2;
  13
  14     /* Local variables */
  15     integer i__, i1, i2, i3, ic, jc, mi, ni, nq;
  16     doublereal aii;
  17     logical left;
  18     extern /* Subroutine */ int dlarf_(char *, integer *, integer *,
  19             doublereal *, integer *, doublereal *, doublereal *, integer *,
  20             doublereal *);
  21     extern logical lsame_(char *, char *);
  22     extern /* Subroutine */ int xerbla_(char *, integer *);
  23     logical notran;
  24
  25
  26 /*  -- LAPACK routine (version 3.1) -- */
  27 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  28 /*     November 2006 */
  29
  30 /*     .. Scalar Arguments .. */
  31 /*     .. */
  32 /*     .. Array Arguments .. */
  33 /*     .. */
  34
  35 /*  Purpose */
  36 /*  ======= */
  37
  38 /*  DORM2R overwrites the general real m by n matrix C with */
  39
  40 /*        Q * C  if SIDE = 'L' and TRANS = 'N', or */
  41
  42 /*        Q'* C  if SIDE = 'L' and TRANS = 'T', or */
  43
  44 /*        C * Q  if SIDE = 'R' and TRANS = 'N', or */
  45
  46 /*        C * Q' if SIDE = 'R' and TRANS = 'T', */
  47
  48 /*  where Q is a real orthogonal matrix defined as the product of k */
  49 /*  elementary reflectors */
  50
  51 /*        Q = H(1) H(2) . . . H(k) */
  52
  53 /*  as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n */
  54 /*  if SIDE = 'R'. */
  55
  56 /*  Arguments */
  57 /*  ========= */
  58
  59 /*  SIDE    (input) CHARACTER*1 */
  60 /*          = 'L': apply Q or Q' from the Left */
  61 /*          = 'R': apply Q or Q' from the Right */
  62
  63 /*  TRANS   (input) CHARACTER*1 */
  64 /*          = 'N': apply Q  (No transpose) */
  65 /*          = 'T': apply Q' (Transpose) */
  66
  67 /*  M       (input) INTEGER */
  68 /*          The number of rows of the matrix C. M >= 0. */
  69
  70 /*  N       (input) INTEGER */
  71 /*          The number of columns of the matrix C. N >= 0. */
  72
  73 /*  K       (input) INTEGER */
  74 /*          The number of elementary reflectors whose product defines */
  75 /*          the matrix Q. */
  76 /*          If SIDE = 'L', M >= K >= 0; */
  77 /*          if SIDE = 'R', N >= K >= 0. */
  78
  79 /*  A       (input) DOUBLE PRECISION array, dimension (LDA,K) */
  80 /*          The i-th column must contain the vector which defines the */
  81 /*          elementary reflector H(i), for i = 1,2,...,k, as returned by */
  82 /*          DGEQRF in the first k columns of its array argument A. */
  83 /*          A is modified by the routine but restored on exit. */
  84
  85 /*  LDA     (input) INTEGER */
  86 /*          The leading dimension of the array A. */
  87 /*          If SIDE = 'L', LDA >= max(1,M); */
  88 /*          if SIDE = 'R', LDA >= max(1,N). */
  89
  90 /*  TAU     (input) DOUBLE PRECISION array, dimension (K) */
  91 /*          TAU(i) must contain the scalar factor of the elementary */
  92 /*          reflector H(i), as returned by DGEQRF. */
  93
  94 /*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
  95 /*          On entry, the m by n matrix C. */
  96 /*          On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. */
  97
  98 /*  LDC     (input) INTEGER */
  99 /*          The leading dimension of the array C. LDC >= max(1,M). */
 100
 101 /*  WORK    (workspace) DOUBLE PRECISION array, dimension */
 102 /*                                   (N) if SIDE = 'L', */
 103 /*                                   (M) if SIDE = 'R' */
 104
 105 /*  INFO    (output) INTEGER */
 106 /*          = 0: successful exit */
 107 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
 108
 109 /*  ===================================================================== */
 110
 111 /*     .. Parameters .. */
 112 /*     .. */
 113 /*     .. Local Scalars .. */
 114 /*     .. */
 115 /*     .. External Functions .. */
 116 /*     .. */
 117 /*     .. External Subroutines .. */
 118 /*     .. */
 119 /*     .. Intrinsic Functions .. */
 120 /*     .. */
 121 /*     .. Executable Statements .. */
 122
 123 /*     Test the input arguments */
 124
 125     /* Parameter adjustments */
 126     a_dim1 = *lda;
 127     a_offset = 1 + a_dim1;
 128     a -= a_offset;
 129     --tau;
 130     c_dim1 = *ldc;
 131     c_offset = 1 + c_dim1;
 132     c__ -= c_offset;
 133     --work;
 134
 135     /* Function Body */
 136     *info = 0;
 137     left = lsame_(side, "L");
 138     notran = lsame_(trans, "N");
 139
 140 /*     NQ is the order of Q */
 141
 142     if (left) {
 143         nq = *m;
 144     } else {
 145         nq = *n;
 146     }
 147     if (! left && ! lsame_(side, "R")) {
 148         *info = -1;
 149     } else if (! notran && ! lsame_(trans, "T")) {
 150         *info = -2;
 151     } else if (*m < 0) {
 152         *info = -3;
 153     } else if (*n < 0) {
 154         *info = -4;
 155     } else if (*k < 0 || *k > nq) {
 156         *info = -5;
 157     } else if (*lda < max(1,nq)) {
 158         *info = -7;
 159     } else if (*ldc < max(1,*m)) {
 160         *info = -10;
 161     }
 162     if (*info != 0) {
 163         i__1 = -(*info);
 164         xerbla_("DORM2R", &i__1);
 165         return 0;
 166     }
 167
 168 /*     Quick return if possible */
 169
 170     if (*m == 0 || *n == 0 || *k == 0) {
 171         return 0;
 172     }
 173
 174     if (left && ! notran || ! left && notran) {
 175         i1 = 1;
 176         i2 = *k;
 177         i3 = 1;
 178     } else {
 179         i1 = *k;
 180         i2 = 1;
 181         i3 = -1;
 182     }
 183
 184     if (left) {
 185         ni = *n;
 186         jc = 1;
 187     } else {
 188         mi = *m;
 189         ic = 1;
 190     }
 191
 192     i__1 = i2;
 193     i__2 = i3;
 194     for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 195         if (left) {
 196
 197 /*           H(i) is applied to C(i:m,1:n) */
 198
 199             mi = *m - i__ + 1;
 200             ic = i__;
 201         } else {
 202
 203 /*           H(i) is applied to C(1:m,i:n) */
 204
 205             ni = *n - i__ + 1;
 206             jc = i__;
 207         }
 208
 209 /*        Apply H(i) */
 210
 211         aii = a[i__ + i__ * a_dim1];
 212         a[i__ + i__ * a_dim1] = 1.;
 213         dlarf_(side, &mi, &ni, &a[i__ + i__ * a_dim1], &c__1, &tau[i__], &c__[
 214                 ic + jc * c_dim1], ldc, &work[1]);
 215         a[i__ + i__ * a_dim1] = aii;
 216 /* L10: */
 217     }
 218     return 0;
 219
 220 /*     End of DORM2R */
 221
 222 } /* dorm2r_ */