1: *> \brief \b ZTPLQT2 computes a LQ factorization of a real or complex "triangular-pentagonal" matrix, which is composed of a triangular block and a pentagonal block, using the compact WY representation for Q.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTPLQT2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztplqt2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztplqt2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztplqt2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTPLQT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, LDB, LDT, N, M, L
25: * ..
26: * .. Array Arguments ..
27: * COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> ZTPLQT2 computes a LQ a factorization of a complex "triangular-pentagonal"
37: *> matrix C, which is composed of a triangular block A and pentagonal block B,
38: *> using the compact WY representation for Q.
39: *> \endverbatim
40: *
41: * Arguments:
42: * ==========
43: *
44: *> \param[in] M
45: *> \verbatim
46: *> M is INTEGER
47: *> The total number of rows of the matrix B.
48: *> M >= 0.
49: *> \endverbatim
50: *>
51: *> \param[in] N
52: *> \verbatim
53: *> N is INTEGER
54: *> The number of columns of the matrix B, and the order of
55: *> the triangular matrix A.
56: *> N >= 0.
57: *> \endverbatim
58: *>
59: *> \param[in] L
60: *> \verbatim
61: *> L is INTEGER
62: *> The number of rows of the lower trapezoidal part of B.
63: *> MIN(M,N) >= L >= 0. See Further Details.
64: *> \endverbatim
65: *>
66: *> \param[in,out] A
67: *> \verbatim
68: *> A is COMPLEX*16 array, dimension (LDA,M)
69: *> On entry, the lower triangular M-by-M matrix A.
70: *> On exit, the elements on and below the diagonal of the array
71: *> contain the lower triangular matrix L.
72: *> \endverbatim
73: *>
74: *> \param[in] LDA
75: *> \verbatim
76: *> LDA is INTEGER
77: *> The leading dimension of the array A. LDA >= max(1,M).
78: *> \endverbatim
79: *>
80: *> \param[in,out] B
81: *> \verbatim
82: *> B is COMPLEX*16 array, dimension (LDB,N)
83: *> On entry, the pentagonal M-by-N matrix B. The first N-L columns
84: *> are rectangular, and the last L columns are lower trapezoidal.
85: *> On exit, B contains the pentagonal matrix V. See Further Details.
86: *> \endverbatim
87: *>
88: *> \param[in] LDB
89: *> \verbatim
90: *> LDB is INTEGER
91: *> The leading dimension of the array B. LDB >= max(1,M).
92: *> \endverbatim
93: *>
94: *> \param[out] T
95: *> \verbatim
96: *> T is COMPLEX*16 array, dimension (LDT,M)
97: *> The N-by-N upper triangular factor T of the block reflector.
98: *> See Further Details.
99: *> \endverbatim
100: *>
101: *> \param[in] LDT
102: *> \verbatim
103: *> LDT is INTEGER
104: *> The leading dimension of the array T. LDT >= max(1,M)
105: *> \endverbatim
106: *>
107: *> \param[out] INFO
108: *> \verbatim
109: *> INFO is INTEGER
110: *> = 0: successful exit
111: *> < 0: if INFO = -i, the i-th argument had an illegal value
112: *> \endverbatim
113: *
114: * Authors:
115: * ========
116: *
117: *> \author Univ. of Tennessee
118: *> \author Univ. of California Berkeley
119: *> \author Univ. of Colorado Denver
120: *> \author NAG Ltd.
121: *
122: *> \ingroup doubleOTHERcomputational
123: *
124: *> \par Further Details:
125: * =====================
126: *>
127: *> \verbatim
128: *>
129: *> The input matrix C is a M-by-(M+N) matrix
130: *>
131: *> C = [ A ][ B ]
132: *>
133: *>
134: *> where A is an lower triangular M-by-M matrix, and B is M-by-N pentagonal
135: *> matrix consisting of a M-by-(N-L) rectangular matrix B1 left of a M-by-L
136: *> upper trapezoidal matrix B2:
137: *>
138: *> B = [ B1 ][ B2 ]
139: *> [ B1 ] <- M-by-(N-L) rectangular
140: *> [ B2 ] <- M-by-L lower trapezoidal.
141: *>
142: *> The lower trapezoidal matrix B2 consists of the first L columns of a
143: *> N-by-N lower triangular matrix, where 0 <= L <= MIN(M,N). If L=0,
144: *> B is rectangular M-by-N; if M=L=N, B is lower triangular.
145: *>
146: *> The matrix W stores the elementary reflectors H(i) in the i-th row
147: *> above the diagonal (of A) in the M-by-(M+N) input matrix C
148: *>
149: *> C = [ A ][ B ]
150: *> [ A ] <- lower triangular M-by-M
151: *> [ B ] <- M-by-N pentagonal
152: *>
153: *> so that W can be represented as
154: *>
155: *> W = [ I ][ V ]
156: *> [ I ] <- identity, M-by-M
157: *> [ V ] <- M-by-N, same form as B.
158: *>
159: *> Thus, all of information needed for W is contained on exit in B, which
160: *> we call V above. Note that V has the same form as B; that is,
161: *>
162: *> W = [ V1 ][ V2 ]
163: *> [ V1 ] <- M-by-(N-L) rectangular
164: *> [ V2 ] <- M-by-L lower trapezoidal.
165: *>
166: *> The rows of V represent the vectors which define the H(i)'s.
167: *> The (M+N)-by-(M+N) block reflector H is then given by
168: *>
169: *> H = I - W**T * T * W
170: *>
171: *> where W^H is the conjugate transpose of W and T is the upper triangular
172: *> factor of the block reflector.
173: *> \endverbatim
174: *>
175: * =====================================================================
176: SUBROUTINE ZTPLQT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
177: *
178: * -- LAPACK computational routine --
179: * -- LAPACK is a software package provided by Univ. of Tennessee, --
180: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181: *
182: * .. Scalar Arguments ..
183: INTEGER INFO, LDA, LDB, LDT, N, M, L
184: * ..
185: * .. Array Arguments ..
186: COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * )
187: * ..
188: *
189: * =====================================================================
190: *
191: * .. Parameters ..
192: COMPLEX*16 ONE, ZERO
193: PARAMETER( ZERO = ( 0.0D+0, 0.0D+0 ),ONE = ( 1.0D+0, 0.0D+0 ) )
194: * ..
195: * .. Local Scalars ..
196: INTEGER I, J, P, MP, NP
197: COMPLEX*16 ALPHA
198: * ..
199: * .. External Subroutines ..
200: EXTERNAL ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA
201: * ..
202: * .. Intrinsic Functions ..
203: INTRINSIC MAX, MIN
204: * ..
205: * .. Executable Statements ..
206: *
207: * Test the input arguments
208: *
209: INFO = 0
210: IF( M.LT.0 ) THEN
211: INFO = -1
212: ELSE IF( N.LT.0 ) THEN
213: INFO = -2
214: ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
215: INFO = -3
216: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
217: INFO = -5
218: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
219: INFO = -7
220: ELSE IF( LDT.LT.MAX( 1, M ) ) THEN
221: INFO = -9
222: END IF
223: IF( INFO.NE.0 ) THEN
224: CALL XERBLA( 'ZTPLQT2', -INFO )
225: RETURN
226: END IF
227: *
228: * Quick return if possible
229: *
230: IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
231: *
232: DO I = 1, M
233: *
234: * Generate elementary reflector H(I) to annihilate B(I,:)
235: *
236: P = N-L+MIN( L, I )
237: CALL ZLARFG( P+1, A( I, I ), B( I, 1 ), LDB, T( 1, I ) )
238: T(1,I)=CONJG(T(1,I))
239: IF( I.LT.M ) THEN
240: DO J = 1, P
241: B( I, J ) = CONJG(B(I,J))
242: END DO
243: *
244: * W(M-I:1) := C(I+1:M,I:N) * C(I,I:N) [use W = T(M,:)]
245: *
246: DO J = 1, M-I
247: T( M, J ) = (A( I+J, I ))
248: END DO
249: CALL ZGEMV( 'N', M-I, P, ONE, B( I+1, 1 ), LDB,
250: $ B( I, 1 ), LDB, ONE, T( M, 1 ), LDT )
251: *
252: * C(I+1:M,I:N) = C(I+1:M,I:N) + alpha * C(I,I:N)*W(M-1:1)^H
253: *
254: ALPHA = -(T( 1, I ))
255: DO J = 1, M-I
256: A( I+J, I ) = A( I+J, I ) + ALPHA*(T( M, J ))
257: END DO
258: CALL ZGERC( M-I, P, (ALPHA), T( M, 1 ), LDT,
259: $ B( I, 1 ), LDB, B( I+1, 1 ), LDB )
260: DO J = 1, P
261: B( I, J ) = CONJG(B(I,J))
262: END DO
263: END IF
264: END DO
265: *
266: DO I = 2, M
267: *
268: * T(I,1:I-1) := C(I:I-1,1:N)**H * (alpha * C(I,I:N))
269: *
270: ALPHA = -(T( 1, I ))
271: DO J = 1, I-1
272: T( I, J ) = ZERO
273: END DO
274: P = MIN( I-1, L )
275: NP = MIN( N-L+1, N )
276: MP = MIN( P+1, M )
277: DO J = 1, N-L+P
278: B(I,J)=CONJG(B(I,J))
279: END DO
280: *
281: * Triangular part of B2
282: *
283: DO J = 1, P
284: T( I, J ) = (ALPHA*B( I, N-L+J ))
285: END DO
286: CALL ZTRMV( 'L', 'N', 'N', P, B( 1, NP ), LDB,
287: $ T( I, 1 ), LDT )
288: *
289: * Rectangular part of B2
290: *
291: CALL ZGEMV( 'N', I-1-P, L, ALPHA, B( MP, NP ), LDB,
292: $ B( I, NP ), LDB, ZERO, T( I,MP ), LDT )
293: *
294: * B1
295:
296: *
297: CALL ZGEMV( 'N', I-1, N-L, ALPHA, B, LDB, B( I, 1 ), LDB,
298: $ ONE, T( I, 1 ), LDT )
299: *
300:
301: *
302: * T(1:I-1,I) := T(1:I-1,1:I-1) * T(I,1:I-1)
303: *
304: DO J = 1, I-1
305: T(I,J)=CONJG(T(I,J))
306: END DO
307: CALL ZTRMV( 'L', 'C', 'N', I-1, T, LDT, T( I, 1 ), LDT )
308: DO J = 1, I-1
309: T(I,J)=CONJG(T(I,J))
310: END DO
311: DO J = 1, N-L+P
312: B(I,J)=CONJG(B(I,J))
313: END DO
314: *
315: * T(I,I) = tau(I)
316: *
317: T( I, I ) = T( 1, I )
318: T( 1, I ) = ZERO
319: END DO
320: DO I=1,M
321: DO J= I+1,M
322: T(I,J)=(T(J,I))
323: T(J,I)=ZERO
324: END DO
325: END DO
326:
327: *
328: * End of ZTPLQT2
329: *
330: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztplqt2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack