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 | |
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')
-rw-r--r-- | SRC/ctpmqrt.f | 264 | ||||
-rw-r--r-- | SRC/ctpqrt.f | 185 | ||||
-rw-r--r-- | SRC/ctpqrt2.f | 224 | ||||
-rw-r--r-- | SRC/ctprfb.f | 708 | ||||
-rw-r--r-- | SRC/dtpmqrt.f | 264 | ||||
-rw-r--r-- | SRC/dtpqrt.f | 185 | ||||
-rw-r--r-- | SRC/dtpqrt2.f | 224 | ||||
-rw-r--r-- | SRC/dtprfb.f | 705 | ||||
-rw-r--r-- | SRC/stpmqrt.f | 264 | ||||
-rw-r--r-- | SRC/stpqrt.f | 185 | ||||
-rw-r--r-- | SRC/stpqrt2.f | 224 | ||||
-rw-r--r-- | SRC/stprfb.f | 705 | ||||
-rw-r--r-- | SRC/ztpmqrt.f | 264 | ||||
-rw-r--r-- | SRC/ztpqrt.f | 185 | ||||
-rw-r--r-- | SRC/ztpqrt2.f | 224 | ||||
-rw-r--r-- | SRC/ztprfb.f | 708 |
16 files changed, 5518 insertions, 0 deletions
diff --git a/SRC/ctpmqrt.f b/SRC/ctpmqrt.f new file mode 100644 index 00000000..1bcd7b7f --- /dev/null +++ b/SRC/ctpmqrt.f @@ -0,0 +1,264 @@ + SUBROUTINE CTPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, + $ A, LDA, B, LDB, 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 .. + CHARACTER SIDE, TRANS + INTEGER INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT +* .. +* .. Array Arguments .. + COMPLEX V( LDV, * ), A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ WORK( * ) +* .. +* +* Purpose +* ======= +* +* CTPMQRT applies a complex orthogonal matrix Q obtained from a +* "triangular-pentagonal" complex block reflector H to a general +* complex matrix C, which consists of two blocks A and B. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply Q or Q**H from the Left; +* = 'R': apply Q or Q**H from the Right. +* +* TRANS (input) CHARACTER*1 +* = 'N': No transpose, apply Q; +* = 'C': Transpose, apply Q**H. +* +* M (input) INTEGER +* The number of rows of the matrix B. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. N >= 0. +* +* K (input) INTEGER +* The number of elementary reflectors whose product defines +* the matrix Q. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* NB (input) INTEGER +* The block size used for the storage of T. K >= NB >= 1. +* This must be the same value of NB used to generate T +* in CTPQRT. +* +* V (input) COMPLEX array, dimension (LDA,K) +* The i-th column must contain the vector which defines the +* elementary reflector H(i), for i = 1,2,...,k, as returned by +* CTPQRT in B. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDA >= max(1,M); +* if SIDE = 'R', LDA >= max(1,N). +* +* T (input) COMPLEX array, dimension (LDT,K) +* The upper triangular factors of the block reflectors +* as returned by CTPQRT, stored as a NB-by-K matrix. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= NB. +* +* A (input/output) COMPLEX array, dimension +* (LDA,N) if SIDE = 'L' or +* (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* Q*C or Q**H*C or C*Q or C*Q**H. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) COMPLEX array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* Q*C or Q**H*C or C*Q or C*Q**H. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace/output) COMPLEX array. The dimension of WORK is +* N*NB if SIDE = 'L', or M*NB if SIDE = 'R'. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The columns of the pentagonal matrix V contain the elementary reflectors +* H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a +* trapezoidal block V2: +* +* V = [V1] +* [V2]. +* +* The size of the trapezoidal block V2 is determined by the parameter L, +* where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L +* rows of a K-by-K upper triangular matrix. If L=K, V2 is upper triangular; +* if L=0, there is no trapezoidal block, hence V = V1 is rectangular. +* +* If SIDE = 'L': C = [A] where A is K-by-N, B is M-by-N and V is M-by-K. +* [B] +* +* If SIDE = 'R': C = [A B] where A is M-by-K, B is M-by-N and V is N-by-K. +* +* The complex orthogonal matrix Q is formed from V and T. +* +* If TRANS='N' and SIDE='L', C is on exit replaced with Q * C. +* +* If TRANS='C' and SIDE='L', C is on exit replaced with Q**H * C. +* +* If TRANS='N' and SIDE='R', C is on exit replaced with C * Q. +* +* If TRANS='C' and SIDE='R', C is on exit replaced with C * Q**H. +* +* ===================================================================== +* +* .. +* .. Local Scalars .. + LOGICAL LEFT, RIGHT, TRAN, NOTRAN + INTEGER I, IB, MB, LB, KF, Q +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, CLARFB +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* .. Test the input arguments .. +* + INFO = 0 + LEFT = LSAME( SIDE, 'L' ) + RIGHT = LSAME( SIDE, 'R' ) + TRAN = LSAME( TRANS, 'C' ) + NOTRAN = LSAME( TRANS, 'N' ) +* + IF( LEFT ) THEN + Q = M + ELSE IF ( RIGHT ) THEN + Q = N + END IF + IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN + INFO = -1 + ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( K.LT.0 ) THEN + INFO = -5 + ELSE IF( L.LT.0 .OR. L.GT.K ) THEN + INFO = -6 + ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN + INFO = -7 + ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN + INFO = -8 + ELSE IF( LDT.LT.NB ) THEN + INFO = -10 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -12 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -12 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CTPMQRT', -INFO ) + RETURN + END IF +* +* .. Quick return if possible .. +* + IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN +* + IF( LEFT .AND. TRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + 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 CTPRFB( 'L', 'C', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. NOTRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL CTPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + ELSE IF( LEFT .AND. NOTRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + 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 CTPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. TRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL CTPRFB( 'R', 'C', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + END IF +* + RETURN +* +* End of CTPMQRT +* + END 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 diff --git a/SRC/ctpqrt2.f b/SRC/ctpqrt2.f new file mode 100644 index 00000000..33ea7b92 --- /dev/null +++ b/SRC/ctpqrt2.f @@ -0,0 +1,224 @@ + SUBROUTINE CTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO ) + IMPLICIT NONE +* +* -- LAPACK routine (version 3.x) -- +* -- 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 +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ) +* .. +* +* Purpose +* ======= +* +* CTPQRT2 computes a 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 total 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. +* +* 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 N-by-N upper triangular factor T of the block reflector. +* See Further Details. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= max(1,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 (M+N)-by-(M+N) block reflector H is then given by +* +* H = I - W * T * W**H +* +* where W**H is the conjugate transpose of W and T is the upper triangular +* factor of the block reflector. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE, ZERO + PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) ) +* .. +* .. Local Scalars .. + INTEGER I, J, P, MP, NP + COMPLEX ALPHA +* .. +* .. External Subroutines .. + EXTERNAL CLARFG, CGEMV, CGERC, CTRMV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -7 + ELSE IF( LDT.LT.MAX( 1, N ) ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CTPQRT2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. M.EQ.0 ) RETURN +* + DO I = 1, N +* +* Generate elementary reflector H(I) to annihilate B(:,I) +* + P = M-L+MIN( L, I ) + CALL CLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) ) + IF( I.LT.N ) THEN +* +* W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)] +* + DO J = 1, N-I + T( J, N ) = CONJG(A( I, I+J )) + END DO + CALL CGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB, + $ B( 1, I ), 1, ONE, T( 1, N ), 1 ) +* +* C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H +* + ALPHA = -CONJG(T( I, 1 )) + DO J = 1, N-I + A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N )) + END DO + CALL CGERC( P, N-I, ALPHA, B( 1, I ), 1, + $ T( 1, N ), 1, B( 1, I+1 ), LDB ) + END IF + END DO +* + DO I = 2, N +* +* T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I)) +* + ALPHA = -T( I, 1 ) + + DO J = 1, I-1 + T( J, I ) = ZERO + END DO + P = MIN( I-1, L ) + MP = MIN( M-L+1, M ) + NP = MIN( P+1, N ) +* +* Triangular part of B2 +* + DO J = 1, P + T( J, I ) = ALPHA*B( M-L+J, I ) + END DO + CALL CTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB, + $ T( 1, I ), 1 ) +* +* Rectangular part of B2 +* + CALL CGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB, + $ B( MP, I ), 1, ZERO, T( NP, I ), 1 ) +* +* B1 +* + CALL CGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, + $ ONE, T( 1, I ), 1 ) +* +* T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I) +* + CALL CTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 ) +* +* T(I,I) = tau(I) +* + T( I, I ) = T( I, 1 ) + T( I, 1 ) = ZERO + END DO + +* +* End of CTPQRT2 +* + END diff --git a/SRC/ctprfb.f b/SRC/ctprfb.f new file mode 100644 index 00000000..9c7980c0 --- /dev/null +++ b/SRC/ctprfb.f @@ -0,0 +1,708 @@ + SUBROUTINE CTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, + $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine (version 3.x) -- +* -- 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 .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ V( LDV, * ), WORK( LDWORK, * ) +* .. +* +* Purpose +* ======= +* +* CTPRFB applies a complex "triangular-pentagonal" block reflector H or its +* conjugate transpose H**H to a complex matrix C, which is composed of two +* blocks A and B, either from the left or right. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply H or H**H from the Left +* = 'R': apply H or H**H from the Right +* +* TRANS (input) CHARACTER*1 +* = 'N': apply H (No transpose) +* = 'C': apply H**H (Conjugate transpose) +* +* DIRECT (input) CHARACTER*1 +* Indicates how H is formed from a product of elementary +* reflectors +* = 'F': H = H(1) H(2) . . . H(k) (Forward) +* = 'B': H = H(k) . . . H(2) H(1) (Backward) +* +* STOREV (input) CHARACTER*1 +* Indicates how the vectors which define the elementary +* reflectors are stored: +* = 'C': Columns +* = 'R': Rows +* +* M (input) INTEGER +* The number of rows of the matrix B. +* M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. +* N >= 0. +* +* K (input) INTEGER +* The order of the matrix T, i.e. the number of elementary +* reflectors whose product defines the block reflector. +* K >= 0. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* V (input) COMPLEX array, dimension +* (LDV,K) if STOREV = 'C' +* (LDV,M) if STOREV = 'R' and SIDE = 'L' +* (LDV,N) if STOREV = 'R' and SIDE = 'R' +* The pentagonal matrix V, which contains the elementary reflectors +* H(1), H(2), ..., H(K). See Further Details. +* +* LDV (input) INTEGER +* The leading dimension of the array V. +* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +* if STOREV = 'R', LDV >= K. +* +* T (input) COMPLEX array, dimension (LDT,K) +* The triangular K-by-K matrix T in the representation of the +* block reflector. +* +* LDT (input) INTEGER +* The leading dimension of the array T. +* LDT >= K. +* +* A (input/output) COMPLEX array, dimension +* (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* H*C or H**H*C or C*H or C*H**H. See Futher Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) COMPLEX array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* H*C or H**H*C or C*H or C*H**H. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace) COMPLEX array, dimension +* (LDWORK,N) if SIDE = 'L', +* (LDWORK,K) if SIDE = 'R'. +* +* LDWORK (input) INTEGER +* The leading dimension of the array WORK. +* If SIDE = 'L', LDWORK >= K; +* if SIDE = 'R', LDWORK >= M. +* +* Further Details +* =============== +* +* The matrix C is a composite matrix formed from blocks A and B. +* The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, +* and if SIDE = 'L', A is of size K-by-N. +* +* If SIDE = 'R' and DIRECT = 'F', C = [A B]. +* +* If SIDE = 'L' and DIRECT = 'F', C = [A] +* [B]. +* +* If SIDE = 'R' and DIRECT = 'B', C = [B A]. +* +* If SIDE = 'L' and DIRECT = 'B', C = [B] +* [A]. +* +* The pentagonal matrix V is composed of a rectangular block V1 and a +* trapezoidal block V2. The size of the trapezoidal block is determined by +* the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; +* if L=0, there is no trapezoidal block, thus V = V1 is rectangular. +* +* If DIRECT = 'F' and STOREV = 'C': V = [V1] +* [V2] +* - V2 is upper trapezoidal (first L rows of K-by-K upper triangular) +* +* If DIRECT = 'F' and STOREV = 'R': V = [V1 V2] +* +* - V2 is lower trapezoidal (first L columns of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'C': V = [V2] +* [V1] +* - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] +* +* - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) +* +* If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. +* +* If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K. +* +* If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L. +* +* If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L. +* +* ========================================================================== +* +* .. Parameters .. + COMPLEX ONE, ZERO + PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,0.0) ) +* .. +* .. Local Scalars .. + INTEGER I, J, MP, NP, KP + LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM +* .. +* .. Intrinsic Functions .. + INTRINSIC CONJG +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN +* + IF( LSAME( STOREV, 'C' ) ) THEN + COLUMN = .TRUE. + ROW = .FALSE. + ELSE IF ( LSAME( STOREV, 'R' ) ) THEN + COLUMN = .FALSE. + ROW = .TRUE. + ELSE + COLUMN = .FALSE. + ROW = .FALSE. + END IF +* + IF( LSAME( SIDE, 'L' ) ) THEN + LEFT = .TRUE. + RIGHT = .FALSE. + ELSE IF( LSAME( SIDE, 'R' ) ) THEN + LEFT = .FALSE. + RIGHT = .TRUE. + ELSE + LEFT = .FALSE. + RIGHT = .FALSE. + END IF +* + IF( LSAME( DIRECT, 'F' ) ) THEN + FORWARD = .TRUE. + BACKWARD = .FALSE. + ELSE IF( LSAME( DIRECT, 'B' ) ) THEN + FORWARD = .FALSE. + BACKWARD = .TRUE. + ELSE + FORWARD = .FALSE. + BACKWARD = .FALSE. + END IF +* +* --------------------------------------------------------------------------- +* + IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (M-by-K) +* +* Form H C or H**H C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - T (A + V**H B) or A = A - T**H (A + V**H B) +* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL CTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + CALL CGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB, + $ ONE, WORK, LDWORK ) + CALL CGEMM( 'C', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL CGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL CTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (N-by-K) +* +* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - (A + B V) T or A = A - (A + B V) T**H +* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL CTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + CALL CGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, + $ V, LDV, ONE, WORK, LDWORK ) + CALL CGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL CGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB ) + CALL CTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (M-by-K) +* [ I ] (K-by-K) +* +* Form H C or H**H C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - T (A + V**H B) or A = A - T**H (A + V**H B) +* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO +* + CALL CTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL CGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL CGEMM( 'C', 'N', K-L, N, M, ONE, V, LDV, + $ B, LDB, ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL CGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL CTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (N-by-K) +* [ I ] (K-by-K) +* +* Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - (A + B V) T or A = A - (A + B V) T**H +* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL CTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL CGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL CGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V, LDV, ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK, + $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB ) + CALL CGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL CTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H**H C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - T (A + V B) or A = A - T**H (A + V B) +* B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL CTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDB ) + CALL CGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, + $ ONE, WORK, LDWORK ) + CALL CGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL CGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL CTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H +* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL CTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + CALL CGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV, + $ ONE, WORK, LDWORK ) + CALL CGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, + $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL CGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL CTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H**H C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - T (A + V B) or A = A - T**H (A + V B) +* B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO + CALL CTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL CGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL CGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL CGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL CTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H +* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL CTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL CGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL CGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, V, LDV, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL CTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL CGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL CGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL CTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* + END IF +* + RETURN +* +* End of CTPRFB +* + END diff --git a/SRC/dtpmqrt.f b/SRC/dtpmqrt.f new file mode 100644 index 00000000..73b257a4 --- /dev/null +++ b/SRC/dtpmqrt.f @@ -0,0 +1,264 @@ + SUBROUTINE DTPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, + $ A, LDA, B, LDB, 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 .. + CHARACTER SIDE, TRANS + INTEGER INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT +* .. +* .. Array Arguments .. + DOUBLE PRECISION V( LDV, * ), A( LDA, * ), B( LDB, * ), + $ T( LDT, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DTPMQRT applies a real orthogonal matrix Q obtained from a +* "triangular-pentagonal" real block reflector H to a general +* real matrix C, which consists of two blocks A and B. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply Q or Q**T from the Left; +* = 'R': apply Q or Q**T from the Right. +* +* TRANS (input) CHARACTER*1 +* = 'N': No transpose, apply Q; +* = 'C': Transpose, apply Q**T. +* +* M (input) INTEGER +* The number of rows of the matrix B. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. N >= 0. +* +* K (input) INTEGER +* The number of elementary reflectors whose product defines +* the matrix Q. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* NB (input) INTEGER +* The block size used for the storage of T. K >= NB >= 1. +* This must be the same value of NB used to generate T +* in CTPQRT. +* +* V (input) DOUBLE PRECISION array, dimension (LDA,K) +* The i-th column must contain the vector which defines the +* elementary reflector H(i), for i = 1,2,...,k, as returned by +* CTPQRT in B. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDA >= max(1,M); +* if SIDE = 'R', LDA >= max(1,N). +* +* T (input) DOUBLE PRECISION array, dimension (LDT,K) +* The upper triangular factors of the block reflectors +* as returned by CTPQRT, stored as a NB-by-K matrix. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= NB. +* +* A (input/output) DOUBLE PRECISION array, dimension +* (LDA,N) if SIDE = 'L' or +* (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* Q*C or Q**T*C or C*Q or C*Q**T. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) DOUBLE PRECISION array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* Q*C or Q**T*C or C*Q or C*Q**T. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace/output) DOUBLE PRECISION array. The dimension of WORK is +* N*NB if SIDE = 'L', or M*NB if SIDE = 'R'. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The columns of the pentagonal matrix V contain the elementary reflectors +* H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a +* trapezoidal block V2: +* +* V = [V1] +* [V2]. +* +* The size of the trapezoidal block V2 is determined by the parameter L, +* where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L +* rows of a K-by-K upper triangular matrix. If L=K, V2 is upper triangular; +* if L=0, there is no trapezoidal block, hence V = V1 is rectangular. +* +* If SIDE = 'L': C = [A] where A is K-by-N, B is M-by-N and V is M-by-K. +* [B] +* +* If SIDE = 'R': C = [A B] where A is M-by-K, B is M-by-N and V is N-by-K. +* +* The real orthogonal matrix Q is formed from V and T. +* +* If TRANS='N' and SIDE='L', C is on exit replaced with Q * C. +* +* If TRANS='T' and SIDE='L', C is on exit replaced with Q**T * C. +* +* If TRANS='N' and SIDE='R', C is on exit replaced with C * Q. +* +* If TRANS='T' and SIDE='R', C is on exit replaced with C * Q**T. +* +* ===================================================================== +* +* .. +* .. Local Scalars .. + LOGICAL LEFT, RIGHT, TRAN, NOTRAN + INTEGER I, IB, MB, LB, KF, Q +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, DLARFB +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* .. Test the input arguments .. +* + INFO = 0 + LEFT = LSAME( SIDE, 'L' ) + RIGHT = LSAME( SIDE, 'R' ) + TRAN = LSAME( TRANS, 'T' ) + NOTRAN = LSAME( TRANS, 'N' ) +* + IF( LEFT ) THEN + Q = M + ELSE IF ( RIGHT ) THEN + Q = N + END IF + IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN + INFO = -1 + ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( K.LT.0 ) THEN + INFO = -5 + ELSE IF( L.LT.0 .OR. L.GT.K ) THEN + INFO = -6 + ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN + INFO = -7 + ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN + INFO = -8 + ELSE IF( LDT.LT.NB ) THEN + INFO = -10 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -12 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -12 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DTPMQRT', -INFO ) + RETURN + END IF +* +* .. Quick return if possible .. +* + IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN +* + IF( LEFT .AND. TRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + 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 DTPRFB( 'L', 'T', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. NOTRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL DTPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + ELSE IF( LEFT .AND. NOTRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + 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 DTPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. TRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL DTPRFB( 'R', 'T', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + END IF +* + RETURN +* +* End of DTPMQRT +* + END diff --git a/SRC/dtpqrt.f b/SRC/dtpqrt.f new file mode 100644 index 00000000..bf1aabb6 --- /dev/null +++ b/SRC/dtpqrt.f @@ -0,0 +1,185 @@ + SUBROUTINE DTPQRT( 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 .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DTPQRT computes a blocked QR factorization of a real +* "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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DTPQRT2, DTPRFB, 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( 'DTPQRT', -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 DTPQRT2( MB, IB, LB, A(I,I), LDA, B( 1, I ), LDB, + $ T(1, I ), LDT, IINFO ) +* +* Update by applying H**T to B(:,I+IB:N) from the left +* + IF( I+IB.LE.N ) THEN + CALL DTPRFB( 'L', 'T', '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 DTPQRT +* + END diff --git a/SRC/dtpqrt2.f b/SRC/dtpqrt2.f new file mode 100644 index 00000000..a6fd36b3 --- /dev/null +++ b/SRC/dtpqrt2.f @@ -0,0 +1,224 @@ + SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO ) + IMPLICIT NONE +* +* -- LAPACK routine (version 3.x) -- +* -- 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 +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ) +* .. +* +* Purpose +* ======= +* +* DTPQRT2 computes a QR factorization of a real "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 total 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. +* +* A (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDT,N) +* The N-by-N upper triangular factor T of the block reflector. +* See Further Details. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= max(1,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 (M+N)-by-(M+N) block reflector H is then given by +* +* H = I - W * T * W**T +* +* where W^H is the conjugate transpose of W and T is the upper triangular +* factor of the block reflector. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER( ONE = 1.0, ZERO = 0.0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, P, MP, NP + DOUBLE PRECISION ALPHA +* .. +* .. External Subroutines .. + EXTERNAL DLARFG, DGEMV, DGER, DTRMV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -7 + ELSE IF( LDT.LT.MAX( 1, N ) ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DTPQRT2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. M.EQ.0 ) RETURN +* + DO I = 1, N +* +* Generate elementary reflector H(I) to annihilate B(:,I) +* + P = M-L+MIN( L, I ) + CALL DLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) ) + IF( I.LT.N ) THEN +* +* W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)] +* + DO J = 1, N-I + T( J, N ) = (A( I, I+J )) + END DO + CALL DGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB, + $ B( 1, I ), 1, ONE, T( 1, N ), 1 ) +* +* C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H +* + ALPHA = -(T( I, 1 )) + DO J = 1, N-I + A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N )) + END DO + CALL DGER( P, N-I, ALPHA, B( 1, I ), 1, + $ T( 1, N ), 1, B( 1, I+1 ), LDB ) + END IF + END DO +* + DO I = 2, N +* +* T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I)) +* + ALPHA = -T( I, 1 ) + + DO J = 1, I-1 + T( J, I ) = ZERO + END DO + P = MIN( I-1, L ) + MP = MIN( M-L+1, M ) + NP = MIN( P+1, N ) +* +* Triangular part of B2 +* + DO J = 1, P + T( J, I ) = ALPHA*B( M-L+J, I ) + END DO + CALL DTRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB, + $ T( 1, I ), 1 ) +* +* Rectangular part of B2 +* + CALL DGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB, + $ B( MP, I ), 1, ZERO, T( NP, I ), 1 ) +* +* B1 +* + CALL DGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, + $ ONE, T( 1, I ), 1 ) +* +* T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I) +* + CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 ) +* +* T(I,I) = tau(I) +* + T( I, I ) = T( I, 1 ) + T( I, 1 ) = ZERO + END DO + +* +* End of DTPQRT2 +* + END diff --git a/SRC/dtprfb.f b/SRC/dtprfb.f new file mode 100644 index 00000000..c7b4d58a --- /dev/null +++ b/SRC/dtprfb.f @@ -0,0 +1,705 @@ + SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, + $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine (version 3.x) -- +* -- 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 .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ V( LDV, * ), WORK( LDWORK, * ) +* .. +* +* Purpose +* ======= +* +* DTPRFB applies a real "triangular-pentagonal" block reflector H or its +* transpose H**T to a real matrix C, which is composed of two +* blocks A and B, either from the left or right. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply H or H**T from the Left +* = 'R': apply H or H**T from the Right +* +* TRANS (input) CHARACTER*1 +* = 'N': apply H (No transpose) +* = 'T': apply H**T (Transpose) +* +* DIRECT (input) CHARACTER*1 +* Indicates how H is formed from a product of elementary +* reflectors +* = 'F': H = H(1) H(2) . . . H(k) (Forward) +* = 'B': H = H(k) . . . H(2) H(1) (Backward) +* +* STOREV (input) CHARACTER*1 +* Indicates how the vectors which define the elementary +* reflectors are stored: +* = 'C': Columns +* = 'R': Rows +* +* M (input) INTEGER +* The number of rows of the matrix B. +* M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. +* N >= 0. +* +* K (input) INTEGER +* The order of the matrix T, i.e. the number of elementary +* reflectors whose product defines the block reflector. +* K >= 0. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* V (input) DOUBLE PRECISION array, dimension +* (LDV,K) if STOREV = 'C' +* (LDV,M) if STOREV = 'R' and SIDE = 'L' +* (LDV,N) if STOREV = 'R' and SIDE = 'R' +* The pentagonal matrix V, which contains the elementary reflectors +* H(1), H(2), ..., H(K). See Further Details. +* +* LDV (input) INTEGER +* The leading dimension of the array V. +* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +* if STOREV = 'R', LDV >= K. +* +* T (input) DOUBLE PRECISION array, dimension (LDT,K) +* The triangular K-by-K matrix T in the representation of the +* block reflector. +* +* LDT (input) INTEGER +* The leading dimension of the array T. +* LDT >= K. +* +* A (input/output) DOUBLE PRECISION array, dimension +* (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* H*C or H**T*C or C*H or C*H**T. See Futher Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) DOUBLE PRECISION array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* H*C or H**T*C or C*H or C*H**T. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace) DOUBLE PRECISION array, dimension +* (LDWORK,N) if SIDE = 'L', +* (LDWORK,K) if SIDE = 'R'. +* +* LDWORK (input) INTEGER +* The leading dimension of the array WORK. +* If SIDE = 'L', LDWORK >= K; +* if SIDE = 'R', LDWORK >= M. +* +* Further Details +* =============== +* +* The matrix C is a composite matrix formed from blocks A and B. +* The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, +* and if SIDE = 'L', A is of size K-by-N. +* +* If SIDE = 'R' and DIRECT = 'F', C = [A B]. +* +* If SIDE = 'L' and DIRECT = 'F', C = [A] +* [B]. +* +* If SIDE = 'R' and DIRECT = 'B', C = [B A]. +* +* If SIDE = 'L' and DIRECT = 'B', C = [B] +* [A]. +* +* The pentagonal matrix V is composed of a rectangular block V1 and a +* trapezoidal block V2. The size of the trapezoidal block is determined by +* the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; +* if L=0, there is no trapezoidal block, thus V = V1 is rectangular. +* +* If DIRECT = 'F' and STOREV = 'C': V = [V1] +* [V2] +* - V2 is upper trapezoidal (first L rows of K-by-K upper triangular) +* +* If DIRECT = 'F' and STOREV = 'R': V = [V1 V2] +* +* - V2 is lower trapezoidal (first L columns of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'C': V = [V2] +* [V1] +* - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] +* +* - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) +* +* If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. +* +* If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K. +* +* If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L. +* +* If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L. +* +* ========================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0, ZERO = 0.0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, MP, NP, KP + LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DGEMM, DTRMM +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN +* + IF( LSAME( STOREV, 'C' ) ) THEN + COLUMN = .TRUE. + ROW = .FALSE. + ELSE IF ( LSAME( STOREV, 'R' ) ) THEN + COLUMN = .FALSE. + ROW = .TRUE. + ELSE + COLUMN = .FALSE. + ROW = .FALSE. + END IF +* + IF( LSAME( SIDE, 'L' ) ) THEN + LEFT = .TRUE. + RIGHT = .FALSE. + ELSE IF( LSAME( SIDE, 'R' ) ) THEN + LEFT = .FALSE. + RIGHT = .TRUE. + ELSE + LEFT = .FALSE. + RIGHT = .FALSE. + END IF +* + IF( LSAME( DIRECT, 'F' ) ) THEN + FORWARD = .TRUE. + BACKWARD = .FALSE. + ELSE IF( LSAME( DIRECT, 'B' ) ) THEN + FORWARD = .FALSE. + BACKWARD = .TRUE. + ELSE + FORWARD = .FALSE. + BACKWARD = .FALSE. + END IF +* +* --------------------------------------------------------------------------- +* + IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (M-by-K) +* +* Form H C or H**T C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W T W**T or H**T = I - W T**T W**T +* +* A = A - T (A + V**T B) or A = A - T**T (A + V**T B) +* B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, + $ ONE, WORK, LDWORK ) + CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (N-by-K) +* +* Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W T W**T or H**T = I - W T**T W**T +* +* A = A - (A + B V) T or A = A - (A + B V) T**T +* B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, + $ V, LDV, ONE, WORK, LDWORK ) + CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB ) + CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (M-by-K) +* [ I ] (K-by-K) +* +* Form H C or H**T C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W T W**T or H**T = I - W T**T W**T +* +* A = A - T (A + V**T B) or A = A - T**T (A + V**T B) +* B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO +* + CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV, + $ B, LDB, ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (N-by-K) +* [ I ] (K-by-K) +* +* Form C H or C H**T where C = [ B A ] (B is M-by-N, A is M-by-K) +* +* H = I - W T W**T or H**T = I - W T**T W**T +* +* A = A - (A + B V) T or A = A - (A + B V) T**T +* B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V, LDV, ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK, + $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB ) + CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H**T C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W**T T W or H**T = I - W**T T**T W +* +* A = A - T (A + V B) or A = A - T**T (A + V B) +* B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDB ) + CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, + $ ONE, WORK, LDWORK ) + CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W**T T W or H**T = I - W**T T**T W +* +* A = A - (A + B V**T) T or A = A - (A + B V**T) T**T +* B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV, + $ ONE, WORK, LDWORK ) + CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, + $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H**T C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W**T T W or H**T = I - W**T T**T W +* +* A = A - T (A + V B) or A = A - T**T (A + V B) +* B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO + CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H**T where C = [ B A ] (A is M-by-K, B is M-by-N) +* +* H = I - W**T T W or H**T = I - W**T T**T W +* +* A = A - (A + B V**T) T or A = A - (A + B V**T) T**T +* B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* + END IF +* + RETURN +* +* End of DTPRFB +* + END diff --git a/SRC/stpmqrt.f b/SRC/stpmqrt.f new file mode 100644 index 00000000..c6c2c858 --- /dev/null +++ b/SRC/stpmqrt.f @@ -0,0 +1,264 @@ + SUBROUTINE STPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, + $ A, LDA, B, LDB, 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 .. + CHARACTER SIDE, TRANS + INTEGER INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT +* .. +* .. Array Arguments .. + REAL V( LDV, * ), A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ WORK( * ) +* .. +* +* Purpose +* ======= +* +* STPMQRT applies a real orthogonal matrix Q obtained from a +* "triangular-pentagonal" real block reflector H to a general +* real matrix C, which consists of two blocks A and B. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply Q or Q^H from the Left; +* = 'R': apply Q or Q^H from the Right. +* +* TRANS (input) CHARACTER*1 +* = 'N': No transpose, apply Q; +* = 'C': Transpose, apply Q^H. +* +* M (input) INTEGER +* The number of rows of the matrix B. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. N >= 0. +* +* K (input) INTEGER +* The number of elementary reflectors whose product defines +* the matrix Q. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* NB (input) INTEGER +* The block size used for the storage of T. K >= NB >= 1. +* This must be the same value of NB used to generate T +* in CTPQRT. +* +* V (input) REAL array, dimension (LDA,K) +* The i-th column must contain the vector which defines the +* elementary reflector H(i), for i = 1,2,...,k, as returned by +* CTPQRT in B. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDA >= max(1,M); +* if SIDE = 'R', LDA >= max(1,N). +* +* T (input) REAL array, dimension (LDT,K) +* The upper triangular factors of the block reflectors +* as returned by CTPQRT, stored as a NB-by-K matrix. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= NB. +* +* A (input/output) REAL array, dimension +* (LDA,N) if SIDE = 'L' or +* (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* Q*C or Q^H*C or C*Q or C*Q^H. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) REAL array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* Q*C or Q^H*C or C*Q or C*Q^H. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace/output) REAL array. The dimension of WORK is +* N*NB if SIDE = 'L', or M*NB if SIDE = 'R'. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The columns of the pentagonal matrix V contain the elementary reflectors +* H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a +* trapezoidal block V2: +* +* V = [V1] +* [V2]. +* +* The size of the trapezoidal block V2 is determined by the parameter L, +* where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L +* rows of a K-by-K upper triangular matrix. If L=K, V2 is upper triangular; +* if L=0, there is no trapezoidal block, hence V = V1 is rectangular. +* +* If SIDE = 'L': C = [A] where A is K-by-N, B is M-by-N and V is M-by-K. +* [B] +* +* If SIDE = 'R': C = [A B] where A is M-by-K, B is M-by-N and V is N-by-K. +* +* The real orthogonal matrix Q is formed from V and T. +* +* If TRANS='N' and SIDE='L', C is on exit replaced with Q * C. +* +* If TRANS='C' and SIDE='L', C is on exit replaced with Q^H * C. +* +* If TRANS='N' and SIDE='R', C is on exit replaced with C * Q. +* +* If TRANS='C' and SIDE='R', C is on exit replaced with C * Q^H. +* +* ===================================================================== +* +* .. +* .. Local Scalars .. + LOGICAL LEFT, RIGHT, TRAN, NOTRAN + INTEGER I, IB, MB, LB, KF, Q +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, SLARFB +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* .. Test the input arguments .. +* + INFO = 0 + LEFT = LSAME( SIDE, 'L' ) + RIGHT = LSAME( SIDE, 'R' ) + TRAN = LSAME( TRANS, 'T' ) + NOTRAN = LSAME( TRANS, 'N' ) +* + IF( LEFT ) THEN + Q = M + ELSE IF ( RIGHT ) THEN + Q = N + END IF + IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN + INFO = -1 + ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( K.LT.0 ) THEN + INFO = -5 + ELSE IF( L.LT.0 .OR. L.GT.K ) THEN + INFO = -6 + ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN + INFO = -7 + ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN + INFO = -8 + ELSE IF( LDT.LT.NB ) THEN + INFO = -10 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -12 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -12 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'STPMQRT', -INFO ) + RETURN + END IF +* +* .. Quick return if possible .. +* + IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN +* + IF( LEFT .AND. TRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + 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 STPRFB( 'L', 'T', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. NOTRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL STPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + ELSE IF( LEFT .AND. NOTRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + 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 STPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. TRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL STPRFB( 'R', 'T', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + END IF +* + RETURN +* +* End of STPMQRT +* + END diff --git a/SRC/stpqrt.f b/SRC/stpqrt.f new file mode 100644 index 00000000..07b9a9a8 --- /dev/null +++ b/SRC/stpqrt.f @@ -0,0 +1,185 @@ + SUBROUTINE STPQRT( 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 .. + REAL A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* STPQRT computes a blocked QR factorization of a real +* "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) REAL 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) REAL 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) REAL 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) REAL 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 STPQRT2, STPRFB, 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( 'STPQRT', -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 STPQRT2( 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 STPRFB( 'L', 'T', '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 STPQRT +* + END diff --git a/SRC/stpqrt2.f b/SRC/stpqrt2.f new file mode 100644 index 00000000..edbed209 --- /dev/null +++ b/SRC/stpqrt2.f @@ -0,0 +1,224 @@ + SUBROUTINE STPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO ) + IMPLICIT NONE +* +* -- LAPACK routine (version 3.x) -- +* -- 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 +* .. +* .. Array Arguments .. + REAL A( LDA, * ), B( LDB, * ), T( LDT, * ) +* .. +* +* Purpose +* ======= +* +* STPQRT2 computes a QR factorization of a real "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 total 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. +* +* A (input/output) REAL 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) REAL 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) REAL array, dimension (LDT,N) +* The N-by-N upper triangular factor T of the block reflector. +* See Further Details. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= max(1,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 (M+N)-by-(M+N) block reflector H is then given by +* +* H = I - W * T * W^H +* +* where W^H is the conjugate transpose of W and T is the upper triangular +* factor of the block reflector. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER( ONE = 1.0, ZERO = 0.0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, P, MP, NP + REAL ALPHA +* .. +* .. External Subroutines .. + EXTERNAL SLARFG, SGEMV, SGER, STRMV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -7 + ELSE IF( LDT.LT.MAX( 1, N ) ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'STPQRT2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. M.EQ.0 ) RETURN +* + DO I = 1, N +* +* Generate elementary reflector H(I) to annihilate B(:,I) +* + P = M-L+MIN( L, I ) + CALL SLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) ) + IF( I.LT.N ) THEN +* +* W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)] +* + DO J = 1, N-I + T( J, N ) = (A( I, I+J )) + END DO + CALL SGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB, + $ B( 1, I ), 1, ONE, T( 1, N ), 1 ) +* +* C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H +* + ALPHA = -(T( I, 1 )) + DO J = 1, N-I + A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N )) + END DO + CALL SGER( P, N-I, ALPHA, B( 1, I ), 1, + $ T( 1, N ), 1, B( 1, I+1 ), LDB ) + END IF + END DO +* + DO I = 2, N +* +* T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I)) +* + ALPHA = -T( I, 1 ) + + DO J = 1, I-1 + T( J, I ) = ZERO + END DO + P = MIN( I-1, L ) + MP = MIN( M-L+1, M ) + NP = MIN( P+1, N ) +* +* Triangular part of B2 +* + DO J = 1, P + T( J, I ) = ALPHA*B( M-L+J, I ) + END DO + CALL STRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB, + $ T( 1, I ), 1 ) +* +* Rectangular part of B2 +* + CALL SGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB, + $ B( MP, I ), 1, ZERO, T( NP, I ), 1 ) +* +* B1 +* + CALL SGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, + $ ONE, T( 1, I ), 1 ) +* +* T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I) +* + CALL STRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 ) +* +* T(I,I) = tau(I) +* + T( I, I ) = T( I, 1 ) + T( I, 1 ) = ZERO + END DO + +* +* End of STPQRT2 +* + END diff --git a/SRC/stprfb.f b/SRC/stprfb.f new file mode 100644 index 00000000..234b8790 --- /dev/null +++ b/SRC/stprfb.f @@ -0,0 +1,705 @@ + SUBROUTINE STPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, + $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine (version 3.x) -- +* -- 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 .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ V( LDV, * ), WORK( LDWORK, * ) +* .. +* +* Purpose +* ======= +* +* STPRFB applies a real "triangular-pentagonal" block reflector H or its +* conjugate transpose H^H to a real matrix C, which is composed of two +* blocks A and B, either from the left or right. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply H or H^H from the Left +* = 'R': apply H or H^H from the Right +* +* TRANS (input) CHARACTER*1 +* = 'N': apply H (No transpose) +* = 'C': apply H^H (Conjugate transpose) +* +* DIRECT (input) CHARACTER*1 +* Indicates how H is formed from a product of elementary +* reflectors +* = 'F': H = H(1) H(2) . . . H(k) (Forward) +* = 'B': H = H(k) . . . H(2) H(1) (Backward) +* +* STOREV (input) CHARACTER*1 +* Indicates how the vectors which define the elementary +* reflectors are stored: +* = 'C': Columns +* = 'R': Rows +* +* M (input) INTEGER +* The number of rows of the matrix B. +* M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. +* N >= 0. +* +* K (input) INTEGER +* The order of the matrix T, i.e. the number of elementary +* reflectors whose product defines the block reflector. +* K >= 0. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* V (input) REAL array, dimension +* (LDV,K) if STOREV = 'C' +* (LDV,M) if STOREV = 'R' and SIDE = 'L' +* (LDV,N) if STOREV = 'R' and SIDE = 'R' +* The pentagonal matrix V, which contains the elementary reflectors +* H(1), H(2), ..., H(K). See Further Details. +* +* LDV (input) INTEGER +* The leading dimension of the array V. +* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +* if STOREV = 'R', LDV >= K. +* +* T (input) REAL array, dimension (LDT,K) +* The triangular K-by-K matrix T in the representation of the +* block reflector. +* +* LDT (input) INTEGER +* The leading dimension of the array T. +* LDT >= K. +* +* A (input/output) REAL array, dimension +* (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* H*C or H^H*C or C*H or C*H^H. See Futher Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) REAL array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* H*C or H^H*C or C*H or C*H^H. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace) REAL array, dimension +* (LDWORK,N) if SIDE = 'L', +* (LDWORK,K) if SIDE = 'R'. +* +* LDWORK (input) INTEGER +* The leading dimension of the array WORK. +* If SIDE = 'L', LDWORK >= K; +* if SIDE = 'R', LDWORK >= M. +* +* Further Details +* =============== +* +* The matrix C is a composite matrix formed from blocks A and B. +* The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, +* and if SIDE = 'L', A is of size K-by-N. +* +* If SIDE = 'R' and DIRECT = 'F', C = [A B]. +* +* If SIDE = 'L' and DIRECT = 'F', C = [A] +* [B]. +* +* If SIDE = 'R' and DIRECT = 'B', C = [B A]. +* +* If SIDE = 'L' and DIRECT = 'B', C = [B] +* [A]. +* +* The pentagonal matrix V is composed of a rectangular block V1 and a +* trapezoidal block V2. The size of the trapezoidal block is determined by +* the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; +* if L=0, there is no trapezoidal block, thus V = V1 is rectangular. +* +* If DIRECT = 'F' and STOREV = 'C': V = [V1] +* [V2] +* - V2 is upper trapezoidal (first L rows of K-by-K upper triangular) +* +* If DIRECT = 'F' and STOREV = 'R': V = [V1 V2] +* +* - V2 is lower trapezoidal (first L columns of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'C': V = [V2] +* [V1] +* - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] +* +* - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) +* +* If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. +* +* If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K. +* +* If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L. +* +* If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L. +* +* ========================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0, ZERO = 0.0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, MP, NP, KP + LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SGEMM, SLACGV, STRMM +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN +* + IF( LSAME( STOREV, 'C' ) ) THEN + COLUMN = .TRUE. + ROW = .FALSE. + ELSE IF ( LSAME( STOREV, 'R' ) ) THEN + COLUMN = .FALSE. + ROW = .TRUE. + ELSE + COLUMN = .FALSE. + ROW = .FALSE. + END IF +* + IF( LSAME( SIDE, 'L' ) ) THEN + LEFT = .TRUE. + RIGHT = .FALSE. + ELSE IF( LSAME( SIDE, 'R' ) ) THEN + LEFT = .FALSE. + RIGHT = .TRUE. + ELSE + LEFT = .FALSE. + RIGHT = .FALSE. + END IF +* + IF( LSAME( DIRECT, 'F' ) ) THEN + FORWARD = .TRUE. + BACKWARD = .FALSE. + ELSE IF( LSAME( DIRECT, 'B' ) ) THEN + FORWARD = .FALSE. + BACKWARD = .TRUE. + ELSE + FORWARD = .FALSE. + BACKWARD = .FALSE. + END IF +* +* --------------------------------------------------------------------------- +* + IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (M-by-K) +* +* Form H C or H^H C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W T W^H or H^H = I - W T^H W^H +* +* A = A - T (A + V^H B) or A = A - T^H (A + V^H B) +* B = B - V T (A + V^H B) or B = B - V T^H (A + V^H B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL STRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + CALL SGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, + $ ONE, WORK, LDWORK ) + CALL SGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL SGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL STRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (N-by-K) +* +* Form C H or C H^H where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W T W^H or H^H = I - W T^H W^H +* +* A = A - (A + B V) T or A = A - (A + B V) T^H +* B = B - (A + B V) T V^H or B = B - (A + B V) T^H V^H +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL STRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + CALL SGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, + $ V, LDV, ONE, WORK, LDWORK ) + CALL SGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL SGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB ) + CALL STRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (M-by-K) +* [ I ] (K-by-K) +* +* Form H C or H^H C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W T W^H or H^H = I - W T^H W^H +* +* A = A - T (A + V^H B) or A = A - T^H (A + V^H B) +* B = B - V T (A + V^H B) or B = B - V T^H (A + V^H B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO +* + CALL STRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL SGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL SGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV, + $ B, LDB, ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL SGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL STRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (N-by-K) +* [ I ] (K-by-K) +* +* Form C H or C H^H where C = [ B A ] (B is M-by-N, A is M-by-K) +* +* H = I - W T W^H or H^H = I - W T^H W^H +* +* A = A - (A + B V) T or A = A - (A + B V) T^H +* B = B - (A + B V) T V^H or B = B - (A + B V) T^H V^H +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL STRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL SGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL SGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V, LDV, ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK, + $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB ) + CALL SGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL STRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H^H C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W^H T W or H^H = I - W^H T^H W +* +* A = A - T (A + V B) or A = A - T^H (A + V B) +* B = B - V^H T (A + V B) or B = B - V^H T^H (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL STRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDB ) + CALL SGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, + $ ONE, WORK, LDWORK ) + CALL SGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL SGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL STRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H^H where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W^H T W or H^H = I - W^H T^H W +* +* A = A - (A + B V^H) T or A = A - (A + B V^H) T^H +* B = B - (A + B V^H) T V or B = B - (A + B V^H) T^H V +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL STRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + CALL SGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV, + $ ONE, WORK, LDWORK ) + CALL SGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, + $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL SGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL STRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H^H C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W^H T W or H^H = I - W^H T^H W +* +* A = A - T (A + V B) or A = A - T^H (A + V B) +* B = B - V^H T (A + V B) or B = B - V^H T^H (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO + CALL STRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL SGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL SGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL SGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL STRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H^H where C = [ B A ] (A is M-by-K, B is M-by-N) +* +* H = I - W^H T W or H^H = I - W^H T^H W +* +* A = A - (A + B V^H) T or A = A - (A + B V^H) T^H +* B = B - (A + B V^H) T V or B = B - (A + B V^H) T^H V +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL STRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL SGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL SGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL STRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL SGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL SGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL STRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* + END IF +* + RETURN +* +* End of STPRFB +* + END diff --git a/SRC/ztpmqrt.f b/SRC/ztpmqrt.f new file mode 100644 index 00000000..ffb0ce25 --- /dev/null +++ b/SRC/ztpmqrt.f @@ -0,0 +1,264 @@ + SUBROUTINE ZTPMQRT( SIDE, TRANS, M, N, K, L, NB, V, LDV, T, LDT, + $ A, LDA, B, LDB, 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 .. + CHARACTER SIDE, TRANS + INTEGER INFO, K, LDV, LDA, LDB, M, N, L, NB, LDT +* .. +* .. Array Arguments .. + COMPLEX*16 V( LDV, * ), A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZTPMQRT applies a complex orthogonal matrix Q obtained from a +* "triangular-pentagonal" complex block reflector H to a general +* complex matrix C, which consists of two blocks A and B. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply Q or Q**H from the Left; +* = 'R': apply Q or Q**H from the Right. +* +* TRANS (input) CHARACTER*1 +* = 'N': No transpose, apply Q; +* = 'C': Transpose, apply Q**H. +* +* M (input) INTEGER +* The number of rows of the matrix B. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. N >= 0. +* +* K (input) INTEGER +* The number of elementary reflectors whose product defines +* the matrix Q. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* NB (input) INTEGER +* The block size used for the storage of T. K >= NB >= 1. +* This must be the same value of NB used to generate T +* in CTPQRT. +* +* V (input) COMPLEX*16 array, dimension (LDA,K) +* The i-th column must contain the vector which defines the +* elementary reflector H(i), for i = 1,2,...,k, as returned by +* CTPQRT in B. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDA >= max(1,M); +* if SIDE = 'R', LDA >= max(1,N). +* +* T (input) COMPLEX*16 array, dimension (LDT,K) +* The upper triangular factors of the block reflectors +* as returned by CTPQRT, stored as a NB-by-K matrix. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= NB. +* +* A (input/output) COMPLEX*16 array, dimension +* (LDA,N) if SIDE = 'L' or +* (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* Q*C or Q**H*C or C*Q or C*Q**H. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) COMPLEX*16 array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* Q*C or Q**H*C or C*Q or C*Q**H. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace/output) COMPLEX*16 array. The dimension of WORK is +* N*NB if SIDE = 'L', or M*NB if SIDE = 'R'. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* Further Details +* =============== +* +* The columns of the pentagonal matrix V contain the elementary reflectors +* H(1), H(2), ..., H(K); V is composed of a rectangular block V1 and a +* trapezoidal block V2: +* +* V = [V1] +* [V2]. +* +* The size of the trapezoidal block V2 is determined by the parameter L, +* where 0 <= L <= K; V2 is upper trapezoidal, consisting of the first L +* rows of a K-by-K upper triangular matrix. If L=K, V2 is upper triangular; +* if L=0, there is no trapezoidal block, hence V = V1 is rectangular. +* +* If SIDE = 'L': C = [A] where A is K-by-N, B is M-by-N and V is M-by-K. +* [B] +* +* If SIDE = 'R': C = [A B] where A is M-by-K, B is M-by-N and V is N-by-K. +* +* The complex orthogonal matrix Q is formed from V and T. +* +* If TRANS='N' and SIDE='L', C is on exit replaced with Q * C. +* +* If TRANS='C' and SIDE='L', C is on exit replaced with Q**H * C. +* +* If TRANS='N' and SIDE='R', C is on exit replaced with C * Q. +* +* If TRANS='C' and SIDE='R', C is on exit replaced with C * Q**H. +* +* ===================================================================== +* +* .. +* .. Local Scalars .. + LOGICAL LEFT, RIGHT, TRAN, NOTRAN + INTEGER I, IB, MB, LB, KF, Q +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZLARFB +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* .. Test the input arguments .. +* + INFO = 0 + LEFT = LSAME( SIDE, 'L' ) + RIGHT = LSAME( SIDE, 'R' ) + TRAN = LSAME( TRANS, 'C' ) + NOTRAN = LSAME( TRANS, 'N' ) +* + IF( LEFT ) THEN + Q = M + ELSE IF ( RIGHT ) THEN + Q = N + END IF + IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN + INFO = -1 + ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( K.LT.0 ) THEN + INFO = -5 + ELSE IF( L.LT.0 .OR. L.GT.K ) THEN + INFO = -6 + ELSE IF( NB.LT.1 .OR. NB.GT.K ) THEN + INFO = -7 + ELSE IF( LDV.LT.MAX( 1, Q ) ) THEN + INFO = -8 + ELSE IF( LDT.LT.NB ) THEN + INFO = -10 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -12 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -12 + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZTPMQRT', -INFO ) + RETURN + END IF +* +* .. Quick return if possible .. +* + IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) RETURN +* + IF( LEFT .AND. TRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + 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 ZTPRFB( 'L', 'C', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. NOTRAN ) THEN +* + DO I = 1, K, NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL ZTPRFB( 'R', 'N', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + ELSE IF( LEFT .AND. NOTRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + 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 ZTPRFB( 'L', 'N', 'F', 'C', MB, N, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( I, 1 ), LDA, B, LDB, WORK, IB ) + END DO +* + ELSE IF( RIGHT .AND. TRAN ) THEN +* + KF = ((K-1)/NB)*NB+1 + DO I = KF, 1, -NB + IB = MIN( NB, K-I+1 ) + MB = MIN( N-L+I+IB-1, N ) + IF( I.GE.L ) THEN + LB = 0 + ELSE + LB = MB-N+L-I+1 + END IF + CALL ZTPRFB( 'R', 'C', 'F', 'C', M, MB, IB, LB, + $ V( 1, I ), LDV, T( 1, I ), LDT, + $ A( 1, I ), LDA, B, LDB, WORK, M ) + END DO +* + END IF +* + RETURN +* +* End of ZTPMQRT +* + END diff --git a/SRC/ztpqrt.f b/SRC/ztpqrt.f new file mode 100644 index 00000000..d0570407 --- /dev/null +++ b/SRC/ztpqrt.f @@ -0,0 +1,185 @@ + SUBROUTINE ZTPQRT( 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*16 A( LDA, * ), B( LDB, * ), T( LDT, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* ZTPQRT 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*16 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*16 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*16 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*16 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 ZTPQRT2, ZTPRFB, 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( 'ZTPQRT', -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 ZTPQRT2( 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 ZTPRFB( '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 ZTPQRT +* + END diff --git a/SRC/ztpqrt2.f b/SRC/ztpqrt2.f new file mode 100644 index 00000000..0d087fd6 --- /dev/null +++ b/SRC/ztpqrt2.f @@ -0,0 +1,224 @@ + SUBROUTINE ZTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO ) + IMPLICIT NONE +* +* -- LAPACK routine (version 3.x) -- +* -- 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 +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ) +* .. +* +* Purpose +* ======= +* +* ZTPQRT2 computes a 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 total 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. +* +* A (input/output) COMPLEX*16 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*16 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*16 array, dimension (LDT,N) +* The N-by-N upper triangular factor T of the block reflector. +* See Further Details. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= max(1,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 (M+N)-by-(M+N) block reflector H is then given by +* +* H = I - W * T * W**H +* +* where W**H is the conjugate transpose of W and T is the upper triangular +* factor of the block reflector. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) ) +* .. +* .. Local Scalars .. + INTEGER I, J, P, MP, NP + COMPLEX*16 ALPHA +* .. +* .. External Subroutines .. + EXTERNAL ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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( LDA.LT.MAX( 1, N ) ) THEN + INFO = -5 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -7 + ELSE IF( LDT.LT.MAX( 1, N ) ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZTPQRT2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 .OR. M.EQ.0 ) RETURN +* + DO I = 1, N +* +* Generate elementary reflector H(I) to annihilate B(:,I) +* + P = M-L+MIN( L, I ) + CALL ZLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) ) + IF( I.LT.N ) THEN +* +* W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)] +* + DO J = 1, N-I + T( J, N ) = CONJG(A( I, I+J )) + END DO + CALL ZGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB, + $ B( 1, I ), 1, ONE, T( 1, N ), 1 ) +* +* C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H +* + ALPHA = -CONJG(T( I, 1 )) + DO J = 1, N-I + A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N )) + END DO + CALL ZGERC( P, N-I, ALPHA, B( 1, I ), 1, + $ T( 1, N ), 1, B( 1, I+1 ), LDB ) + END IF + END DO +* + DO I = 2, N +* +* T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I)) +* + ALPHA = -T( I, 1 ) + + DO J = 1, I-1 + T( J, I ) = ZERO + END DO + P = MIN( I-1, L ) + MP = MIN( M-L+1, M ) + NP = MIN( P+1, N ) +* +* Triangular part of B2 +* + DO J = 1, P + T( J, I ) = ALPHA*B( M-L+J, I ) + END DO + CALL ZTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB, + $ T( 1, I ), 1 ) +* +* Rectangular part of B2 +* + CALL ZGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB, + $ B( MP, I ), 1, ZERO, T( NP, I ), 1 ) +* +* B1 +* + CALL ZGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1, + $ ONE, T( 1, I ), 1 ) +* +* T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I) +* + CALL ZTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 ) +* +* T(I,I) = tau(I) +* + T( I, I ) = T( I, 1 ) + T( I, 1 ) = ZERO + END DO + +* +* End of ZTPQRT2 +* + END diff --git a/SRC/ztprfb.f b/SRC/ztprfb.f new file mode 100644 index 00000000..c97c35d1 --- /dev/null +++ b/SRC/ztprfb.f @@ -0,0 +1,708 @@ + SUBROUTINE ZTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, + $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine (version 3.x) -- +* -- 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 .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ), + $ V( LDV, * ), WORK( LDWORK, * ) +* .. +* +* Purpose +* ======= +* +* ZTPRFB applies a complex "triangular-pentagonal" block reflector H or its +* conjugate transpose H**H to a complex matrix C, which is composed of two +* blocks A and B, either from the left or right. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply H or H**H from the Left +* = 'R': apply H or H**H from the Right +* +* TRANS (input) CHARACTER*1 +* = 'N': apply H (No transpose) +* = 'C': apply H**H (Conjugate transpose) +* +* DIRECT (input) CHARACTER*1 +* Indicates how H is formed from a product of elementary +* reflectors +* = 'F': H = H(1) H(2) . . . H(k) (Forward) +* = 'B': H = H(k) . . . H(2) H(1) (Backward) +* +* STOREV (input) CHARACTER*1 +* Indicates how the vectors which define the elementary +* reflectors are stored: +* = 'C': Columns +* = 'R': Rows +* +* M (input) INTEGER +* The number of rows of the matrix B. +* M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix B. +* N >= 0. +* +* K (input) INTEGER +* The order of the matrix T, i.e. the number of elementary +* reflectors whose product defines the block reflector. +* K >= 0. +* +* L (input) INTEGER +* The order of the trapezoidal part of V. +* K >= L >= 0. See Further Details. +* +* V (input) COMPLEX*16 array, dimension +* (LDV,K) if STOREV = 'C' +* (LDV,M) if STOREV = 'R' and SIDE = 'L' +* (LDV,N) if STOREV = 'R' and SIDE = 'R' +* The pentagonal matrix V, which contains the elementary reflectors +* H(1), H(2), ..., H(K). See Further Details. +* +* LDV (input) INTEGER +* The leading dimension of the array V. +* If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +* if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +* if STOREV = 'R', LDV >= K. +* +* T (input) COMPLEX*16 array, dimension (LDT,K) +* The triangular K-by-K matrix T in the representation of the +* block reflector. +* +* LDT (input) INTEGER +* The leading dimension of the array T. +* LDT >= K. +* +* A (input/output) COMPLEX*16 array, dimension +* (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' +* On entry, the K-by-N or M-by-K matrix A. +* On exit, A is overwritten by the corresponding block of +* H*C or H**H*C or C*H or C*H**H. See Futher Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDC >= max(1,K); +* If SIDE = 'R', LDC >= max(1,M). +* +* B (input/output) COMPLEX*16 array, dimension (LDB,N) +* On entry, the M-by-N matrix B. +* On exit, B is overwritten by the corresponding block of +* H*C or H**H*C or C*H or C*H**H. See Further Details. +* +* LDB (input) INTEGER +* The leading dimension of the array B. +* LDB >= max(1,M). +* +* WORK (workspace) COMPLEX*16 array, dimension +* (LDWORK,N) if SIDE = 'L', +* (LDWORK,K) if SIDE = 'R'. +* +* LDWORK (input) INTEGER +* The leading dimension of the array WORK. +* If SIDE = 'L', LDWORK >= K; +* if SIDE = 'R', LDWORK >= M. +* +* Further Details +* =============== +* +* The matrix C is a composite matrix formed from blocks A and B. +* The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, +* and if SIDE = 'L', A is of size K-by-N. +* +* If SIDE = 'R' and DIRECT = 'F', C = [A B]. +* +* If SIDE = 'L' and DIRECT = 'F', C = [A] +* [B]. +* +* If SIDE = 'R' and DIRECT = 'B', C = [B A]. +* +* If SIDE = 'L' and DIRECT = 'B', C = [B] +* [A]. +* +* The pentagonal matrix V is composed of a rectangular block V1 and a +* trapezoidal block V2. The size of the trapezoidal block is determined by +* the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; +* if L=0, there is no trapezoidal block, thus V = V1 is rectangular. +* +* If DIRECT = 'F' and STOREV = 'C': V = [V1] +* [V2] +* - V2 is upper trapezoidal (first L rows of K-by-K upper triangular) +* +* If DIRECT = 'F' and STOREV = 'R': V = [V1 V2] +* +* - V2 is lower trapezoidal (first L columns of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'C': V = [V2] +* [V1] +* - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) +* +* If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] +* +* - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) +* +* If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. +* +* If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K. +* +* If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L. +* +* If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L. +* +* ========================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,0.0) ) +* .. +* .. Local Scalars .. + INTEGER I, J, MP, NP, KP + LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL ZCOPY, ZGEMM, ZLACGV, ZTRMM +* .. +* .. Intrinsic Functions .. + INTRINSIC CONJG +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN +* + IF( LSAME( STOREV, 'C' ) ) THEN + COLUMN = .TRUE. + ROW = .FALSE. + ELSE IF ( LSAME( STOREV, 'R' ) ) THEN + COLUMN = .FALSE. + ROW = .TRUE. + ELSE + COLUMN = .FALSE. + ROW = .FALSE. + END IF +* + IF( LSAME( SIDE, 'L' ) ) THEN + LEFT = .TRUE. + RIGHT = .FALSE. + ELSE IF( LSAME( SIDE, 'R' ) ) THEN + LEFT = .FALSE. + RIGHT = .TRUE. + ELSE + LEFT = .FALSE. + RIGHT = .FALSE. + END IF +* + IF( LSAME( DIRECT, 'F' ) ) THEN + FORWARD = .TRUE. + BACKWARD = .FALSE. + ELSE IF( LSAME( DIRECT, 'B' ) ) THEN + FORWARD = .FALSE. + BACKWARD = .TRUE. + ELSE + FORWARD = .FALSE. + BACKWARD = .FALSE. + END IF +* +* --------------------------------------------------------------------------- +* + IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (M-by-K) +* +* Form H C or H**H C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - T (A + V**H B) or A = A - T**H (A + V**H B) +* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB, + $ ONE, WORK, LDWORK ) + CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I ] (K-by-K) +* [ V ] (N-by-K) +* +* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - (A + B V) T or A = A - (A + B V) T**H +* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, + $ V, LDV, ONE, WORK, LDWORK ) + CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB ) + CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( NP, 1 ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (M-by-K) +* [ I ] (K-by-K) +* +* Form H C or H**H C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - T (A + V**H B) or A = A - T**H (A + V**H B) +* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO +* + CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V, LDV, + $ B, LDB, ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V ] (N-by-K) +* [ I ] (K-by-K) +* +* Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K) +* +* H = I - W T W**H or H**H = I - W T**H W**H +* +* A = A - (A + B V) T or A = A - (A + B V) T**H +* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, + $ V, LDV, ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK, + $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB ) + CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, KP ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H**H C where C = [ A ] (K-by-N) +* [ B ] (M-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - T (A + V B) or A = A - T**H (A + V B) +* B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( M-L+1, M ) + KP = MIN( L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( I, J ) = B( M-L+I, J ) + END DO + END DO + CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDB ) + CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, + $ ONE, WORK, LDWORK ) + CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, + $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, B, LDB ) + CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, + $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) + CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, MP ), LDV, + $ WORK, LDWORK ) + DO J = 1, N + DO I = 1, L + B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ I V ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H +* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V +* +* --------------------------------------------------------------------------- +* + NP = MIN( N-L+1, N ) + KP = MIN( L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, J ) = B( I, N-L+J ) + END DO + END DO + CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV, + $ ONE, WORK, LDWORK ) + CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, + $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL ZGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, + $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, + $ WORK, LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-M ) +* +* Form H C or H**H C where C = [ B ] (M-by-N) +* [ A ] (K-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - T (A + V B) or A = A - T**H (A + V B) +* B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B) +* +* --------------------------------------------------------------------------- +* + MP = MIN( L+1, M ) + KP = MIN( K-L+1, K ) +* + DO J = 1, N + DO I = 1, L + WORK( K-L+I, J ) = B( I, J ) + END DO + END DO + CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV, + $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) + CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, N + DO I = 1, K + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, + $ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) + CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV, + $ WORK, LDWORK, ONE, B, LDB ) + CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( KP, 1 ), LDV, + $ WORK( KP, 1 ), LDWORK ) + DO J = 1, N + DO I = 1, L + B( I, J ) = B( I, J ) - WORK( K-L+I, J ) + END DO + END DO +* +* --------------------------------------------------------------------------- +* + ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN +* +* --------------------------------------------------------------------------- +* +* Let W = [ V I ] ( I is K-by-K, V is K-by-N ) +* +* Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N) +* +* H = I - W**H T W or H**H = I - W**H T**H W +* +* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H +* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V +* +* --------------------------------------------------------------------------- +* + NP = MIN( L+1, N ) + KP = MIN( K-L+1, K ) +* + DO J = 1, L + DO I = 1, M + WORK( I, K-L+J ) = B( I, J ) + END DO + END DO + CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB, + $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) + CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, V, LDV, + $ ZERO, WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + WORK( I, J ) = WORK( I, J ) + A( I, J ) + END DO + END DO +* + CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, + $ WORK, LDWORK ) +* + DO J = 1, K + DO I = 1, M + A( I, J ) = A( I, J ) - WORK( I, J ) + END DO + END DO +* + CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, + $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) + CALL ZGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, + $ V, LDV, ONE, B, LDB ) + CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, + $ WORK( 1, KP ), LDWORK ) + DO J = 1, L + DO I = 1, M + B( I, J ) = B( I, J ) - WORK( I, K-L+J ) + END DO + END DO +* + END IF +* + RETURN +* +* End of ZTPRFB +* + END |