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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dorgl2_(integer *m, integer *n, integer *k, doublereal *
   4         a, integer *lda, doublereal *tau, doublereal *work, integer *info)
   5 {
   6     /* System generated locals */
   7     integer a_dim1, a_offset, i__1, i__2;
   8     doublereal d__1;
   9
  10     /* Local variables */
  11     integer i__, j, l;
  12     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  13             integer *), dlarf_(char *, integer *, integer *, doublereal *,
  14             integer *, doublereal *, doublereal *, integer *, doublereal *), xerbla_(char *, integer *);
  15
  16
  17 /*  -- LAPACK routine (version 3.1) -- */
  18 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  19 /*     November 2006 */
  20
  21 /*     .. Scalar Arguments .. */
  22 /*     .. */
  23 /*     .. Array Arguments .. */
  24 /*     .. */
  25
  26 /*  Purpose */
  27 /*  ======= */
  28
  29 /*  DORGL2 generates an m by n real matrix Q with orthonormal rows, */
  30 /*  which is defined as the first m rows of a product of k elementary */
  31 /*  reflectors of order n */
  32
  33 /*        Q  =  H(k) . . . H(2) H(1) */
  34
  35 /*  as returned by DGELQF. */
  36
  37 /*  Arguments */
  38 /*  ========= */
  39
  40 /*  M       (input) INTEGER */
  41 /*          The number of rows of the matrix Q. M >= 0. */
  42
  43 /*  N       (input) INTEGER */
  44 /*          The number of columns of the matrix Q. N >= M. */
  45
  46 /*  K       (input) INTEGER */
  47 /*          The number of elementary reflectors whose product defines the */
  48 /*          matrix Q. M >= K >= 0. */
  49
  50 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
  51 /*          On entry, the i-th row must contain the vector which defines */
  52 /*          the elementary reflector H(i), for i = 1,2,...,k, as returned */
  53 /*          by DGELQF in the first k rows of its array argument A. */
  54 /*          On exit, the m-by-n matrix Q. */
  55
  56 /*  LDA     (input) INTEGER */
  57 /*          The first dimension of the array A. LDA >= max(1,M). */
  58
  59 /*  TAU     (input) DOUBLE PRECISION array, dimension (K) */
  60 /*          TAU(i) must contain the scalar factor of the elementary */
  61 /*          reflector H(i), as returned by DGELQF. */
  62
  63 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (M) */
  64
  65 /*  INFO    (output) INTEGER */
  66 /*          = 0: successful exit */
  67 /*          < 0: if INFO = -i, the i-th argument has an illegal value */
  68
  69 /*  ===================================================================== */
  70
  71 /*     .. Parameters .. */
  72 /*     .. */
  73 /*     .. Local Scalars .. */
  74 /*     .. */
  75 /*     .. External Subroutines .. */
  76 /*     .. */
  77 /*     .. Intrinsic Functions .. */
  78 /*     .. */
  79 /*     .. Executable Statements .. */
  80
  81 /*     Test the input arguments */
  82
  83     /* Parameter adjustments */
  84     a_dim1 = *lda;
  85     a_offset = 1 + a_dim1;
  86     a -= a_offset;
  87     --tau;
  88     --work;
  89
  90     /* Function Body */
  91     *info = 0;
  92     if (*m < 0) {
  93         *info = -1;
  94     } else if (*n < *m) {
  95         *info = -2;
  96     } else if (*k < 0 || *k > *m) {
  97         *info = -3;
  98     } else if (*lda < max(1,*m)) {
  99         *info = -5;
 100     }
 101     if (*info != 0) {
 102         i__1 = -(*info);
 103         xerbla_("DORGL2", &i__1);
 104         return 0;
 105     }
 106
 107 /*     Quick return if possible */
 108
 109     if (*m <= 0) {
 110         return 0;
 111     }
 112
 113     if (*k < *m) {
 114
 115 /*        Initialise rows k+1:m to rows of the unit matrix */
 116
 117         i__1 = *n;
 118         for (j = 1; j <= i__1; ++j) {
 119             i__2 = *m;
 120             for (l = *k + 1; l <= i__2; ++l) {
 121                 a[l + j * a_dim1] = 0.;
 122 /* L10: */
 123             }
 124             if (j > *k && j <= *m) {
 125                 a[j + j * a_dim1] = 1.;
 126             }
 127 /* L20: */
 128         }
 129     }
 130
 131     for (i__ = *k; i__ >= 1; --i__) {
 132
 133 /*        Apply H(i) to A(i:m,i:n) from the right */
 134
 135         if (i__ < *n) {
 136             if (i__ < *m) {
 137                 a[i__ + i__ * a_dim1] = 1.;
 138                 i__1 = *m - i__;
 139                 i__2 = *n - i__ + 1;
 140                 dlarf_("Right", &i__1, &i__2, &a[i__ + i__ * a_dim1], lda, &
 141                         tau[i__], &a[i__ + 1 + i__ * a_dim1], lda, &work[1]);
 142             }
 143             i__1 = *n - i__;
 144             d__1 = -tau[i__];
 145             dscal_(&i__1, &d__1, &a[i__ + (i__ + 1) * a_dim1], lda);
 146         }
 147         a[i__ + i__ * a_dim1] = 1. - tau[i__];
 148
 149 /*        Set A(i,1:i-1) to zero */
 150
 151         i__1 = i__ - 1;
 152         for (l = 1; l <= i__1; ++l) {
 153             a[i__ + l * a_dim1] = 0.;
 154 /* L30: */
 155         }
 156 /* L40: */
 157     }
 158     return 0;
 159
 160 /*     End of DORGL2 */
 161
 162 } /* dorgl2_ */