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

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