00001 /* slarzt.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static real c_b8 = 0.f; 00019 static integer c__1 = 1; 00020 00021 /* Subroutine */ int slarzt_(char *direct, char *storev, integer *n, integer * 00022 k, real *v, integer *ldv, real *tau, real *t, integer *ldt) 00023 { 00024 /* System generated locals */ 00025 integer t_dim1, t_offset, v_dim1, v_offset, i__1; 00026 real r__1; 00027 00028 /* Local variables */ 00029 integer i__, j, info; 00030 extern logical lsame_(char *, char *); 00031 extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 00032 real *, integer *, real *, integer *, real *, real *, integer *), strmv_(char *, char *, char *, integer *, real *, 00033 integer *, real *, integer *), xerbla_( 00034 char *, integer *); 00035 00036 00037 /* -- LAPACK routine (version 3.2) -- */ 00038 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00039 /* November 2006 */ 00040 00041 /* .. Scalar Arguments .. */ 00042 /* .. */ 00043 /* .. Array Arguments .. */ 00044 /* .. */ 00045 00046 /* Purpose */ 00047 /* ======= */ 00048 00049 /* SLARZT forms the triangular factor T of a real block reflector */ 00050 /* H of order > n, which is defined as a product of k elementary */ 00051 /* reflectors. */ 00052 00053 /* If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; */ 00054 00055 /* If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. */ 00056 00057 /* If STOREV = 'C', the vector which defines the elementary reflector */ 00058 /* H(i) is stored in the i-th column of the array V, and */ 00059 00060 /* H = I - V * T * V' */ 00061 00062 /* If STOREV = 'R', the vector which defines the elementary reflector */ 00063 /* H(i) is stored in the i-th row of the array V, and */ 00064 00065 /* H = I - V' * T * V */ 00066 00067 /* Currently, only STOREV = 'R' and DIRECT = 'B' are supported. */ 00068 00069 /* Arguments */ 00070 /* ========= */ 00071 00072 /* DIRECT (input) CHARACTER*1 */ 00073 /* Specifies the order in which the elementary reflectors are */ 00074 /* multiplied to form the block reflector: */ 00075 /* = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet) */ 00076 /* = 'B': H = H(k) . . . H(2) H(1) (Backward) */ 00077 00078 /* STOREV (input) CHARACTER*1 */ 00079 /* Specifies how the vectors which define the elementary */ 00080 /* reflectors are stored (see also Further Details): */ 00081 /* = 'C': columnwise (not supported yet) */ 00082 /* = 'R': rowwise */ 00083 00084 /* N (input) INTEGER */ 00085 /* The order of the block reflector H. N >= 0. */ 00086 00087 /* K (input) INTEGER */ 00088 /* The order of the triangular factor T (= the number of */ 00089 /* elementary reflectors). K >= 1. */ 00090 00091 /* V (input/output) REAL array, dimension */ 00092 /* (LDV,K) if STOREV = 'C' */ 00093 /* (LDV,N) if STOREV = 'R' */ 00094 /* The matrix V. See further details. */ 00095 00096 /* LDV (input) INTEGER */ 00097 /* The leading dimension of the array V. */ 00098 /* If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. */ 00099 00100 /* TAU (input) REAL array, dimension (K) */ 00101 /* TAU(i) must contain the scalar factor of the elementary */ 00102 /* reflector H(i). */ 00103 00104 /* T (output) REAL array, dimension (LDT,K) */ 00105 /* The k by k triangular factor T of the block reflector. */ 00106 /* If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is */ 00107 /* lower triangular. The rest of the array is not used. */ 00108 00109 /* LDT (input) INTEGER */ 00110 /* The leading dimension of the array T. LDT >= K. */ 00111 00112 /* Further Details */ 00113 /* =============== */ 00114 00115 /* Based on contributions by */ 00116 /* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */ 00117 00118 /* The shape of the matrix V and the storage of the vectors which define */ 00119 /* the H(i) is best illustrated by the following example with n = 5 and */ 00120 /* k = 3. The elements equal to 1 are not stored; the corresponding */ 00121 /* array elements are modified but restored on exit. The rest of the */ 00122 /* array is not used. */ 00123 00124 /* DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': */ 00125 00126 /* ______V_____ */ 00127 /* ( v1 v2 v3 ) / \ */ 00128 /* ( v1 v2 v3 ) ( v1 v1 v1 v1 v1 . . . . 1 ) */ 00129 /* V = ( v1 v2 v3 ) ( v2 v2 v2 v2 v2 . . . 1 ) */ 00130 /* ( v1 v2 v3 ) ( v3 v3 v3 v3 v3 . . 1 ) */ 00131 /* ( v1 v2 v3 ) */ 00132 /* . . . */ 00133 /* . . . */ 00134 /* 1 . . */ 00135 /* 1 . */ 00136 /* 1 */ 00137 00138 /* DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': */ 00139 00140 /* ______V_____ */ 00141 /* 1 / \ */ 00142 /* . 1 ( 1 . . . . v1 v1 v1 v1 v1 ) */ 00143 /* . . 1 ( . 1 . . . v2 v2 v2 v2 v2 ) */ 00144 /* . . . ( . . 1 . . v3 v3 v3 v3 v3 ) */ 00145 /* . . . */ 00146 /* ( v1 v2 v3 ) */ 00147 /* ( v1 v2 v3 ) */ 00148 /* V = ( v1 v2 v3 ) */ 00149 /* ( v1 v2 v3 ) */ 00150 /* ( v1 v2 v3 ) */ 00151 00152 /* ===================================================================== */ 00153 00154 /* .. Parameters .. */ 00155 /* .. */ 00156 /* .. Local Scalars .. */ 00157 /* .. */ 00158 /* .. External Subroutines .. */ 00159 /* .. */ 00160 /* .. External Functions .. */ 00161 /* .. */ 00162 /* .. Executable Statements .. */ 00163 00164 /* Check for currently supported options */ 00165 00166 /* Parameter adjustments */ 00167 v_dim1 = *ldv; 00168 v_offset = 1 + v_dim1; 00169 v -= v_offset; 00170 --tau; 00171 t_dim1 = *ldt; 00172 t_offset = 1 + t_dim1; 00173 t -= t_offset; 00174 00175 /* Function Body */ 00176 info = 0; 00177 if (! lsame_(direct, "B")) { 00178 info = -1; 00179 } else if (! lsame_(storev, "R")) { 00180 info = -2; 00181 } 00182 if (info != 0) { 00183 i__1 = -info; 00184 xerbla_("SLARZT", &i__1); 00185 return 0; 00186 } 00187 00188 for (i__ = *k; i__ >= 1; --i__) { 00189 if (tau[i__] == 0.f) { 00190 00191 /* H(i) = I */ 00192 00193 i__1 = *k; 00194 for (j = i__; j <= i__1; ++j) { 00195 t[j + i__ * t_dim1] = 0.f; 00196 /* L10: */ 00197 } 00198 } else { 00199 00200 /* general case */ 00201 00202 if (i__ < *k) { 00203 00204 /* T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)' */ 00205 00206 i__1 = *k - i__; 00207 r__1 = -tau[i__]; 00208 sgemv_("No transpose", &i__1, n, &r__1, &v[i__ + 1 + v_dim1], 00209 ldv, &v[i__ + v_dim1], ldv, &c_b8, &t[i__ + 1 + i__ * 00210 t_dim1], &c__1); 00211 00212 /* T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i) */ 00213 00214 i__1 = *k - i__; 00215 strmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ + 1 00216 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ * t_dim1] 00217 , &c__1); 00218 } 00219 t[i__ + i__ * t_dim1] = tau[i__]; 00220 } 00221 /* L20: */ 00222 } 00223 return 0; 00224 00225 /* End of SLARZT */ 00226 00227 } /* slarzt_ */