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

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