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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static doublereal c_b8 = 0.;
   7
   8 /* Subroutine */ int dlarft_(char *direct, char *storev, integer *n, integer *
   9         k, doublereal *v, integer *ldv, doublereal *tau, doublereal *t,
  10         integer *ldt)
  11 {
  12     /* System generated locals */
  13     integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3;
  14     doublereal d__1;
  15
  16     /* Local variables */
  17     integer i__, j;
  18     doublereal vii;
  19     extern logical lsame_(char *, char *);
  20     extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
  21             doublereal *, doublereal *, integer *, doublereal *, integer *,
  22             doublereal *, doublereal *, integer *), dtrmv_(char *,
  23             char *, char *, integer *, doublereal *, integer *, doublereal *,
  24             integer *);
  25
  26
  27 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  28 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  29 /*     November 2006 */
  30
  31 /*     .. Scalar Arguments .. */
  32 /*     .. */
  33 /*     .. Array Arguments .. */
  34 /*     .. */
  35
  36 /*  Purpose */
  37 /*  ======= */
  38
  39 /*  DLARFT forms the triangular factor T of a real block reflector H */
  40 /*  of order n, which is defined as a product of k elementary reflectors. */
  41
  42 /*  If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; */
  43
  44 /*  If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. */
  45
  46 /*  If STOREV = 'C', the vector which defines the elementary reflector */
  47 /*  H(i) is stored in the i-th column of the array V, and */
  48
  49 /*     H  =  I - V * T * V' */
  50
  51 /*  If STOREV = 'R', the vector which defines the elementary reflector */
  52 /*  H(i) is stored in the i-th row of the array V, and */
  53
  54 /*     H  =  I - V' * T * V */
  55
  56 /*  Arguments */
  57 /*  ========= */
  58
  59 /*  DIRECT  (input) CHARACTER*1 */
  60 /*          Specifies the order in which the elementary reflectors are */
  61 /*          multiplied to form the block reflector: */
  62 /*          = 'F': H = H(1) H(2) . . . H(k) (Forward) */
  63 /*          = 'B': H = H(k) . . . H(2) H(1) (Backward) */
  64
  65 /*  STOREV  (input) CHARACTER*1 */
  66 /*          Specifies how the vectors which define the elementary */
  67 /*          reflectors are stored (see also Further Details): */
  68 /*          = 'C': columnwise */
  69 /*          = 'R': rowwise */
  70
  71 /*  N       (input) INTEGER */
  72 /*          The order of the block reflector H. N >= 0. */
  73
  74 /*  K       (input) INTEGER */
  75 /*          The order of the triangular factor T (= the number of */
  76 /*          elementary reflectors). K >= 1. */
  77
  78 /*  V       (input/output) DOUBLE PRECISION array, dimension */
  79 /*                               (LDV,K) if STOREV = 'C' */
  80 /*                               (LDV,N) if STOREV = 'R' */
  81 /*          The matrix V. See further details. */
  82
  83 /*  LDV     (input) INTEGER */
  84 /*          The leading dimension of the array V. */
  85 /*          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. */
  86
  87 /*  TAU     (input) DOUBLE PRECISION array, dimension (K) */
  88 /*          TAU(i) must contain the scalar factor of the elementary */
  89 /*          reflector H(i). */
  90
  91 /*  T       (output) DOUBLE PRECISION array, dimension (LDT,K) */
  92 /*          The k by k triangular factor T of the block reflector. */
  93 /*          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is */
  94 /*          lower triangular. The rest of the array is not used. */
  95
  96 /*  LDT     (input) INTEGER */
  97 /*          The leading dimension of the array T. LDT >= K. */
  98
  99 /*  Further Details */
 100 /*  =============== */
 101
 102 /*  The shape of the matrix V and the storage of the vectors which define */
 103 /*  the H(i) is best illustrated by the following example with n = 5 and */
 104 /*  k = 3. The elements equal to 1 are not stored; the corresponding */
 105 /*  array elements are modified but restored on exit. The rest of the */
 106 /*  array is not used. */
 107
 108 /*  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R': */
 109
 110 /*               V = (  1       )                 V = (  1 v1 v1 v1 v1 ) */
 111 /*                   ( v1  1    )                     (     1 v2 v2 v2 ) */
 112 /*                   ( v1 v2  1 )                     (        1 v3 v3 ) */
 113 /*                   ( v1 v2 v3 ) */
 114 /*                   ( v1 v2 v3 ) */
 115
 116 /*  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R': */
 117
 118 /*               V = ( v1 v2 v3 )                 V = ( v1 v1  1       ) */
 119 /*                   ( v1 v2 v3 )                     ( v2 v2 v2  1    ) */
 120 /*                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 ) */
 121 /*                   (     1 v3 ) */
 122 /*                   (        1 ) */
 123
 124 /*  ===================================================================== */
 125
 126 /*     .. Parameters .. */
 127 /*     .. */
 128 /*     .. Local Scalars .. */
 129 /*     .. */
 130 /*     .. External Subroutines .. */
 131 /*     .. */
 132 /*     .. External Functions .. */
 133 /*     .. */
 134 /*     .. Executable Statements .. */
 135
 136 /*     Quick return if possible */
 137
 138     /* Parameter adjustments */
 139     v_dim1 = *ldv;
 140     v_offset = 1 + v_dim1;
 141     v -= v_offset;
 142     --tau;
 143     t_dim1 = *ldt;
 144     t_offset = 1 + t_dim1;
 145     t -= t_offset;
 146
 147     /* Function Body */
 148     if (*n == 0) {
 149         return 0;
 150     }
 151
 152     if (lsame_(direct, "F")) {
 153         i__1 = *k;
 154         for (i__ = 1; i__ <= i__1; ++i__) {
 155             if (tau[i__] == 0.) {
 156
 157 /*              H(i)  =  I */
 158
 159                 i__2 = i__;
 160                 for (j = 1; j <= i__2; ++j) {
 161                     t[j + i__ * t_dim1] = 0.;
 162 /* L10: */
 163                 }
 164             } else {
 165
 166 /*              general case */
 167
 168                 vii = v[i__ + i__ * v_dim1];
 169                 v[i__ + i__ * v_dim1] = 1.;
 170                 if (lsame_(storev, "C")) {
 171
 172 /*                 T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) */
 173
 174                     i__2 = *n - i__ + 1;
 175                     i__3 = i__ - 1;
 176                     d__1 = -tau[i__];
 177                     dgemv_("Transpose", &i__2, &i__3, &d__1, &v[i__ + v_dim1],
 178                              ldv, &v[i__ + i__ * v_dim1], &c__1, &c_b8, &t[
 179                             i__ * t_dim1 + 1], &c__1);
 180                 } else {
 181
 182 /*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' */
 183
 184                     i__2 = i__ - 1;
 185                     i__3 = *n - i__ + 1;
 186                     d__1 = -tau[i__];
 187                     dgemv_("No transpose", &i__2, &i__3, &d__1, &v[i__ *
 188                             v_dim1 + 1], ldv, &v[i__ + i__ * v_dim1], ldv, &
 189                             c_b8, &t[i__ * t_dim1 + 1], &c__1);
 190                 }
 191                 v[i__ + i__ * v_dim1] = vii;
 192
 193 /*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) */
 194
 195                 i__2 = i__ - 1;
 196                 dtrmv_("Upper", "No transpose", "Non-unit", &i__2, &t[
 197                         t_offset], ldt, &t[i__ * t_dim1 + 1], &c__1);
 198                 t[i__ + i__ * t_dim1] = tau[i__];
 199             }
 200 /* L20: */
 201         }
 202     } else {
 203         for (i__ = *k; i__ >= 1; --i__) {
 204             if (tau[i__] == 0.) {
 205
 206 /*              H(i)  =  I */
 207
 208                 i__1 = *k;
 209                 for (j = i__; j <= i__1; ++j) {
 210                     t[j + i__ * t_dim1] = 0.;
 211 /* L30: */
 212                 }
 213             } else {
 214
 215 /*              general case */
 216
 217                 if (i__ < *k) {
 218                     if (lsame_(storev, "C")) {
 219                         vii = v[*n - *k + i__ + i__ * v_dim1];
 220                         v[*n - *k + i__ + i__ * v_dim1] = 1.;
 221
 222 /*                    T(i+1:k,i) := */
 223 /*                            - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) */
 224
 225                         i__1 = *n - *k + i__;
 226                         i__2 = *k - i__;
 227                         d__1 = -tau[i__];
 228                         dgemv_("Transpose", &i__1, &i__2, &d__1, &v[(i__ + 1)
 229                                 * v_dim1 + 1], ldv, &v[i__ * v_dim1 + 1], &
 230                                 c__1, &c_b8, &t[i__ + 1 + i__ * t_dim1], &
 231                                 c__1);
 232                         v[*n - *k + i__ + i__ * v_dim1] = vii;
 233                     } else {
 234                         vii = v[i__ + (*n - *k + i__) * v_dim1];
 235                         v[i__ + (*n - *k + i__) * v_dim1] = 1.;
 236
 237 /*                    T(i+1:k,i) := */
 238 /*                            - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' */
 239
 240                         i__1 = *k - i__;
 241                         i__2 = *n - *k + i__;
 242                         d__1 = -tau[i__];
 243                         dgemv_("No transpose", &i__1, &i__2, &d__1, &v[i__ +
 244                                 1 + v_dim1], ldv, &v[i__ + v_dim1], ldv, &
 245                                 c_b8, &t[i__ + 1 + i__ * t_dim1], &c__1);
 246                         v[i__ + (*n - *k + i__) * v_dim1] = vii;
 247                     }
 248
 249 /*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) */
 250
 251                     i__1 = *k - i__;
 252                     dtrmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__
 253                             + 1 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ *
 254                              t_dim1], &c__1)
 255                             ;
 256                 }
 257                 t[i__ + i__ * t_dim1] = tau[i__];
 258             }
 259 /* L40: */
 260         }
 261     }
 262     return 0;
 263
 264 /*     End of DLARFT */
 265
 266 } /* dlarft_ */