Annotation of rpl/lapack/lapack/ztzrzf.f, revision 1.3
1.1 bertrand 1: SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, LDA, LWORK, M, N
10: * ..
11: * .. Array Arguments ..
12: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
13: * ..
14: *
15: * Purpose
16: * =======
17: *
18: * ZTZRZF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A
19: * to upper triangular form by means of unitary transformations.
20: *
21: * The upper trapezoidal matrix A is factored as
22: *
23: * A = ( R 0 ) * Z,
24: *
25: * where Z is an N-by-N unitary matrix and R is an M-by-M upper
26: * triangular matrix.
27: *
28: * Arguments
29: * =========
30: *
31: * M (input) INTEGER
32: * The number of rows of the matrix A. M >= 0.
33: *
34: * N (input) INTEGER
35: * The number of columns of the matrix A. N >= M.
36: *
37: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
38: * On entry, the leading M-by-N upper trapezoidal part of the
39: * array A must contain the matrix to be factorized.
40: * On exit, the leading M-by-M upper triangular part of A
41: * contains the upper triangular matrix R, and elements M+1 to
42: * N of the first M rows of A, with the array TAU, represent the
43: * unitary matrix Z as a product of M elementary reflectors.
44: *
45: * LDA (input) INTEGER
46: * The leading dimension of the array A. LDA >= max(1,M).
47: *
48: * TAU (output) COMPLEX*16 array, dimension (M)
49: * The scalar factors of the elementary reflectors.
50: *
51: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
52: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
53: *
54: * LWORK (input) INTEGER
55: * The dimension of the array WORK. LWORK >= max(1,M).
56: * For optimum performance LWORK >= M*NB, where NB is
57: * the optimal blocksize.
58: *
59: * If LWORK = -1, then a workspace query is assumed; the routine
60: * only calculates the optimal size of the WORK array, returns
61: * this value as the first entry of the WORK array, and no error
62: * message related to LWORK is issued by XERBLA.
63: *
64: * INFO (output) INTEGER
65: * = 0: successful exit
66: * < 0: if INFO = -i, the i-th argument had an illegal value
67: *
68: * Further Details
69: * ===============
70: *
71: * Based on contributions by
72: * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
73: *
74: * The factorization is obtained by Householder's method. The kth
75: * transformation matrix, Z( k ), which is used to introduce zeros into
76: * the ( m - k + 1 )th row of A, is given in the form
77: *
78: * Z( k ) = ( I 0 ),
79: * ( 0 T( k ) )
80: *
81: * where
82: *
83: * T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ),
84: * ( 0 )
85: * ( z( k ) )
86: *
87: * tau is a scalar and z( k ) is an ( n - m ) element vector.
88: * tau and z( k ) are chosen to annihilate the elements of the kth row
89: * of X.
90: *
91: * The scalar tau is returned in the kth element of TAU and the vector
92: * u( k ) in the kth row of A, such that the elements of z( k ) are
93: * in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
94: * the upper triangular part of A.
95: *
96: * Z is given by
97: *
98: * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ).
99: *
100: * =====================================================================
101: *
102: * .. Parameters ..
103: COMPLEX*16 ZERO
104: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
105: * ..
106: * .. Local Scalars ..
107: LOGICAL LQUERY
108: INTEGER I, IB, IWS, KI, KK, LDWORK, LWKOPT, M1, MU, NB,
109: $ NBMIN, NX
110: * ..
111: * .. External Subroutines ..
112: EXTERNAL XERBLA, ZLARZB, ZLARZT, ZLATRZ
113: * ..
114: * .. Intrinsic Functions ..
115: INTRINSIC MAX, MIN
116: * ..
117: * .. External Functions ..
118: INTEGER ILAENV
119: EXTERNAL ILAENV
120: * ..
121: * .. Executable Statements ..
122: *
123: * Test the input arguments
124: *
125: INFO = 0
126: LQUERY = ( LWORK.EQ.-1 )
127: IF( M.LT.0 ) THEN
128: INFO = -1
129: ELSE IF( N.LT.M ) THEN
130: INFO = -2
131: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
132: INFO = -4
133: END IF
134: *
135: IF( INFO.EQ.0 ) THEN
136: IF( M.EQ.0 .OR. M.EQ.N ) THEN
137: LWKOPT = 1
138: ELSE
139: *
140: * Determine the block size.
141: *
142: NB = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
143: LWKOPT = M*NB
144: END IF
145: WORK( 1 ) = LWKOPT
146: *
147: IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
148: INFO = -7
149: END IF
150: END IF
151: *
152: IF( INFO.NE.0 ) THEN
153: CALL XERBLA( 'ZTZRZF', -INFO )
154: RETURN
155: ELSE IF( LQUERY ) THEN
156: RETURN
157: END IF
158: *
159: * Quick return if possible
160: *
161: IF( M.EQ.0 ) THEN
162: RETURN
163: ELSE IF( M.EQ.N ) THEN
164: DO 10 I = 1, N
165: TAU( I ) = ZERO
166: 10 CONTINUE
167: RETURN
168: END IF
169: *
170: NBMIN = 2
171: NX = 1
172: IWS = M
173: IF( NB.GT.1 .AND. NB.LT.M ) THEN
174: *
175: * Determine when to cross over from blocked to unblocked code.
176: *
177: NX = MAX( 0, ILAENV( 3, 'ZGERQF', ' ', M, N, -1, -1 ) )
178: IF( NX.LT.M ) THEN
179: *
180: * Determine if workspace is large enough for blocked code.
181: *
182: LDWORK = M
183: IWS = LDWORK*NB
184: IF( LWORK.LT.IWS ) THEN
185: *
186: * Not enough workspace to use optimal NB: reduce NB and
187: * determine the minimum value of NB.
188: *
189: NB = LWORK / LDWORK
190: NBMIN = MAX( 2, ILAENV( 2, 'ZGERQF', ' ', M, N, -1,
191: $ -1 ) )
192: END IF
193: END IF
194: END IF
195: *
196: IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
197: *
198: * Use blocked code initially.
199: * The last kk rows are handled by the block method.
200: *
201: M1 = MIN( M+1, N )
202: KI = ( ( M-NX-1 ) / NB )*NB
203: KK = MIN( M, KI+NB )
204: *
205: DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
206: IB = MIN( M-I+1, NB )
207: *
208: * Compute the TZ factorization of the current block
209: * A(i:i+ib-1,i:n)
210: *
211: CALL ZLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
212: $ WORK )
213: IF( I.GT.1 ) THEN
214: *
215: * Form the triangular factor of the block reflector
216: * H = H(i+ib-1) . . . H(i+1) H(i)
217: *
218: CALL ZLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
219: $ LDA, TAU( I ), WORK, LDWORK )
220: *
221: * Apply H to A(1:i-1,i:n) from the right
222: *
223: CALL ZLARZB( 'Right', 'No transpose', 'Backward',
224: $ 'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
225: $ LDA, WORK, LDWORK, A( 1, I ), LDA,
226: $ WORK( IB+1 ), LDWORK )
227: END IF
228: 20 CONTINUE
229: MU = I + NB - 1
230: ELSE
231: MU = M
232: END IF
233: *
234: * Use unblocked code to factor the last or only block
235: *
236: IF( MU.GT.0 )
237: $ CALL ZLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
238: *
239: WORK( 1 ) = LWKOPT
240: *
241: RETURN
242: *
243: * End of ZTZRZF
244: *
245: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztzrzf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack