Annotation of rpl/lapack/lapack/dtpqrt2.f, revision 1.5
1.3 bertrand 1: *> \brief \b DTPQRT2 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.
1.1 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTPQRT2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtpqrt2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtpqrt2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtpqrt2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, LDB, LDT, N, M, L
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DTPQRT2 computes a QR factorization of a real "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.
39: *> \endverbatim
40: *
41: * Arguments:
42: * ==========
43: *
44: *> \param[in] M
45: *> \verbatim
46: *> M is INTEGER
47: *> The total number of rows of the matrix B.
48: *> M >= 0.
49: *> \endverbatim
50: *>
51: *> \param[in] N
52: *> \verbatim
53: *> N is INTEGER
54: *> The number of columns of the matrix B, and the order of
55: *> the triangular matrix A.
56: *> N >= 0.
57: *> \endverbatim
58: *>
59: *> \param[in] L
60: *> \verbatim
61: *> L is INTEGER
62: *> The number of rows of the upper trapezoidal part of B.
63: *> MIN(M,N) >= L >= 0. See Further Details.
64: *> \endverbatim
65: *>
66: *> \param[in,out] A
67: *> \verbatim
68: *> A is DOUBLE PRECISION 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.
72: *> \endverbatim
73: *>
74: *> \param[in] LDA
75: *> \verbatim
76: *> LDA is INTEGER
77: *> The leading dimension of the array A. LDA >= max(1,N).
78: *> \endverbatim
79: *>
80: *> \param[in,out] B
81: *> \verbatim
82: *> B is DOUBLE PRECISION 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.
86: *> \endverbatim
87: *>
88: *> \param[in] LDB
89: *> \verbatim
90: *> LDB is INTEGER
91: *> The leading dimension of the array B. LDB >= max(1,M).
92: *> \endverbatim
93: *>
94: *> \param[out] T
95: *> \verbatim
96: *> T is DOUBLE PRECISION array, dimension (LDT,N)
97: *> The N-by-N upper triangular factor T of the block reflector.
98: *> See Further Details.
99: *> \endverbatim
100: *>
101: *> \param[in] LDT
102: *> \verbatim
103: *> LDT is INTEGER
104: *> The leading dimension of the array T. LDT >= max(1,N)
105: *> \endverbatim
106: *>
107: *> \param[out] INFO
108: *> \verbatim
109: *> INFO is INTEGER
110: *> = 0: successful exit
111: *> < 0: if INFO = -i, the i-th argument had an illegal value
112: *> \endverbatim
113: *
114: * Authors:
115: * ========
116: *
117: *> \author Univ. of Tennessee
118: *> \author Univ. of California Berkeley
119: *> \author Univ. of Colorado Denver
120: *> \author NAG Ltd.
121: *
1.3 bertrand 122: *> \date September 2012
1.1 bertrand 123: *
124: *> \ingroup doubleOTHERcomputational
125: *
126: *> \par Further Details:
127: * =====================
128: *>
129: *> \verbatim
130: *>
131: *> The input matrix C is a (N+M)-by-N matrix
132: *>
133: *> C = [ A ]
134: *> [ B ]
135: *>
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:
139: *>
140: *> B = [ B1 ] <- (M-L)-by-N rectangular
141: *> [ B2 ] <- L-by-N upper trapezoidal.
142: *>
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.
146: *>
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
149: *>
150: *> C = [ A ] <- upper triangular N-by-N
151: *> [ B ] <- M-by-N pentagonal
152: *>
153: *> so that W can be represented as
154: *>
155: *> W = [ I ] <- identity, N-by-N
156: *> [ V ] <- M-by-N, same form as B.
157: *>
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,
160: *>
161: *> V = [ V1 ] <- (M-L)-by-N rectangular
162: *> [ V2 ] <- L-by-N upper trapezoidal.
163: *>
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
166: *>
167: *> H = I - W * T * W**T
168: *>
169: *> where W^H is the conjugate transpose of W and T is the upper triangular
170: *> factor of the block reflector.
171: *> \endverbatim
172: *>
173: * =====================================================================
174: SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
175: *
1.3 bertrand 176: * -- LAPACK computational routine (version 3.4.2) --
1.1 bertrand 177: * -- LAPACK is a software package provided by Univ. of Tennessee, --
178: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.3 bertrand 179: * September 2012
1.1 bertrand 180: *
181: * .. Scalar Arguments ..
182: INTEGER INFO, LDA, LDB, LDT, N, M, L
183: * ..
184: * .. Array Arguments ..
185: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * )
186: * ..
187: *
188: * =====================================================================
189: *
190: * .. Parameters ..
191: DOUBLE PRECISION ONE, ZERO
192: PARAMETER( ONE = 1.0, ZERO = 0.0 )
193: * ..
194: * .. Local Scalars ..
195: INTEGER I, J, P, MP, NP
196: DOUBLE PRECISION ALPHA
197: * ..
198: * .. External Subroutines ..
199: EXTERNAL DLARFG, DGEMV, DGER, DTRMV, XERBLA
200: * ..
201: * .. Intrinsic Functions ..
202: INTRINSIC MAX, MIN
203: * ..
204: * .. Executable Statements ..
205: *
206: * Test the input arguments
207: *
208: INFO = 0
209: IF( M.LT.0 ) THEN
210: INFO = -1
211: ELSE IF( N.LT.0 ) THEN
212: INFO = -2
213: ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
214: INFO = -3
215: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
216: INFO = -5
217: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
218: INFO = -7
219: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
220: INFO = -9
221: END IF
222: IF( INFO.NE.0 ) THEN
223: CALL XERBLA( 'DTPQRT2', -INFO )
224: RETURN
225: END IF
226: *
227: * Quick return if possible
228: *
229: IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
230: *
231: DO I = 1, N
232: *
233: * Generate elementary reflector H(I) to annihilate B(:,I)
234: *
235: P = M-L+MIN( L, I )
236: CALL DLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
237: IF( I.LT.N ) THEN
238: *
239: * W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)]
240: *
241: DO J = 1, N-I
242: T( J, N ) = (A( I, I+J ))
243: END DO
244: CALL DGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB,
245: $ B( 1, I ), 1, ONE, T( 1, N ), 1 )
246: *
247: * C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H
248: *
249: ALPHA = -(T( I, 1 ))
250: DO J = 1, N-I
251: A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N ))
252: END DO
253: CALL DGER( P, N-I, ALPHA, B( 1, I ), 1,
254: $ T( 1, N ), 1, B( 1, I+1 ), LDB )
255: END IF
256: END DO
257: *
258: DO I = 2, N
259: *
260: * T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I))
261: *
262: ALPHA = -T( I, 1 )
263:
264: DO J = 1, I-1
265: T( J, I ) = ZERO
266: END DO
267: P = MIN( I-1, L )
268: MP = MIN( M-L+1, M )
269: NP = MIN( P+1, N )
270: *
271: * Triangular part of B2
272: *
273: DO J = 1, P
274: T( J, I ) = ALPHA*B( M-L+J, I )
275: END DO
276: CALL DTRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB,
277: $ T( 1, I ), 1 )
278: *
279: * Rectangular part of B2
280: *
281: CALL DGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB,
282: $ B( MP, I ), 1, ZERO, T( NP, I ), 1 )
283: *
284: * B1
285: *
286: CALL DGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1,
287: $ ONE, T( 1, I ), 1 )
288: *
289: * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
290: *
291: CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
292: *
293: * T(I,I) = tau(I)
294: *
295: T( I, I ) = T( I, 1 )
296: T( I, 1 ) = ZERO
297: END DO
298:
299: *
300: * End of DTPQRT2
301: *
302: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtpqrt2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack