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