1: *> \brief \b ZGEQRT2
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGEQRT2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeqrt2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeqrt2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeqrt2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGEQRT2( M, N, A, LDA, T, LDT, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, LDT, M, N
25: * ..
26: * .. Array Arguments ..
27: * COMPLEX*16 A( LDA, * ), T( LDT, * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> ZGEQRT2 computes a QR factorization of a complex M-by-N matrix A,
37: *> using the compact WY representation of Q.
38: *> \endverbatim
39: *
40: * Arguments:
41: * ==========
42: *
43: *> \param[in] M
44: *> \verbatim
45: *> M is INTEGER
46: *> The number of rows of the matrix A. M >= N.
47: *> \endverbatim
48: *>
49: *> \param[in] N
50: *> \verbatim
51: *> N is INTEGER
52: *> The number of columns of the matrix A. N >= 0.
53: *> \endverbatim
54: *>
55: *> \param[in,out] A
56: *> \verbatim
57: *> A is COMPLEX*16 array, dimension (LDA,N)
58: *> On entry, the complex M-by-N matrix A. On exit, the elements on and
59: *> above the diagonal contain the N-by-N upper triangular matrix R; the
60: *> elements below the diagonal are the columns of V. See below for
61: *> further details.
62: *> \endverbatim
63: *>
64: *> \param[in] LDA
65: *> \verbatim
66: *> LDA is INTEGER
67: *> The leading dimension of the array A. LDA >= max(1,M).
68: *> \endverbatim
69: *>
70: *> \param[out] T
71: *> \verbatim
72: *> T is COMPLEX*16 array, dimension (LDT,N)
73: *> The N-by-N upper triangular factor of the block reflector.
74: *> The elements on and above the diagonal contain the block
75: *> reflector T; the elements below the diagonal are not used.
76: *> See below for further details.
77: *> \endverbatim
78: *>
79: *> \param[in] LDT
80: *> \verbatim
81: *> LDT is INTEGER
82: *> The leading dimension of the array T. LDT >= max(1,N).
83: *> \endverbatim
84: *>
85: *> \param[out] INFO
86: *> \verbatim
87: *> INFO is INTEGER
88: *> = 0: successful exit
89: *> < 0: if INFO = -i, the i-th argument had an illegal value
90: *> \endverbatim
91: *
92: * Authors:
93: * ========
94: *
95: *> \author Univ. of Tennessee
96: *> \author Univ. of California Berkeley
97: *> \author Univ. of Colorado Denver
98: *> \author NAG Ltd.
99: *
100: *> \date November 2011
101: *
102: *> \ingroup complex16GEcomputational
103: *
104: *> \par Further Details:
105: * =====================
106: *>
107: *> \verbatim
108: *>
109: *> The matrix V stores the elementary reflectors H(i) in the i-th column
110: *> below the diagonal. For example, if M=5 and N=3, the matrix V is
111: *>
112: *> V = ( 1 )
113: *> ( v1 1 )
114: *> ( v1 v2 1 )
115: *> ( v1 v2 v3 )
116: *> ( v1 v2 v3 )
117: *>
118: *> where the vi's represent the vectors which define H(i), which are returned
119: *> in the matrix A. The 1's along the diagonal of V are not stored in A. The
120: *> block reflector H is then given by
121: *>
122: *> H = I - V * T * V**H
123: *>
124: *> where V**H is the conjugate transpose of V.
125: *> \endverbatim
126: *>
127: * =====================================================================
128: SUBROUTINE ZGEQRT2( M, N, A, LDA, T, LDT, INFO )
129: *
130: * -- LAPACK computational routine (version 3.4.0) --
131: * -- LAPACK is a software package provided by Univ. of Tennessee, --
132: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133: * November 2011
134: *
135: * .. Scalar Arguments ..
136: INTEGER INFO, LDA, LDT, M, N
137: * ..
138: * .. Array Arguments ..
139: COMPLEX*16 A( LDA, * ), T( LDT, * )
140: * ..
141: *
142: * =====================================================================
143: *
144: * .. Parameters ..
145: COMPLEX*16 ONE, ZERO
146: PARAMETER( ONE = (1.0D+00,0.0D+00), ZERO = (0.0D+00,0.0D+00) )
147: * ..
148: * .. Local Scalars ..
149: INTEGER I, K
150: COMPLEX*16 AII, ALPHA
151: * ..
152: * .. External Subroutines ..
153: EXTERNAL ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA
154: * ..
155: * .. Executable Statements ..
156: *
157: * Test the input arguments
158: *
159: INFO = 0
160: IF( M.LT.0 ) THEN
161: INFO = -1
162: ELSE IF( N.LT.0 ) THEN
163: INFO = -2
164: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
165: INFO = -4
166: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
167: INFO = -6
168: END IF
169: IF( INFO.NE.0 ) THEN
170: CALL XERBLA( 'ZGEQRT2', -INFO )
171: RETURN
172: END IF
173: *
174: K = MIN( M, N )
175: *
176: DO I = 1, K
177: *
178: * Generate elem. refl. H(i) to annihilate A(i+1:m,i), tau(I) -> T(I,1)
179: *
180: CALL ZLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
181: $ T( I, 1 ) )
182: IF( I.LT.N ) THEN
183: *
184: * Apply H(i) to A(I:M,I+1:N) from the left
185: *
186: AII = A( I, I )
187: A( I, I ) = ONE
188: *
189: * W(1:N-I) := A(I:M,I+1:N)^H * A(I:M,I) [W = T(:,N)]
190: *
191: CALL ZGEMV( 'C',M-I+1, N-I, ONE, A( I, I+1 ), LDA,
192: $ A( I, I ), 1, ZERO, T( 1, N ), 1 )
193: *
194: * A(I:M,I+1:N) = A(I:m,I+1:N) + alpha*A(I:M,I)*W(1:N-1)^H
195: *
196: ALPHA = -CONJG(T( I, 1 ))
197: CALL ZGERC( M-I+1, N-I, ALPHA, A( I, I ), 1,
198: $ T( 1, N ), 1, A( I, I+1 ), LDA )
199: A( I, I ) = AII
200: END IF
201: END DO
202: *
203: DO I = 2, N
204: AII = A( I, I )
205: A( I, I ) = ONE
206: *
207: * T(1:I-1,I) := alpha * A(I:M,1:I-1)**H * A(I:M,I)
208: *
209: ALPHA = -T( I, 1 )
210: CALL ZGEMV( 'C', M-I+1, I-1, ALPHA, A( I, 1 ), LDA,
211: $ A( I, I ), 1, ZERO, T( 1, I ), 1 )
212: A( I, I ) = AII
213: *
214: * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
215: *
216: CALL ZTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
217: *
218: * T(I,I) = tau(I)
219: *
220: T( I, I ) = T( I, 1 )
221: T( I, 1) = ZERO
222: END DO
223:
224: *
225: * End of ZGEQRT2
226: *
227: END
CVSweb interface <joel.bertrand@systella.fr>