1 *> \brief \b CTPQRT2 computes a QR factorization of a real or complex "triangular-pentagonal" matrix, which is composed of a triangular block and a pentagonal block, using the compact WY representation for Q.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CTPQRT2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctpqrt2.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctpqrt2.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctpqrt2.f">
21 * SUBROUTINE CTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
23 * .. Scalar Arguments ..
24 * INTEGER INFO, LDA, LDB, LDT, N, M, L
26 * .. Array Arguments ..
27 * COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * )
36 *> CTPQRT2 computes a QR factorization of a complex "triangular-pentagonal"
37 *> matrix C, which is composed of a triangular block A and pentagonal block B,
38 *> using the compact WY representation for Q.
47 *> The total number of rows of the matrix B.
54 *> The number of columns of the matrix B, and the order of
55 *> the triangular matrix A.
62 *> The number of rows of the upper trapezoidal part of B.
63 *> MIN(M,N) >= L >= 0. See Further Details.
68 *> A is COMPLEX array, dimension (LDA,N)
69 *> On entry, the upper triangular N-by-N matrix A.
70 *> On exit, the elements on and above the diagonal of the array
71 *> contain the upper triangular matrix R.
77 *> The leading dimension of the array A. LDA >= max(1,N).
82 *> B is COMPLEX array, dimension (LDB,N)
83 *> On entry, the pentagonal M-by-N matrix B. The first M-L rows
84 *> are rectangular, and the last L rows are upper trapezoidal.
85 *> On exit, B contains the pentagonal matrix V. See Further Details.
91 *> The leading dimension of the array B. LDB >= max(1,M).
96 *> T is COMPLEX array, dimension (LDT,N)
97 *> The N-by-N upper triangular factor T of the block reflector.
98 *> See Further Details.
104 *> The leading dimension of the array T. LDT >= max(1,N)
110 *> = 0: successful exit
111 *> < 0: if INFO = -i, the i-th argument had an illegal value
117 *> \author Univ. of Tennessee
118 *> \author Univ. of California Berkeley
119 *> \author Univ. of Colorado Denver
122 *> \date September 2012
124 *> \ingroup complexOTHERcomputational
126 *> \par Further Details:
127 * =====================
131 *> The input matrix C is a (N+M)-by-N matrix
136 *> where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
137 *> matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
138 *> upper trapezoidal matrix B2:
140 *> B = [ B1 ] <- (M-L)-by-N rectangular
141 *> [ B2 ] <- L-by-N upper trapezoidal.
143 *> The upper trapezoidal matrix B2 consists of the first L rows of a
144 *> N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N). If L=0,
145 *> B is rectangular M-by-N; if M=L=N, B is upper triangular.
147 *> The matrix W stores the elementary reflectors H(i) in the i-th column
148 *> below the diagonal (of A) in the (N+M)-by-N input matrix C
150 *> C = [ A ] <- upper triangular N-by-N
151 *> [ B ] <- M-by-N pentagonal
153 *> so that W can be represented as
155 *> W = [ I ] <- identity, N-by-N
156 *> [ V ] <- M-by-N, same form as B.
158 *> Thus, all of information needed for W is contained on exit in B, which
159 *> we call V above. Note that V has the same form as B; that is,
161 *> V = [ V1 ] <- (M-L)-by-N rectangular
162 *> [ V2 ] <- L-by-N upper trapezoidal.
164 *> The columns of V represent the vectors which define the H(i)'s.
165 *> The (M+N)-by-(M+N) block reflector H is then given by
167 *> H = I - W * T * W**H
169 *> where W**H is the conjugate transpose of W and T is the upper triangular
170 *> factor of the block reflector.
173 * =====================================================================
174 SUBROUTINE CTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
176 * -- LAPACK computational routine (version 3.4.2) --
177 * -- LAPACK is a software package provided by Univ. of Tennessee, --
178 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
181 * .. Scalar Arguments ..
182 INTEGER INFO, LDA, LDB, LDT, N, M, L
184 * .. Array Arguments ..
185 COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * )
188 * =====================================================================
192 PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
194 * .. Local Scalars ..
195 INTEGER I, J, P, MP, NP
198 * .. External Subroutines ..
199 EXTERNAL CLARFG, CGEMV, CGERC, CTRMV, XERBLA
201 * .. Intrinsic Functions ..
204 * .. Executable Statements ..
206 * Test the input arguments
211 ELSE IF( N.LT.0 ) THEN
213 ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
215 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
217 ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
219 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
223 CALL XERBLA( 'CTPQRT2', -INFO )
227 * Quick return if possible
229 IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
233 * Generate elementary reflector H(I) to annihilate B(:,I)
236 CALL CLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
239 * W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)]
242 T( J, N ) = CONJG(A( I, I+J ))
244 CALL CGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB,
245 $ B( 1, I ), 1, ONE, T( 1, N ), 1 )
247 * C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H
249 ALPHA = -CONJG(T( I, 1 ))
251 A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N ))
253 CALL CGERC( P, N-I, ALPHA, B( 1, I ), 1,
254 $ T( 1, N ), 1, B( 1, I+1 ), LDB )
260 * T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I))
271 * Triangular part of B2
274 T( J, I ) = ALPHA*B( M-L+J, I )
276 CALL CTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB,
279 * Rectangular part of B2
281 CALL CGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB,
282 $ B( MP, I ), 1, ZERO, T( NP, I ), 1 )
286 CALL CGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1,
287 $ ONE, T( 1, I ), 1 )
289 * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
291 CALL CTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
295 T( I, I ) = T( I, 1 )