1: *> \brief \b ZGELQT3 recursively computes a LQ factorization of a general real or complex matrix using the compact WY representation of Q.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DGEQRT3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelqt3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelqt3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelqt3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * RECURSIVE SUBROUTINE ZGELQT3( M, N, A, LDA, T, LDT, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, M, N, LDT
25: * ..
26: * .. Array Arguments ..
27: * COMPLEX*16 A( LDA, * ), T( LDT, * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DGELQT3 recursively computes a LQ factorization of a complex M-by-N
37: *> matrix A, using the compact WY representation of Q.
38: *>
39: *> Based on the algorithm of Elmroth and Gustavson,
40: *> IBM J. Res. Develop. Vol 44 No. 4 July 2000.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] M
47: *> \verbatim
48: *> M is INTEGER
49: *> The number of rows of the matrix A. M =< N.
50: *> \endverbatim
51: *>
52: *> \param[in] N
53: *> \verbatim
54: *> N is INTEGER
55: *> The number of columns of the matrix A. N >= 0.
56: *> \endverbatim
57: *>
58: *> \param[in,out] A
59: *> \verbatim
60: *> A is COMPLEX*16 array, dimension (LDA,N)
61: *> On entry, the real M-by-N matrix A. On exit, the elements on and
62: *> below the diagonal contain the N-by-N lower triangular matrix L; the
63: *> elements above the diagonal are the rows of V. See below for
64: *> further details.
65: *> \endverbatim
66: *>
67: *> \param[in] LDA
68: *> \verbatim
69: *> LDA is INTEGER
70: *> The leading dimension of the array A. LDA >= max(1,M).
71: *> \endverbatim
72: *>
73: *> \param[out] T
74: *> \verbatim
75: *> T is COMPLEX*16 array, dimension (LDT,N)
76: *> The N-by-N upper triangular factor of the block reflector.
77: *> The elements on and above the diagonal contain the block
78: *> reflector T; the elements below the diagonal are not used.
79: *> See below for further details.
80: *> \endverbatim
81: *>
82: *> \param[in] LDT
83: *> \verbatim
84: *> LDT is INTEGER
85: *> The leading dimension of the array T. LDT >= max(1,N).
86: *> \endverbatim
87: *>
88: *> \param[out] INFO
89: *> \verbatim
90: *> INFO is INTEGER
91: *> = 0: successful exit
92: *> < 0: if INFO = -i, the i-th argument had an illegal value
93: *> \endverbatim
94: *
95: * Authors:
96: * ========
97: *
98: *> \author Univ. of Tennessee
99: *> \author Univ. of California Berkeley
100: *> \author Univ. of Colorado Denver
101: *> \author NAG Ltd.
102: *
103: *> \date December 2016
104: *
105: *> \ingroup doubleGEcomputational
106: *
107: *> \par Further Details:
108: * =====================
109: *>
110: *> \verbatim
111: *>
112: *> The matrix V stores the elementary reflectors H(i) in the i-th column
113: *> below the diagonal. For example, if M=5 and N=3, the matrix V is
114: *>
115: *> V = ( 1 v1 v1 v1 v1 )
116: *> ( 1 v2 v2 v2 )
117: *> ( 1 v3 v3 v3 )
118: *>
119: *>
120: *> where the vi's represent the vectors which define H(i), which are returned
121: *> in the matrix A. The 1's along the diagonal of V are not stored in A. The
122: *> block reflector H is then given by
123: *>
124: *> H = I - V * T * V**T
125: *>
126: *> where V**T is the transpose of V.
127: *>
128: *> For details of the algorithm, see Elmroth and Gustavson (cited above).
129: *> \endverbatim
130: *>
131: * =====================================================================
132: RECURSIVE SUBROUTINE ZGELQT3( M, N, A, LDA, T, LDT, INFO )
133: *
134: * -- LAPACK computational routine (version 3.7.0) --
135: * -- LAPACK is a software package provided by Univ. of Tennessee, --
136: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
137: * December 2016
138: *
139: * .. Scalar Arguments ..
140: INTEGER INFO, LDA, M, N, LDT
141: * ..
142: * .. Array Arguments ..
143: COMPLEX*16 A( LDA, * ), T( LDT, * )
144: * ..
145: *
146: * =====================================================================
147: *
148: * .. Parameters ..
149: COMPLEX*16 ONE, ZERO
150: PARAMETER ( ONE = (1.0D+00,0.0D+00) )
151: PARAMETER ( ZERO = (0.0D+00,0.0D+00))
152: * ..
153: * .. Local Scalars ..
154: INTEGER I, I1, J, J1, M1, M2, N1, N2, IINFO
155: * ..
156: * .. External Subroutines ..
157: EXTERNAL ZLARFG, ZTRMM, ZGEMM, XERBLA
158: * ..
159: * .. Executable Statements ..
160: *
161: INFO = 0
162: IF( M .LT. 0 ) THEN
163: INFO = -1
164: ELSE IF( N .LT. M ) THEN
165: INFO = -2
166: ELSE IF( LDA .LT. MAX( 1, M ) ) THEN
167: INFO = -4
168: ELSE IF( LDT .LT. MAX( 1, M ) ) THEN
169: INFO = -6
170: END IF
171: IF( INFO.NE.0 ) THEN
172: CALL XERBLA( 'ZGELQT3', -INFO )
173: RETURN
174: END IF
175: *
176: IF( M.EQ.1 ) THEN
177: *
178: * Compute Householder transform when N=1
179: *
180: CALL ZLARFG( N, A, A( 1, MIN( 2, N ) ), LDA, T )
181: T(1,1)=CONJG(T(1,1))
182: *
183: ELSE
184: *
185: * Otherwise, split A into blocks...
186: *
187: M1 = M/2
188: M2 = M-M1
189: I1 = MIN( M1+1, M )
190: J1 = MIN( M+1, N )
191: *
192: * Compute A(1:M1,1:N) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1^H
193: *
194: CALL ZGELQT3( M1, N, A, LDA, T, LDT, IINFO )
195: *
196: * Compute A(J1:M,1:N) = A(J1:M,1:N) Q1^H [workspace: T(1:N1,J1:N)]
197: *
198: DO I=1,M2
199: DO J=1,M1
200: T( I+M1, J ) = A( I+M1, J )
201: END DO
202: END DO
203: CALL ZTRMM( 'R', 'U', 'C', 'U', M2, M1, ONE,
204: & A, LDA, T( I1, 1 ), LDT )
205: *
206: CALL ZGEMM( 'N', 'C', M2, M1, N-M1, ONE, A( I1, I1 ), LDA,
207: & A( 1, I1 ), LDA, ONE, T( I1, 1 ), LDT)
208: *
209: CALL ZTRMM( 'R', 'U', 'N', 'N', M2, M1, ONE,
210: & T, LDT, T( I1, 1 ), LDT )
211: *
212: CALL ZGEMM( 'N', 'N', M2, N-M1, M1, -ONE, T( I1, 1 ), LDT,
213: & A( 1, I1 ), LDA, ONE, A( I1, I1 ), LDA )
214: *
215: CALL ZTRMM( 'R', 'U', 'N', 'U', M2, M1 , ONE,
216: & A, LDA, T( I1, 1 ), LDT )
217: *
218: DO I=1,M2
219: DO J=1,M1
220: A( I+M1, J ) = A( I+M1, J ) - T( I+M1, J )
221: T( I+M1, J )= ZERO
222: END DO
223: END DO
224: *
225: * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2^H
226: *
227: CALL ZGELQT3( M2, N-M1, A( I1, I1 ), LDA,
228: & T( I1, I1 ), LDT, IINFO )
229: *
230: * Compute T3 = T(J1:N1,1:N) = -T1 Y1^H Y2 T2
231: *
232: DO I=1,M2
233: DO J=1,M1
234: T( J, I+M1 ) = (A( J, I+M1 ))
235: END DO
236: END DO
237: *
238: CALL ZTRMM( 'R', 'U', 'C', 'U', M1, M2, ONE,
239: & A( I1, I1 ), LDA, T( 1, I1 ), LDT )
240: *
241: CALL ZGEMM( 'N', 'C', M1, M2, N-M, ONE, A( 1, J1 ), LDA,
242: & A( I1, J1 ), LDA, ONE, T( 1, I1 ), LDT )
243: *
244: CALL ZTRMM( 'L', 'U', 'N', 'N', M1, M2, -ONE, T, LDT,
245: & T( 1, I1 ), LDT )
246: *
247: CALL ZTRMM( 'R', 'U', 'N', 'N', M1, M2, ONE,
248: & T( I1, I1 ), LDT, T( 1, I1 ), LDT )
249: *
250: *
251: *
252: * Y = (Y1,Y2); L = [ L1 0 ]; T = [T1 T3]
253: * [ A(1:N1,J1:N) L2 ] [ 0 T2]
254: *
255: END IF
256: *
257: RETURN
258: *
259: * End of ZGELQT3
260: *
261: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelqt3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack