File:
[local] /
rpl /
lapack /
lapack /
dgeqrt2.f
Revision
1.5:
download - view:
text,
annotated -
select for diffs -
revision graph
Mon Jan 27 09:28:16 2014 UTC (10 years, 4 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_24,
rpl-4_1_23,
rpl-4_1_22,
rpl-4_1_21,
rpl-4_1_20,
rpl-4_1_19,
rpl-4_1_18,
rpl-4_1_17,
HEAD
Cohérence.
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: *> \date September 2012
101: *
102: *> \ingroup doubleGEcomputational
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**T
123: *>
124: *> where V**T is the transpose of V.
125: *> \endverbatim
126: *>
127: * =====================================================================
128: SUBROUTINE DGEQRT2( M, N, A, LDA, T, LDT, INFO )
129: *
130: * -- LAPACK computational routine (version 3.4.2) --
131: * -- LAPACK is a software package provided by Univ. of Tennessee, --
132: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133: * September 2012
134: *
135: * .. Scalar Arguments ..
136: INTEGER INFO, LDA, LDT, M, N
137: * ..
138: * .. Array Arguments ..
139: DOUBLE PRECISION A( LDA, * ), T( LDT, * )
140: * ..
141: *
142: * =====================================================================
143: *
144: * .. Parameters ..
145: DOUBLE PRECISION ONE, ZERO
146: PARAMETER( ONE = 1.0D+00, ZERO = 0.0D+00 )
147: * ..
148: * .. Local Scalars ..
149: INTEGER I, K
150: DOUBLE PRECISION AII, ALPHA
151: * ..
152: * .. External Subroutines ..
153: EXTERNAL DLARFG, DGEMV, DGER, DTRMV, 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( 'DGEQRT2', -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 DLARFG( 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 DGEMV( 'T',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 = -(T( I, 1 ))
197: CALL DGER( 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)**T * A(I:M,I)
208: *
209: ALPHA = -T( I, 1 )
210: CALL DGEMV( 'T', 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 DTRMV( '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 DGEQRT2
226: *
227: END
CVSweb interface <joel.bertrand@systella.fr>