diff options
author | james <james@8a072113-8704-0410-8d35-dd094bca7971> | 2011-08-08 22:07:49 +0000 |
---|---|---|
committer | james <james@8a072113-8704-0410-8d35-dd094bca7971> | 2011-08-08 22:07:49 +0000 |
commit | 1cb631b23a3241a10238242418956cac81043a2f (patch) | |
tree | 8a332287b2a06c5b971e741603d08c64730667c4 /SRC/ctpqrt.f | |
parent | ccebfeaeba01e127c3e80ab8f52e4d243a62d33d (diff) | |
download | lapack-1cb631b23a3241a10238242418956cac81043a2f.tar.gz lapack-1cb631b23a3241a10238242418956cac81043a2f.tar.bz2 lapack-1cb631b23a3241a10238242418956cac81043a2f.zip |
QR factorization for triangular-pentagonal matrices
(generalization of triangle-over-square and triangle-over-triangle)
Diffstat (limited to 'SRC/ctpqrt.f')
-rw-r--r-- | SRC/ctpqrt.f | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/SRC/ctpqrt.f b/SRC/ctpqrt.f new file mode 100644 index 00000000..ada345b4 --- /dev/null +++ b/SRC/ctpqrt.f @@ -0,0 +1,185 @@ + SUBROUTINE CTPQRT( M, N, L, NB, A, LDA, B, LDB, T, LDT, WORK, + $ INFO ) + IMPLICIT NONE +* +* -- LAPACK routine (version 3.?) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* -- July 2011 -- +* +* .. Scalar Arguments .. + INTEGER INFO, LDA, LDB, LDT, N, M, L, NB +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* CTPQRT computes a blocked QR factorization of a complex +* "triangular-pentagonal" matrix C, which is composed of a +* triangular block A and pentagonal block B, using the compact +* WY representation for Q. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix B. +* M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B, and the order of the +* triangular matrix A. +* N >= 0. +* +* L (input) INTEGER +* The number of rows of the upper trapezoidal part of B. +* MIN(M,N) >= L >= 0. See Further Details. +* +* NB (input) INTEGER +* The block size to be used in the blocked QR. N >= NB >= 1. +* +* A (input/output) COMPLEX array, dimension (LDA,N) +* On entry, the upper triangular N-by-N matrix A. +* On exit, the elements on and above the diagonal of the array +* contain the upper triangular matrix R. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* B (input/output) COMPLEX array, dimension (LDB,N) +* On entry, the pentagonal M-by-N matrix B. The first M-L rows +* are rectangular, and the last L rows are upper trapezoidal. +* On exit, B contains the pentagonal matrix V. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. LDB >= max(1,M). +* +* T (output) COMPLEX array, dimension (LDT,N) +* The upper triangular block reflectors stored in compact form +* as a sequence of upper triangular blocks. See Further Details. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= NB. +* +* WORK (workspace) COMPLEX array, dimension (NB*N) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The input matrix C is a (N+M)-by-N matrix +* +* C = [ A ] +* [ B ] +* +* where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal +* matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N +* upper trapezoidal matrix B2: +* +* B = [ B1 ] <- (M-L)-by-N rectangular +* [ B2 ] <- L-by-N upper trapezoidal. +* +* The upper trapezoidal matrix B2 consists of the first L rows of a +* N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N). If L=0, +* B is rectangular M-by-N; if M=L=N, B is upper triangular. +* +* The matrix W stores the elementary reflectors H(i) in the i-th column +* below the diagonal (of A) in the (N+M)-by-N input matrix C +* +* C = [ A ] <- upper triangular N-by-N +* [ B ] <- M-by-N pentagonal +* +* so that W can be represented as +* +* W = [ I ] <- identity, N-by-N +* [ V ] <- M-by-N, same form as B. +* +* Thus, all of information needed for W is contained on exit in B, which +* we call V above. Note that V has the same form as B; that is, +* +* V = [ V1 ] <- (M-L)-by-N rectangular +* [ V2 ] <- L-by-N upper trapezoidal. +* +* The columns of V represent the vectors which define the H(i)'s. +* +* The number of blocks is B = ceiling(N/NB), where each +* block is of order NB except for the last block, which is of order +* IB = N - (B-1)*NB. For each of the B blocks, a upper triangular block +* reflector factor is computed: T1, T2, ..., TB. The NB-by-NB (and IB-by-IB +* for the last block) T's are stored in the NB-by-N matrix T as +* +* T = [T1 T2 ... TB]. +* +* ===================================================================== +* +* .. +* .. Local Scalars .. + INTEGER I, IB, LB, MB, IINFO +* .. +* .. External Subroutines .. + EXTERNAL CTPQRT2, CTPRFB, XERBLA +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( M.LT.0 ) THEN + INFO = -2 + ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN + INFO = -3 + ELSE IF( NB.LT.1 .OR. NB.GT.N ) THEN + INFO = -4 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -6 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -8 + ELSE IF( LDT.LT.NB ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CTPQRT', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.EQ.0 .OR. N.EQ.0 ) RETURN +* + DO I = 1, N, NB +* +* Compute the QR factorization of the current block +* + IB = MIN( N-I+1, NB ) + MB = MIN( M-L+I+IB-1, M ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-M+L-I+1 + END IF +* + CALL CTPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB, + $ T(1, I ), LDT, IINFO ) +* +* Update by applying H**H to B(:,I+IB:N) from the left +* + IF( I+IB.LE.N ) THEN + CALL CTPRFB( 'L', 'C', 'F', 'C', MB, N-I-IB+1, IB, LB, + $ B( 1, I ), LDB, T( 1, I ), LDT, + $ A( I, I+IB ), LDA, B( 1, I+IB ), LDB, + $ WORK, IB ) + END IF + END DO + RETURN +* +* End of CTPQRT +* + END |