00001 /* clarzt.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 complex c_b1 = {0.f,0.f}; 00019 static integer c__1 = 1; 00020 00021 /* Subroutine */ int clarzt_(char *direct, char *storev, integer *n, integer * 00022 k, complex *v, integer *ldv, complex *tau, complex *t, integer *ldt) 00023 { 00024 /* System generated locals */ 00025 integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2; 00026 complex q__1; 00027 00028 /* Local variables */ 00029 integer i__, j, info; 00030 extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex * 00031 , complex *, integer *, complex *, integer *, complex *, complex * 00032 , integer *); 00033 extern logical lsame_(char *, char *); 00034 extern /* Subroutine */ int ctrmv_(char *, char *, char *, integer *, 00035 complex *, integer *, complex *, integer *), clacgv_(integer *, complex *, integer *), xerbla_(char *, 00036 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 /* CLARZT forms the triangular factor T of a complex 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) COMPLEX 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) COMPLEX array, dimension (K) */ 00103 /* TAU(i) must contain the scalar factor of the elementary */ 00104 /* reflector H(i). */ 00105 00106 /* T (output) COMPLEX 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_("CLARZT", &i__1); 00187 return 0; 00188 } 00189 00190 for (i__ = *k; i__ >= 1; --i__) { 00191 i__1 = i__; 00192 if (tau[i__1].r == 0.f && tau[i__1].i == 0.f) { 00193 00194 /* H(i) = I */ 00195 00196 i__1 = *k; 00197 for (j = i__; j <= i__1; ++j) { 00198 i__2 = j + i__ * t_dim1; 00199 t[i__2].r = 0.f, t[i__2].i = 0.f; 00200 /* L10: */ 00201 } 00202 } else { 00203 00204 /* general case */ 00205 00206 if (i__ < *k) { 00207 00208 /* T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)' */ 00209 00210 clacgv_(n, &v[i__ + v_dim1], ldv); 00211 i__1 = *k - i__; 00212 i__2 = i__; 00213 q__1.r = -tau[i__2].r, q__1.i = -tau[i__2].i; 00214 cgemv_("No transpose", &i__1, n, &q__1, &v[i__ + 1 + v_dim1], 00215 ldv, &v[i__ + v_dim1], ldv, &c_b1, &t[i__ + 1 + i__ * 00216 t_dim1], &c__1); 00217 clacgv_(n, &v[i__ + v_dim1], ldv); 00218 00219 /* T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i) */ 00220 00221 i__1 = *k - i__; 00222 ctrmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ + 1 00223 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ * t_dim1] 00224 , &c__1); 00225 } 00226 i__1 = i__ + i__ * t_dim1; 00227 i__2 = i__; 00228 t[i__1].r = tau[i__2].r, t[i__1].i = tau[i__2].i; 00229 } 00230 /* L20: */ 00231 } 00232 return 0; 00233 00234 /* End of CLARZT */ 00235 00236 } /* clarzt_ */