1 /* ./src_f77/slarzt.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static real c_b8 = 0.f;
11 static integer c__1 = 1;
12
slarzt_(char * direct,char * storev,integer * n,integer * k,real * v,integer * ldv,real * tau,real * t,integer * ldt,ftnlen direct_len,ftnlen storev_len)13 /* Subroutine */ int slarzt_(char *direct, char *storev, integer *n, integer *
14 k, real *v, integer *ldv, real *tau, real *t, integer *ldt, ftnlen
15 direct_len, ftnlen storev_len)
16 {
17 /* System generated locals */
18 integer t_dim1, t_offset, v_dim1, v_offset, i__1;
19 real r__1;
20
21 /* Local variables */
22 static integer i__, j, info;
23 extern logical lsame_(char *, char *, ftnlen, ftnlen);
24 extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *,
25 real *, integer *, real *, integer *, real *, real *, integer *,
26 ftnlen), strmv_(char *, char *, char *, integer *, real *,
27 integer *, real *, integer *, ftnlen, ftnlen, ftnlen), xerbla_(
28 char *, integer *, ftnlen);
29
30
31 /* -- LAPACK routine (version 3.0) -- */
32 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
33 /* Courant Institute, Argonne National Lab, and Rice University */
34 /* June 30, 1999 */
35
36 /* .. Scalar Arguments .. */
37 /* .. */
38 /* .. Array Arguments .. */
39 /* .. */
40
41 /* Purpose */
42 /* ======= */
43
44 /* SLARZT forms the triangular factor T of a real block reflector */
45 /* H of order > n, which is defined as a product of k elementary */
46 /* reflectors. */
47
48 /* If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; */
49
50 /* If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. */
51
52 /* If STOREV = 'C', the vector which defines the elementary reflector */
53 /* H(i) is stored in the i-th column of the array V, and */
54
55 /* H = I - V * T * V' */
56
57 /* If STOREV = 'R', the vector which defines the elementary reflector */
58 /* H(i) is stored in the i-th row of the array V, and */
59
60 /* H = I - V' * T * V */
61
62 /* Currently, only STOREV = 'R' and DIRECT = 'B' are supported. */
63
64 /* Arguments */
65 /* ========= */
66
67 /* DIRECT (input) CHARACTER*1 */
68 /* Specifies the order in which the elementary reflectors are */
69 /* multiplied to form the block reflector: */
70 /* = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet) */
71 /* = 'B': H = H(k) . . . H(2) H(1) (Backward) */
72
73 /* STOREV (input) CHARACTER*1 */
74 /* Specifies how the vectors which define the elementary */
75 /* reflectors are stored (see also Further Details): */
76 /* = 'C': columnwise (not supported yet) */
77 /* = 'R': rowwise */
78
79 /* N (input) INTEGER */
80 /* The order of the block reflector H. N >= 0. */
81
82 /* K (input) INTEGER */
83 /* The order of the triangular factor T (= the number of */
84 /* elementary reflectors). K >= 1. */
85
86 /* V (input/output) REAL array, dimension */
87 /* (LDV,K) if STOREV = 'C' */
88 /* (LDV,N) if STOREV = 'R' */
89 /* The matrix V. See further details. */
90
91 /* LDV (input) INTEGER */
92 /* The leading dimension of the array V. */
93 /* If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. */
94
95 /* TAU (input) REAL array, dimension (K) */
96 /* TAU(i) must contain the scalar factor of the elementary */
97 /* reflector H(i). */
98
99 /* T (output) REAL array, dimension (LDT,K) */
100 /* The k by k triangular factor T of the block reflector. */
101 /* If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is */
102 /* lower triangular. The rest of the array is not used. */
103
104 /* LDT (input) INTEGER */
105 /* The leading dimension of the array T. LDT >= K. */
106
107 /* Further Details */
108 /* =============== */
109
110 /* Based on contributions by */
111 /* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */
112
113 /* The shape of the matrix V and the storage of the vectors which define */
114 /* the H(i) is best illustrated by the following example with n = 5 and */
115 /* k = 3. The elements equal to 1 are not stored; the corresponding */
116 /* array elements are modified but restored on exit. The rest of the */
117 /* array is not used. */
118
119 /* DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': */
120
121 /* ______V_____ */
122 /* ( v1 v2 v3 ) / \ */
123 /* ( v1 v2 v3 ) ( v1 v1 v1 v1 v1 . . . . 1 ) */
124 /* V = ( v1 v2 v3 ) ( v2 v2 v2 v2 v2 . . . 1 ) */
125 /* ( v1 v2 v3 ) ( v3 v3 v3 v3 v3 . . 1 ) */
126 /* ( v1 v2 v3 ) */
127 /* . . . */
128 /* . . . */
129 /* 1 . . */
130 /* 1 . */
131 /* 1 */
132
133 /* DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': */
134
135 /* ______V_____ */
136 /* 1 / \ */
137 /* . 1 ( 1 . . . . v1 v1 v1 v1 v1 ) */
138 /* . . 1 ( . 1 . . . v2 v2 v2 v2 v2 ) */
139 /* . . . ( . . 1 . . v3 v3 v3 v3 v3 ) */
140 /* . . . */
141 /* ( v1 v2 v3 ) */
142 /* ( v1 v2 v3 ) */
143 /* V = ( v1 v2 v3 ) */
144 /* ( v1 v2 v3 ) */
145 /* ( v1 v2 v3 ) */
146
147 /* ===================================================================== */
148
149 /* .. Parameters .. */
150 /* .. */
151 /* .. Local Scalars .. */
152 /* .. */
153 /* .. External Subroutines .. */
154 /* .. */
155 /* .. External Functions .. */
156 /* .. */
157 /* .. Executable Statements .. */
158
159 /* Check for currently supported options */
160
161 /* Parameter adjustments */
162 v_dim1 = *ldv;
163 v_offset = 1 + v_dim1;
164 v -= v_offset;
165 --tau;
166 t_dim1 = *ldt;
167 t_offset = 1 + t_dim1;
168 t -= t_offset;
169
170 /* Function Body */
171 info = 0;
172 if (! lsame_(direct, "B", (ftnlen)1, (ftnlen)1)) {
173 info = -1;
174 } else if (! lsame_(storev, "R", (ftnlen)1, (ftnlen)1)) {
175 info = -2;
176 }
177 if (info != 0) {
178 i__1 = -info;
179 xerbla_("SLARZT", &i__1, (ftnlen)6);
180 return 0;
181 }
182
183 for (i__ = *k; i__ >= 1; --i__) {
184 if (tau[i__] == 0.f) {
185
186 /* H(i) = I */
187
188 i__1 = *k;
189 for (j = i__; j <= i__1; ++j) {
190 t[j + i__ * t_dim1] = 0.f;
191 /* L10: */
192 }
193 } else {
194
195 /* general case */
196
197 if (i__ < *k) {
198
199 /* T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)' */
200
201 i__1 = *k - i__;
202 r__1 = -tau[i__];
203 sgemv_("No transpose", &i__1, n, &r__1, &v[i__ + 1 + v_dim1],
204 ldv, &v[i__ + v_dim1], ldv, &c_b8, &t[i__ + 1 + i__ *
205 t_dim1], &c__1, (ftnlen)12);
206
207 /* T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i) */
208
209 i__1 = *k - i__;
210 strmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ + 1
211 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ * t_dim1]
212 , &c__1, (ftnlen)5, (ftnlen)12, (ftnlen)8);
213 }
214 t[i__ + i__ * t_dim1] = tau[i__];
215 }
216 /* L20: */
217 }
218 return 0;
219
220 /* End of SLARZT */
221
222 } /* slarzt_ */
223
224