1 | *> \brief \b STZRZF |
---|
2 | * |
---|
3 | * =========== DOCUMENTATION =========== |
---|
4 | * |
---|
5 | * Online html documentation available at |
---|
6 | * http://www.netlib.org/lapack/explore-html/ |
---|
7 | * |
---|
8 | *> \htmlonly |
---|
9 | *> Download STZRZF + dependencies |
---|
10 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stzrzf.f"> |
---|
11 | *> [TGZ]</a> |
---|
12 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stzrzf.f"> |
---|
13 | *> [ZIP]</a> |
---|
14 | *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stzrzf.f"> |
---|
15 | *> [TXT]</a> |
---|
16 | *> \endhtmlonly |
---|
17 | * |
---|
18 | * Definition: |
---|
19 | * =========== |
---|
20 | * |
---|
21 | * SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) |
---|
22 | * |
---|
23 | * .. Scalar Arguments .. |
---|
24 | * INTEGER INFO, LDA, LWORK, M, N |
---|
25 | * .. |
---|
26 | * .. Array Arguments .. |
---|
27 | * REAL A( LDA, * ), TAU( * ), WORK( * ) |
---|
28 | * .. |
---|
29 | * |
---|
30 | * |
---|
31 | *> \par Purpose: |
---|
32 | * ============= |
---|
33 | *> |
---|
34 | *> \verbatim |
---|
35 | *> |
---|
36 | *> STZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A |
---|
37 | *> to upper triangular form by means of orthogonal transformations. |
---|
38 | *> |
---|
39 | *> The upper trapezoidal matrix A is factored as |
---|
40 | *> |
---|
41 | *> A = ( R 0 ) * Z, |
---|
42 | *> |
---|
43 | *> where Z is an N-by-N orthogonal matrix and R is an M-by-M upper |
---|
44 | *> triangular matrix. |
---|
45 | *> \endverbatim |
---|
46 | * |
---|
47 | * Arguments: |
---|
48 | * ========== |
---|
49 | * |
---|
50 | *> \param[in] M |
---|
51 | *> \verbatim |
---|
52 | *> M is INTEGER |
---|
53 | *> The number of rows of the matrix A. M >= 0. |
---|
54 | *> \endverbatim |
---|
55 | *> |
---|
56 | *> \param[in] N |
---|
57 | *> \verbatim |
---|
58 | *> N is INTEGER |
---|
59 | *> The number of columns of the matrix A. N >= M. |
---|
60 | *> \endverbatim |
---|
61 | *> |
---|
62 | *> \param[in,out] A |
---|
63 | *> \verbatim |
---|
64 | *> A is REAL array, dimension (LDA,N) |
---|
65 | *> On entry, the leading M-by-N upper trapezoidal part of the |
---|
66 | *> array A must contain the matrix to be factorized. |
---|
67 | *> On exit, the leading M-by-M upper triangular part of A |
---|
68 | *> contains the upper triangular matrix R, and elements M+1 to |
---|
69 | *> N of the first M rows of A, with the array TAU, represent the |
---|
70 | *> orthogonal matrix Z as a product of M elementary reflectors. |
---|
71 | *> \endverbatim |
---|
72 | *> |
---|
73 | *> \param[in] LDA |
---|
74 | *> \verbatim |
---|
75 | *> LDA is INTEGER |
---|
76 | *> The leading dimension of the array A. LDA >= max(1,M). |
---|
77 | *> \endverbatim |
---|
78 | *> |
---|
79 | *> \param[out] TAU |
---|
80 | *> \verbatim |
---|
81 | *> TAU is REAL array, dimension (M) |
---|
82 | *> The scalar factors of the elementary reflectors. |
---|
83 | *> \endverbatim |
---|
84 | *> |
---|
85 | *> \param[out] WORK |
---|
86 | *> \verbatim |
---|
87 | *> WORK is REAL array, dimension (MAX(1,LWORK)) |
---|
88 | *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
---|
89 | *> \endverbatim |
---|
90 | *> |
---|
91 | *> \param[in] LWORK |
---|
92 | *> \verbatim |
---|
93 | *> LWORK is INTEGER |
---|
94 | *> The dimension of the array WORK. LWORK >= max(1,M). |
---|
95 | *> For optimum performance LWORK >= M*NB, where NB is |
---|
96 | *> the optimal blocksize. |
---|
97 | *> |
---|
98 | *> If LWORK = -1, then a workspace query is assumed; the routine |
---|
99 | *> only calculates the optimal size of the WORK array, returns |
---|
100 | *> this value as the first entry of the WORK array, and no error |
---|
101 | *> message related to LWORK is issued by XERBLA. |
---|
102 | *> \endverbatim |
---|
103 | *> |
---|
104 | *> \param[out] INFO |
---|
105 | *> \verbatim |
---|
106 | *> INFO is INTEGER |
---|
107 | *> = 0: successful exit |
---|
108 | *> < 0: if INFO = -i, the i-th argument had an illegal value |
---|
109 | *> \endverbatim |
---|
110 | * |
---|
111 | * Authors: |
---|
112 | * ======== |
---|
113 | * |
---|
114 | *> \author Univ. of Tennessee |
---|
115 | *> \author Univ. of California Berkeley |
---|
116 | *> \author Univ. of Colorado Denver |
---|
117 | *> \author NAG Ltd. |
---|
118 | * |
---|
119 | *> \date November 2011 |
---|
120 | * |
---|
121 | *> \ingroup realOTHERcomputational |
---|
122 | * |
---|
123 | *> \par Contributors: |
---|
124 | * ================== |
---|
125 | *> |
---|
126 | *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA |
---|
127 | * |
---|
128 | *> \par Further Details: |
---|
129 | * ===================== |
---|
130 | *> |
---|
131 | *> \verbatim |
---|
132 | *> |
---|
133 | *> The factorization is obtained by Householder's method. The kth |
---|
134 | *> transformation matrix, Z( k ), which is used to introduce zeros into |
---|
135 | *> the ( m - k + 1 )th row of A, is given in the form |
---|
136 | *> |
---|
137 | *> Z( k ) = ( I 0 ), |
---|
138 | *> ( 0 T( k ) ) |
---|
139 | *> |
---|
140 | *> where |
---|
141 | *> |
---|
142 | *> T( k ) = I - tau*u( k )*u( k )**T, u( k ) = ( 1 ), |
---|
143 | *> ( 0 ) |
---|
144 | *> ( z( k ) ) |
---|
145 | *> |
---|
146 | *> tau is a scalar and z( k ) is an ( n - m ) element vector. |
---|
147 | *> tau and z( k ) are chosen to annihilate the elements of the kth row |
---|
148 | *> of X. |
---|
149 | *> |
---|
150 | *> The scalar tau is returned in the kth element of TAU and the vector |
---|
151 | *> u( k ) in the kth row of A, such that the elements of z( k ) are |
---|
152 | *> in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in |
---|
153 | *> the upper triangular part of A. |
---|
154 | *> |
---|
155 | *> Z is given by |
---|
156 | *> |
---|
157 | *> Z = Z( 1 ) * Z( 2 ) * ... * Z( m ). |
---|
158 | *> \endverbatim |
---|
159 | *> |
---|
160 | * ===================================================================== |
---|
161 | SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) |
---|
162 | * |
---|
163 | * -- LAPACK computational routine (version 3.4.0) -- |
---|
164 | * -- LAPACK is a software package provided by Univ. of Tennessee, -- |
---|
165 | * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
---|
166 | * November 2011 |
---|
167 | * |
---|
168 | * .. Scalar Arguments .. |
---|
169 | INTEGER INFO, LDA, LWORK, M, N |
---|
170 | * .. |
---|
171 | * .. Array Arguments .. |
---|
172 | REAL A( LDA, * ), TAU( * ), WORK( * ) |
---|
173 | * .. |
---|
174 | * |
---|
175 | * ===================================================================== |
---|
176 | * |
---|
177 | * .. Parameters .. |
---|
178 | REAL ZERO |
---|
179 | PARAMETER ( ZERO = 0.0E+0 ) |
---|
180 | * .. |
---|
181 | * .. Local Scalars .. |
---|
182 | LOGICAL LQUERY |
---|
183 | INTEGER I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT, |
---|
184 | $ M1, MU, NB, NBMIN, NX |
---|
185 | * .. |
---|
186 | * .. External Subroutines .. |
---|
187 | EXTERNAL XERBLA, SLARZB, SLARZT, SLATRZ |
---|
188 | * .. |
---|
189 | * .. Intrinsic Functions .. |
---|
190 | INTRINSIC MAX, MIN |
---|
191 | * .. |
---|
192 | * .. External Functions .. |
---|
193 | INTEGER ILAENV |
---|
194 | EXTERNAL ILAENV |
---|
195 | * .. |
---|
196 | * .. Executable Statements .. |
---|
197 | * |
---|
198 | * Test the input arguments |
---|
199 | * |
---|
200 | INFO = 0 |
---|
201 | LQUERY = ( LWORK.EQ.-1 ) |
---|
202 | IF( M.LT.0 ) THEN |
---|
203 | INFO = -1 |
---|
204 | ELSE IF( N.LT.M ) THEN |
---|
205 | INFO = -2 |
---|
206 | ELSE IF( LDA.LT.MAX( 1, M ) ) THEN |
---|
207 | INFO = -4 |
---|
208 | END IF |
---|
209 | * |
---|
210 | IF( INFO.EQ.0 ) THEN |
---|
211 | IF( M.EQ.0 .OR. M.EQ.N ) THEN |
---|
212 | LWKOPT = 1 |
---|
213 | LWKMIN = 1 |
---|
214 | ELSE |
---|
215 | * |
---|
216 | * Determine the block size. |
---|
217 | * |
---|
218 | NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 ) |
---|
219 | LWKOPT = M*NB |
---|
220 | LWKMIN = MAX( 1, M ) |
---|
221 | END IF |
---|
222 | WORK( 1 ) = LWKOPT |
---|
223 | * |
---|
224 | IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN |
---|
225 | INFO = -7 |
---|
226 | END IF |
---|
227 | END IF |
---|
228 | * |
---|
229 | IF( INFO.NE.0 ) THEN |
---|
230 | CALL XERBLA( 'STZRZF', -INFO ) |
---|
231 | RETURN |
---|
232 | ELSE IF( LQUERY ) THEN |
---|
233 | RETURN |
---|
234 | END IF |
---|
235 | * |
---|
236 | * Quick return if possible |
---|
237 | * |
---|
238 | IF( M.EQ.0 ) THEN |
---|
239 | RETURN |
---|
240 | ELSE IF( M.EQ.N ) THEN |
---|
241 | DO 10 I = 1, N |
---|
242 | TAU( I ) = ZERO |
---|
243 | 10 CONTINUE |
---|
244 | RETURN |
---|
245 | END IF |
---|
246 | * |
---|
247 | NBMIN = 2 |
---|
248 | NX = 1 |
---|
249 | IWS = M |
---|
250 | IF( NB.GT.1 .AND. NB.LT.M ) THEN |
---|
251 | * |
---|
252 | * Determine when to cross over from blocked to unblocked code. |
---|
253 | * |
---|
254 | NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) ) |
---|
255 | IF( NX.LT.M ) THEN |
---|
256 | * |
---|
257 | * Determine if workspace is large enough for blocked code. |
---|
258 | * |
---|
259 | LDWORK = M |
---|
260 | IWS = LDWORK*NB |
---|
261 | IF( LWORK.LT.IWS ) THEN |
---|
262 | * |
---|
263 | * Not enough workspace to use optimal NB: reduce NB and |
---|
264 | * determine the minimum value of NB. |
---|
265 | * |
---|
266 | NB = LWORK / LDWORK |
---|
267 | NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1, |
---|
268 | $ -1 ) ) |
---|
269 | END IF |
---|
270 | END IF |
---|
271 | END IF |
---|
272 | * |
---|
273 | IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN |
---|
274 | * |
---|
275 | * Use blocked code initially. |
---|
276 | * The last kk rows are handled by the block method. |
---|
277 | * |
---|
278 | M1 = MIN( M+1, N ) |
---|
279 | KI = ( ( M-NX-1 ) / NB )*NB |
---|
280 | KK = MIN( M, KI+NB ) |
---|
281 | * |
---|
282 | DO 20 I = M - KK + KI + 1, M - KK + 1, -NB |
---|
283 | IB = MIN( M-I+1, NB ) |
---|
284 | * |
---|
285 | * Compute the TZ factorization of the current block |
---|
286 | * A(i:i+ib-1,i:n) |
---|
287 | * |
---|
288 | CALL SLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ), |
---|
289 | $ WORK ) |
---|
290 | IF( I.GT.1 ) THEN |
---|
291 | * |
---|
292 | * Form the triangular factor of the block reflector |
---|
293 | * H = H(i+ib-1) . . . H(i+1) H(i) |
---|
294 | * |
---|
295 | CALL SLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ), |
---|
296 | $ LDA, TAU( I ), WORK, LDWORK ) |
---|
297 | * |
---|
298 | * Apply H to A(1:i-1,i:n) from the right |
---|
299 | * |
---|
300 | CALL SLARZB( 'Right', 'No transpose', 'Backward', |
---|
301 | $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ), |
---|
302 | $ LDA, WORK, LDWORK, A( 1, I ), LDA, |
---|
303 | $ WORK( IB+1 ), LDWORK ) |
---|
304 | END IF |
---|
305 | 20 CONTINUE |
---|
306 | MU = I + NB - 1 |
---|
307 | ELSE |
---|
308 | MU = M |
---|
309 | END IF |
---|
310 | * |
---|
311 | * Use unblocked code to factor the last or only block |
---|
312 | * |
---|
313 | IF( MU.GT.0 ) |
---|
314 | $ CALL SLATRZ( MU, N, N-M, A, LDA, TAU, WORK ) |
---|
315 | * |
---|
316 | WORK( 1 ) = LWKOPT |
---|
317 | * |
---|
318 | RETURN |
---|
319 | * |
---|
320 | * End of STZRZF |
---|
321 | * |
---|
322 | END |
---|