Annotation of rpl/lapack/lapack/ztzrzf.f, revision 1.15
1.9 bertrand 1: *> \brief \b ZTZRZF
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTZRZF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztzrzf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztzrzf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztzrzf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, LWORK, M, N
25: * ..
26: * .. Array Arguments ..
27: * COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> ZTZRZF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A
37: *> to upper triangular form by means of unitary 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 unitary 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 COMPLEX*16 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: *> unitary 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 COMPLEX*16 array, dimension (M)
82: *> The scalar factors of the elementary reflectors.
83: *> \endverbatim
84: *>
85: *> \param[out] WORK
86: *> \verbatim
87: *> WORK is COMPLEX*16 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: *
1.11 bertrand 119: *> \date April 2012
1.9 bertrand 120: *
121: *> \ingroup complex16OTHERcomputational
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: *>
1.11 bertrand 133: *> The N-by-N matrix Z can be computed by
1.9 bertrand 134: *>
1.11 bertrand 135: *> Z = Z(1)*Z(2)* ... *Z(M)
1.9 bertrand 136: *>
1.11 bertrand 137: *> where each N-by-N Z(k) is given by
1.9 bertrand 138: *>
1.11 bertrand 139: *> Z(k) = I - tau(k)*v(k)*v(k)**H
1.9 bertrand 140: *>
1.11 bertrand 141: *> with v(k) is the kth row vector of the M-by-N matrix
1.9 bertrand 142: *>
1.11 bertrand 143: *> V = ( I A(:,M+1:N) )
1.9 bertrand 144: *>
1.11 bertrand 145: *> I is the M-by-M identity matrix, A(:,M+1:N)
146: *> is the output stored in A on exit from DTZRZF,
147: *> and tau(k) is the kth element of the array TAU.
1.9 bertrand 148: *>
149: *> \endverbatim
150: *>
151: * =====================================================================
1.1 bertrand 152: SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
153: *
1.11 bertrand 154: * -- LAPACK computational routine (version 3.4.1) --
1.1 bertrand 155: * -- LAPACK is a software package provided by Univ. of Tennessee, --
156: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 157: * April 2012
1.1 bertrand 158: *
159: * .. Scalar Arguments ..
160: INTEGER INFO, LDA, LWORK, M, N
161: * ..
162: * .. Array Arguments ..
163: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
164: * ..
165: *
166: * =====================================================================
167: *
168: * .. Parameters ..
169: COMPLEX*16 ZERO
170: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
171: * ..
172: * .. Local Scalars ..
173: LOGICAL LQUERY
1.8 bertrand 174: INTEGER I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
175: $ M1, MU, NB, NBMIN, NX
1.1 bertrand 176: * ..
177: * .. External Subroutines ..
178: EXTERNAL XERBLA, ZLARZB, ZLARZT, ZLATRZ
179: * ..
180: * .. Intrinsic Functions ..
181: INTRINSIC MAX, MIN
182: * ..
183: * .. External Functions ..
184: INTEGER ILAENV
185: EXTERNAL ILAENV
186: * ..
187: * .. Executable Statements ..
188: *
189: * Test the input arguments
190: *
191: INFO = 0
192: LQUERY = ( LWORK.EQ.-1 )
193: IF( M.LT.0 ) THEN
194: INFO = -1
195: ELSE IF( N.LT.M ) THEN
196: INFO = -2
197: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
198: INFO = -4
199: END IF
200: *
201: IF( INFO.EQ.0 ) THEN
202: IF( M.EQ.0 .OR. M.EQ.N ) THEN
203: LWKOPT = 1
1.8 bertrand 204: LWKMIN = 1
1.1 bertrand 205: ELSE
206: *
207: * Determine the block size.
208: *
209: NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
210: LWKOPT = M*NB
1.8 bertrand 211: LWKMIN = MAX( 1, M )
1.1 bertrand 212: END IF
213: WORK( 1 ) = LWKOPT
214: *
1.8 bertrand 215: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
1.1 bertrand 216: INFO = -7
217: END IF
218: END IF
219: *
220: IF( INFO.NE.0 ) THEN
221: CALL XERBLA( 'ZTZRZF', -INFO )
222: RETURN
223: ELSE IF( LQUERY ) THEN
224: RETURN
225: END IF
226: *
227: * Quick return if possible
228: *
229: IF( M.EQ.0 ) THEN
230: RETURN
231: ELSE IF( M.EQ.N ) THEN
232: DO 10 I = 1, N
233: TAU( I ) = ZERO
234: 10 CONTINUE
235: RETURN
236: END IF
237: *
238: NBMIN = 2
239: NX = 1
240: IWS = M
241: IF( NB.GT.1 .AND. NB.LT.M ) THEN
242: *
243: * Determine when to cross over from blocked to unblocked code.
244: *
245: NX = MAX( 0, ILAENV( 3, 'ZGERQF', ' ', M, N, -1, -1 ) )
246: IF( NX.LT.M ) THEN
247: *
248: * Determine if workspace is large enough for blocked code.
249: *
250: LDWORK = M
251: IWS = LDWORK*NB
252: IF( LWORK.LT.IWS ) THEN
253: *
254: * Not enough workspace to use optimal NB: reduce NB and
255: * determine the minimum value of NB.
256: *
257: NB = LWORK / LDWORK
258: NBMIN = MAX( 2, ILAENV( 2, 'ZGERQF', ' ', M, N, -1,
259: $ -1 ) )
260: END IF
261: END IF
262: END IF
263: *
264: IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
265: *
266: * Use blocked code initially.
267: * The last kk rows are handled by the block method.
268: *
269: M1 = MIN( M+1, N )
270: KI = ( ( M-NX-1 ) / NB )*NB
271: KK = MIN( M, KI+NB )
272: *
273: DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
274: IB = MIN( M-I+1, NB )
275: *
276: * Compute the TZ factorization of the current block
277: * A(i:i+ib-1,i:n)
278: *
279: CALL ZLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
280: $ WORK )
281: IF( I.GT.1 ) THEN
282: *
283: * Form the triangular factor of the block reflector
284: * H = H(i+ib-1) . . . H(i+1) H(i)
285: *
286: CALL ZLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
287: $ LDA, TAU( I ), WORK, LDWORK )
288: *
289: * Apply H to A(1:i-1,i:n) from the right
290: *
291: CALL ZLARZB( 'Right', 'No transpose', 'Backward',
292: $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
293: $ LDA, WORK, LDWORK, A( 1, I ), LDA,
294: $ WORK( IB+1 ), LDWORK )
295: END IF
296: 20 CONTINUE
297: MU = I + NB - 1
298: ELSE
299: MU = M
300: END IF
301: *
302: * Use unblocked code to factor the last or only block
303: *
304: IF( MU.GT.0 )
305: $ CALL ZLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
306: *
307: WORK( 1 ) = LWKOPT
308: *
309: RETURN
310: *
311: * End of ZTZRZF
312: *
313: END
CVSweb interface <joel.bertrand@systella.fr>