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

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