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