Annotation of rpl/lapack/lapack/dgeqrt3.f, revision 1.11
1.3 bertrand 1: *> \brief \b DGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q.
1.1 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.8 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: *> \htmlonly
1.8 bertrand 9: *> Download DGEQRT3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrt3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrt3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrt3.f">
1.1 bertrand 15: *> [TXT]</a>
1.8 bertrand 16: *> \endhtmlonly
1.1 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * RECURSIVE SUBROUTINE DGEQRT3( M, N, A, LDA, T, LDT, INFO )
1.8 bertrand 22: *
1.1 bertrand 23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, M, N, LDT
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), T( LDT, * )
28: * ..
1.8 bertrand 29: *
1.1 bertrand 30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
1.8 bertrand 36: *> DGEQRT3 recursively computes a QR factorization of a real M-by-N
37: *> matrix A, using the compact WY representation of Q.
1.1 bertrand 38: *>
1.8 bertrand 39: *> Based on the algorithm of Elmroth and Gustavson,
1.1 bertrand 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 DOUBLE PRECISION array, dimension (LDA,N)
61: *> On entry, the real M-by-N matrix A. On exit, the elements on and
62: *> above the diagonal contain the N-by-N upper triangular matrix R; the
63: *> elements below the diagonal are the columns 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 DOUBLE PRECISION 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: *
1.8 bertrand 98: *> \author Univ. of Tennessee
99: *> \author Univ. of California Berkeley
100: *> \author Univ. of Colorado Denver
101: *> \author NAG Ltd.
1.1 bertrand 102: *
103: *> \ingroup doubleGEcomputational
104: *
105: *> \par Further Details:
106: * =====================
107: *>
108: *> \verbatim
109: *>
110: *> The matrix V stores the elementary reflectors H(i) in the i-th column
111: *> below the diagonal. For example, if M=5 and N=3, the matrix V is
112: *>
113: *> V = ( 1 )
114: *> ( v1 1 )
115: *> ( v1 v2 1 )
116: *> ( v1 v2 v3 )
117: *> ( v1 v2 v3 )
118: *>
119: *> where the vi's represent the vectors which define H(i), which are returned
120: *> in the matrix A. The 1's along the diagonal of V are not stored in A. The
121: *> block reflector H is then given by
122: *>
123: *> H = I - V * T * V**T
124: *>
125: *> where V**T is the transpose of V.
126: *>
127: *> For details of the algorithm, see Elmroth and Gustavson (cited above).
128: *> \endverbatim
129: *>
130: * =====================================================================
131: RECURSIVE SUBROUTINE DGEQRT3( M, N, A, LDA, T, LDT, INFO )
132: *
1.11 ! bertrand 133: * -- LAPACK computational routine --
1.1 bertrand 134: * -- LAPACK is a software package provided by Univ. of Tennessee, --
135: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136: *
137: * .. Scalar Arguments ..
138: INTEGER INFO, LDA, M, N, LDT
139: * ..
140: * .. Array Arguments ..
141: DOUBLE PRECISION A( LDA, * ), T( LDT, * )
142: * ..
143: *
144: * =====================================================================
145: *
146: * .. Parameters ..
147: DOUBLE PRECISION ONE
148: PARAMETER ( ONE = 1.0D+00 )
149: * ..
150: * .. Local Scalars ..
151: INTEGER I, I1, J, J1, N1, N2, IINFO
152: * ..
153: * .. External Subroutines ..
154: EXTERNAL DLARFG, DTRMM, DGEMM, XERBLA
155: * ..
156: * .. Executable Statements ..
157: *
158: INFO = 0
159: IF( N .LT. 0 ) THEN
160: INFO = -2
161: ELSE IF( M .LT. N ) THEN
162: INFO = -1
163: ELSE IF( LDA .LT. MAX( 1, M ) ) THEN
164: INFO = -4
165: ELSE IF( LDT .LT. MAX( 1, N ) ) THEN
166: INFO = -6
167: END IF
168: IF( INFO.NE.0 ) THEN
169: CALL XERBLA( 'DGEQRT3', -INFO )
170: RETURN
171: END IF
172: *
173: IF( N.EQ.1 ) THEN
174: *
175: * Compute Householder transform when N=1
176: *
1.6 bertrand 177: CALL DLARFG( M, A(1,1), A( MIN( 2, M ), 1 ), 1, T(1,1) )
1.8 bertrand 178: *
1.1 bertrand 179: ELSE
180: *
181: * Otherwise, split A into blocks...
182: *
183: N1 = N/2
184: N2 = N-N1
185: J1 = MIN( N1+1, N )
186: I1 = MIN( N+1, M )
187: *
188: * Compute A(1:M,1:N1) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1^H
189: *
190: CALL DGEQRT3( M, N1, A, LDA, T, LDT, IINFO )
191: *
192: * Compute A(1:M,J1:N) = Q1^H A(1:M,J1:N) [workspace: T(1:N1,J1:N)]
193: *
194: DO J=1,N2
195: DO I=1,N1
196: T( I, J+N1 ) = A( I, J+N1 )
197: END DO
198: END DO
1.8 bertrand 199: CALL DTRMM( 'L', 'L', 'T', 'U', N1, N2, ONE,
1.1 bertrand 200: & A, LDA, T( 1, J1 ), LDT )
201: *
202: CALL DGEMM( 'T', 'N', N1, N2, M-N1, ONE, A( J1, 1 ), LDA,
203: & A( J1, J1 ), LDA, ONE, T( 1, J1 ), LDT)
204: *
205: CALL DTRMM( 'L', 'U', 'T', 'N', N1, N2, ONE,
206: & T, LDT, T( 1, J1 ), LDT )
207: *
1.8 bertrand 208: CALL DGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( J1, 1 ), LDA,
1.1 bertrand 209: & T( 1, J1 ), LDT, ONE, A( J1, J1 ), LDA )
210: *
211: CALL DTRMM( 'L', 'L', 'N', 'U', N1, N2, ONE,
212: & A, LDA, T( 1, J1 ), LDT )
213: *
214: DO J=1,N2
215: DO I=1,N1
216: A( I, J+N1 ) = A( I, J+N1 ) - T( I, J+N1 )
217: END DO
218: END DO
219: *
220: * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2^H
221: *
1.8 bertrand 222: CALL DGEQRT3( M-N1, N2, A( J1, J1 ), LDA,
1.1 bertrand 223: & T( J1, J1 ), LDT, IINFO )
224: *
225: * Compute T3 = T(1:N1,J1:N) = -T1 Y1^H Y2 T2
226: *
227: DO I=1,N1
228: DO J=1,N2
229: T( I, J+N1 ) = (A( J+N1, I ))
230: END DO
231: END DO
232: *
233: CALL DTRMM( 'R', 'L', 'N', 'U', N1, N2, ONE,
234: & A( J1, J1 ), LDA, T( 1, J1 ), LDT )
235: *
1.8 bertrand 236: CALL DGEMM( 'T', 'N', N1, N2, M-N, ONE, A( I1, 1 ), LDA,
1.1 bertrand 237: & A( I1, J1 ), LDA, ONE, T( 1, J1 ), LDT )
238: *
1.8 bertrand 239: CALL DTRMM( 'L', 'U', 'N', 'N', N1, N2, -ONE, T, LDT,
1.1 bertrand 240: & T( 1, J1 ), LDT )
241: *
1.8 bertrand 242: CALL DTRMM( 'R', 'U', 'N', 'N', N1, N2, ONE,
1.1 bertrand 243: & T( J1, J1 ), LDT, T( 1, J1 ), LDT )
244: *
245: * Y = (Y1,Y2); R = [ R1 A(1:N1,J1:N) ]; T = [T1 T3]
246: * [ 0 R2 ] [ 0 T2]
247: *
248: END IF
249: *
250: RETURN
251: *
252: * End of DGEQRT3
253: *
254: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgeqrt3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack