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