1: *> \brief \b DGEQRT2 computes a QR 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 DGEQRT2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrt2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrt2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrt2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEQRT2( M, N, A, LDA, T, LDT, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, LDT, M, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), T( LDT, * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DGEQRT2 computes a QR factorization of a real 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 DOUBLE PRECISION array, dimension (LDA,N)
58: *> On entry, the real 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 DOUBLE PRECISION 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: *> \ingroup doubleGEcomputational
101: *
102: *> \par Further Details:
103: * =====================
104: *>
105: *> \verbatim
106: *>
107: *> The matrix V stores the elementary reflectors H(i) in the i-th column
108: *> below the diagonal. For example, if M=5 and N=3, the matrix V is
109: *>
110: *> V = ( 1 )
111: *> ( v1 1 )
112: *> ( v1 v2 1 )
113: *> ( v1 v2 v3 )
114: *> ( v1 v2 v3 )
115: *>
116: *> where the vi's represent the vectors which define H(i), which are returned
117: *> in the matrix A. The 1's along the diagonal of V are not stored in A. The
118: *> block reflector H is then given by
119: *>
120: *> H = I - V * T * V**T
121: *>
122: *> where V**T is the transpose of V.
123: *> \endverbatim
124: *>
125: * =====================================================================
126: SUBROUTINE DGEQRT2( M, N, A, LDA, T, LDT, INFO )
127: *
128: * -- LAPACK computational routine --
129: * -- LAPACK is a software package provided by Univ. of Tennessee, --
130: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
131: *
132: * .. Scalar Arguments ..
133: INTEGER INFO, LDA, LDT, M, N
134: * ..
135: * .. Array Arguments ..
136: DOUBLE PRECISION A( LDA, * ), T( LDT, * )
137: * ..
138: *
139: * =====================================================================
140: *
141: * .. Parameters ..
142: DOUBLE PRECISION ONE, ZERO
143: PARAMETER( ONE = 1.0D+00, ZERO = 0.0D+00 )
144: * ..
145: * .. Local Scalars ..
146: INTEGER I, K
147: DOUBLE PRECISION AII, ALPHA
148: * ..
149: * .. External Subroutines ..
150: EXTERNAL DLARFG, DGEMV, DGER, DTRMV, XERBLA
151: * ..
152: * .. Executable Statements ..
153: *
154: * Test the input arguments
155: *
156: INFO = 0
157: IF( N.LT.0 ) THEN
158: INFO = -2
159: ELSE IF( M.LT.N ) THEN
160: INFO = -1
161: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
162: INFO = -4
163: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
164: INFO = -6
165: END IF
166: IF( INFO.NE.0 ) THEN
167: CALL XERBLA( 'DGEQRT2', -INFO )
168: RETURN
169: END IF
170: *
171: K = MIN( M, N )
172: *
173: DO I = 1, K
174: *
175: * Generate elem. refl. H(i) to annihilate A(i+1:m,i), tau(I) -> T(I,1)
176: *
177: CALL DLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
178: $ T( I, 1 ) )
179: IF( I.LT.N ) THEN
180: *
181: * Apply H(i) to A(I:M,I+1:N) from the left
182: *
183: AII = A( I, I )
184: A( I, I ) = ONE
185: *
186: * W(1:N-I) := A(I:M,I+1:N)^H * A(I:M,I) [W = T(:,N)]
187: *
188: CALL DGEMV( 'T',M-I+1, N-I, ONE, A( I, I+1 ), LDA,
189: $ A( I, I ), 1, ZERO, T( 1, N ), 1 )
190: *
191: * A(I:M,I+1:N) = A(I:m,I+1:N) + alpha*A(I:M,I)*W(1:N-1)^H
192: *
193: ALPHA = -(T( I, 1 ))
194: CALL DGER( M-I+1, N-I, ALPHA, A( I, I ), 1,
195: $ T( 1, N ), 1, A( I, I+1 ), LDA )
196: A( I, I ) = AII
197: END IF
198: END DO
199: *
200: DO I = 2, N
201: AII = A( I, I )
202: A( I, I ) = ONE
203: *
204: * T(1:I-1,I) := alpha * A(I:M,1:I-1)**T * A(I:M,I)
205: *
206: ALPHA = -T( I, 1 )
207: CALL DGEMV( 'T', M-I+1, I-1, ALPHA, A( I, 1 ), LDA,
208: $ A( I, I ), 1, ZERO, T( 1, I ), 1 )
209: A( I, I ) = AII
210: *
211: * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
212: *
213: CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
214: *
215: * T(I,I) = tau(I)
216: *
217: T( I, I ) = T( I, 1 )
218: T( I, 1) = ZERO
219: END DO
220:
221: *
222: * End of DGEQRT2
223: *
224: END
CVSweb interface <joel.bertrand@systella.fr>