1 *> \brief \b DGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DGEQRT3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrt3.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrt3.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrt3.f">
21 * RECURSIVE SUBROUTINE DGEQRT3( M, N, A, LDA, T, LDT, INFO )
23 * .. Scalar Arguments ..
24 * INTEGER INFO, LDA, M, N, LDT
26 * .. Array Arguments ..
27 * DOUBLE PRECISION A( LDA, * ), T( LDT, * )
36 *> DGEQRT3 recursively computes a QR factorization of a real M-by-N
37 *> matrix A, using the compact WY representation of Q.
39 *> Based on the algorithm of Elmroth and Gustavson,
40 *> IBM J. Res. Develop. Vol 44 No. 4 July 2000.
49 *> The number of rows of the matrix A. M >= N.
55 *> The number of columns of the matrix A. N >= 0.
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
70 *> The leading dimension of the array A. LDA >= max(1,M).
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.
85 *> The leading dimension of the array T. LDT >= max(1,N).
91 *> = 0: successful exit
92 *> < 0: if INFO = -i, the i-th argument had an illegal value
98 *> \author Univ. of Tennessee
99 *> \author Univ. of California Berkeley
100 *> \author Univ. of Colorado Denver
105 *> \ingroup doubleGEcomputational
107 *> \par Further Details:
108 * =====================
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
121 *> where the vi's represent the vectors which define H(i), which are returned
122 *> in the matrix A. The 1's along the diagonal of V are not stored in A. The
123 *> block reflector H is then given by
125 *> H = I - V * T * V**T
127 *> where V**T is the transpose of V.
129 *> For details of the algorithm, see Elmroth and Gustavson (cited above).
132 * =====================================================================
133 RECURSIVE SUBROUTINE DGEQRT3( M, N, A, LDA, T, LDT, INFO )
135 * -- LAPACK computational routine (version 3.6.1) --
136 * -- LAPACK is a software package provided by Univ. of Tennessee, --
137 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140 * .. Scalar Arguments ..
141 INTEGER INFO, LDA, M, N, LDT
143 * .. Array Arguments ..
144 DOUBLE PRECISION A( LDA, * ), T( LDT, * )
147 * =====================================================================
151 PARAMETER ( ONE = 1.0D+00 )
153 * .. Local Scalars ..
154 INTEGER I, I1, J, J1, N1, N2, IINFO
156 * .. External Subroutines ..
157 EXTERNAL DLARFG, DTRMM, DGEMM, XERBLA
159 * .. Executable Statements ..
164 ELSE IF( M .LT. N ) THEN
166 ELSE IF( LDA .LT. MAX( 1, M ) ) THEN
168 ELSE IF( LDT .LT. MAX( 1, N ) ) THEN
172 CALL XERBLA( 'DGEQRT3', -INFO )
178 * Compute Householder transform when N=1
180 CALL DLARFG( M, A(1,1), A( MIN( 2, M ), 1 ), 1, T(1,1) )
184 * Otherwise, split A into blocks...
191 * Compute A(1:M,1:N1) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1^H
193 CALL DGEQRT3( M, N1, A, LDA, T, LDT, IINFO )
195 * Compute A(1:M,J1:N) = Q1^H A(1:M,J1:N) [workspace: T(1:N1,J1:N)]
199 T( I, J+N1 ) = A( I, J+N1 )
202 CALL DTRMM( 'L', 'L', 'T', 'U', N1, N2, ONE,
203 & A, LDA, T( 1, J1 ), LDT )
205 CALL DGEMM( 'T', 'N', N1, N2, M-N1, ONE, A( J1, 1 ), LDA,
206 & A( J1, J1 ), LDA, ONE, T( 1, J1 ), LDT)
208 CALL DTRMM( 'L', 'U', 'T', 'N', N1, N2, ONE,
209 & T, LDT, T( 1, J1 ), LDT )
211 CALL DGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( J1, 1 ), LDA,
212 & T( 1, J1 ), LDT, ONE, A( J1, J1 ), LDA )
214 CALL DTRMM( 'L', 'L', 'N', 'U', N1, N2, ONE,
215 & A, LDA, T( 1, J1 ), LDT )
219 A( I, J+N1 ) = A( I, J+N1 ) - T( I, J+N1 )
223 * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2^H
225 CALL DGEQRT3( M-N1, N2, A( J1, J1 ), LDA,
226 & T( J1, J1 ), LDT, IINFO )
228 * Compute T3 = T(1:N1,J1:N) = -T1 Y1^H Y2 T2
232 T( I, J+N1 ) = (A( J+N1, I ))
236 CALL DTRMM( 'R', 'L', 'N', 'U', N1, N2, ONE,
237 & A( J1, J1 ), LDA, T( 1, J1 ), LDT )
239 CALL DGEMM( 'T', 'N', N1, N2, M-N, ONE, A( I1, 1 ), LDA,
240 & A( I1, J1 ), LDA, ONE, T( 1, J1 ), LDT )
242 CALL DTRMM( 'L', 'U', 'N', 'N', N1, N2, -ONE, T, LDT,
245 CALL DTRMM( 'R', 'U', 'N', 'N', N1, N2, ONE,
246 & T( J1, J1 ), LDT, T( 1, J1 ), LDT )
248 * Y = (Y1,Y2); R = [ R1 A(1:N1,J1:N) ]; T = [T1 T3]