From 4ca2feaf79883558f849f792f6813819da97a821 Mon Sep 17 00:00:00 2001 From: brian Date: Wed, 3 Nov 2010 23:02:29 +0000 Subject: Added CS decomposition source files to SRC/ --- SRC/cuncsd.f | 466 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 466 insertions(+) create mode 100644 SRC/cuncsd.f (limited to 'SRC/cuncsd.f') diff --git a/SRC/cuncsd.f b/SRC/cuncsd.f new file mode 100644 index 00000000..278564ae --- /dev/null +++ b/SRC/cuncsd.f @@ -0,0 +1,466 @@ + RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, + $ SIGNS, M, P, Q, X11, LDX11, X12, + $ LDX12, X21, LDX21, X22, LDX22, THETA, + $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, + $ LDV2T, WORK, LWORK, RWORK, LRWORK, + $ IWORK, INFO ) + IMPLICIT NONE +* +* Brian Sutton +* Randolph-Macon College +* July 2010 +* +* .. Scalar Arguments .. + CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS + INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, + $ LDX21, LDX22, LRWORK, LWORK, M, P, Q +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL THETA( * ) + REAL RWORK( * ) + COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), + $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), + $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, + $ * ) +* .. +* +* Purpose +* ======= +* +* CUNCSD computes the CS decomposition of an M-by-M partitioned +* unitary matrix X: +* +* [ I 0 0 | 0 0 0 ] +* [ 0 C 0 | 0 -S 0 ] +* [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H +* X = [-----------] = [---------] [---------------------] [---------] . +* [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ] +* [ 0 S 0 | 0 C 0 ] +* [ 0 0 I | 0 0 0 ] +* +* X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P, +* (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are +* R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in +* which R = MIN(P,M-P,Q,M-Q). +* +* Arguments +* ========= +* +* JOBU1 (input) CHARACTER +* = 'Y': U1 is computed; +* otherwise: U1 is not computed. +* +* JOBU2 (input) CHARACTER +* = 'Y': U2 is computed; +* otherwise: U2 is not computed. +* +* JOBV1T (input) CHARACTER +* = 'Y': V1T is computed; +* otherwise: V1T is not computed. +* +* JOBV2T (input) CHARACTER +* = 'Y': V2T is computed; +* otherwise: V2T is not computed. +* +* TRANS (input) CHARACTER +* = 'T': X, U1, U2, V1T, and V2T are stored in row-major +* order; +* otherwise: X, U1, U2, V1T, and V2T are stored in column- +* major order. +* +* SIGNS (input) CHARACTER +* = 'O': The lower-left block is made nonpositive (the +* "other" convention); +* otherwise: The upper-right block is made nonpositive (the +* "default" convention). +* +* M (input) INTEGER +* The number of rows and columns in X. +* +* P (input) INTEGER +* The number of rows in X11 and X12. 0 <= P <= M. +* +* Q (input) INTEGER +* The number of columns in X11 and X21. 0 <= Q <= M. +* +* X (input/workspace) COMPLEX array, dimension (LDX,M) +* On entry, the unitary matrix whose CSD is desired. +* +* LDX (input) INTEGER +* The leading dimension of X. LDX >= MAX(1,M). +* +* THETA (output) REAL array, dimension (R), in which R = +* MIN(P,M-P,Q,M-Q). +* C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and +* S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). +* +* U1 (output) COMPLEX array, dimension (P) +* If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1. +* +* LDU1 (input) INTEGER +* The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= +* MAX(1,P). +* +* U2 (output) COMPLEX array, dimension (M-P) +* If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary +* matrix U2. +* +* LDU2 (input) INTEGER +* The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= +* MAX(1,M-P). +* +* V1T (output) COMPLEX array, dimension (Q) +* If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary +* matrix V1**H. +* +* LDV1T (input) INTEGER +* The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= +* MAX(1,Q). +* +* V2T (output) COMPLEX array, dimension (M-Q) +* If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary +* matrix V2**H. +* +* LDV2T (input) INTEGER +* The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >= +* MAX(1,M-Q). +* +* WORK (workspace) COMPLEX array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the work array, and no error +* message related to LWORK is issued by XERBLA. +* +* RWORK (workspace) REAL array, dimension MAX(1,LRWORK) +* On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. +* If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1), +* ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), +* define the matrix in intermediate bidiagonal-block form +* remaining after nonconvergence. INFO specifies the number +* of nonzero PHI's. +* +* LRWORK (input) INTEGER +* The dimension of the array RWORK. +* +* If LRWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the RWORK array, returns +* this value as the first entry of the work array, and no error +* message related to LRWORK is issued by XERBLA. +* +* IWORK (workspace) INTEGER array, dimension (M-Q) +* +* INFO (output) INTEGER +* = 0: successful exit. +* < 0: if INFO = -i, the i-th argument had an illegal value. +* > 0: CBBCSD did not converge. See the description of RWORK +* above for details. +* +* Reference +* ========= +* +* [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. +* Algorithms, 50(1):33-65, 2009. +* +* =================================================================== +* +* .. Parameters .. + REAL REALONE + PARAMETER ( REALONE = 1.0E0 ) + COMPLEX NEGONE, ONE, PIOVER2, ZERO + PARAMETER ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0), + $ PIOVER2 = 1.57079632679489662E0, + $ ZERO = (0.0E0,0.0E0) ) +* .. +* .. Local Scalars .. + CHARACTER TRANST, SIGNST + INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, + $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, + $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, + $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, + $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, + $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, + $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, + $ LORGQRWORKOPT, LWORKMIN, LWORKOPT + LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, + $ WANTV1T, WANTV2T + INTEGER LRWORKMIN, LRWORKOPT + LOGICAL LRQUERY +* .. +* .. External Subroutines .. + EXTERNAL CBBCSD, CLACPY, CLAPMR, CLAPMT, CLASCL, CLASET, + $ CUNBDB, CUNGLQ, CUNGQR, XERBLA +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Intrinsic Functions + INTRINSIC COS, INT, MAX, MIN, SIN +* .. +* .. Executable Statements .. +* +* Test input arguments +* + INFO = 0 + WANTU1 = LSAME( JOBU1, 'Y' ) + WANTU2 = LSAME( JOBU2, 'Y' ) + WANTV1T = LSAME( JOBV1T, 'Y' ) + WANTV2T = LSAME( JOBV2T, 'Y' ) + COLMAJOR = .NOT. LSAME( TRANS, 'T' ) + DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) + LQUERY = LWORK .EQ. -1 + LRQUERY = LRWORK .EQ. -1 + IF( M .LT. 0 ) THEN + INFO = -7 + ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN + INFO = -8 + ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN + INFO = -9 + ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR. + $ ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN + INFO = -11 + ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN + INFO = -14 + ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN + INFO = -16 + ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN + INFO = -18 + ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN + INFO = -20 + END IF +* +* Work with transpose if convenient +* + IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN + IF( COLMAJOR ) THEN + TRANST = 'T' + ELSE + TRANST = 'N' + END IF + IF( DEFAULTSIGNS ) THEN + SIGNST = 'O' + ELSE + SIGNST = 'D' + END IF + CALL CUNCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, + $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, + $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, + $ U2, LDU2, WORK, LWORK, RWORK, LRWORK, IWORK, + $ INFO ) + RETURN + END IF +* +* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if +* convenient +* + IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN + IF( DEFAULTSIGNS ) THEN + SIGNST = 'O' + ELSE + SIGNST = 'D' + END IF + CALL CUNCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, + $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, + $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, + $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO ) + RETURN + END IF +* +* Compute workspace +* + IF( INFO .EQ. 0 ) THEN +* +* Real workspace +* + IPHI = 2 + IB11D = IPHI + MAX( 1, Q - 1 ) + IB11E = IB11D + MAX( 1, Q ) + IB12D = IB11E + MAX( 1, Q - 1 ) + IB12E = IB12D + MAX( 1, Q ) + IB21D = IB12E + MAX( 1, Q - 1 ) + IB21E = IB21D + MAX( 1, Q ) + IB22D = IB21E + MAX( 1, Q - 1 ) + IB22E = IB22D + MAX( 1, Q ) + IBBCSD = IB22E + MAX( 1, Q - 1 ) + CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0, + $ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0, + $ 0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO ) + LBBCSDWORKOPT = INT( RWORK(1) ) + LBBCSDWORKMIN = LBBCSDWORKOPT + LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1 + LRWORKMIN = IBBCSD + LBBCSDWORKMIN - 1 + RWORK(1) = LRWORKOPT +* +* Complex workspace +* + ITAUP1 = 2 + ITAUP2 = ITAUP1 + MAX( 1, P ) + ITAUQ1 = ITAUP2 + MAX( 1, M - P ) + ITAUQ2 = ITAUQ1 + MAX( 1, Q ) + IORGQR = ITAUQ2 + MAX( 1, M - Q ) + CALL CUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, + $ CHILDINFO ) + LORGQRWORKOPT = INT( WORK(1) ) + LORGQRWORKMIN = MAX( 1, M - Q ) + IORGLQ = ITAUQ2 + MAX( 1, M - Q ) + CALL CUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, + $ CHILDINFO ) + LORGLQWORKOPT = INT( WORK(1) ) + LORGLQWORKMIN = MAX( 1, M - Q ) + IORBDB = ITAUQ2 + MAX( 1, M - Q ) + CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, + $ X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK, + $ -1, CHILDINFO ) + LORBDBWORKOPT = INT( WORK(1) ) + LORBDBWORKMIN = LORBDBWORKOPT + LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, + $ IORBDB + LORBDBWORKOPT ) - 1 + LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, + $ IORBDB + LORBDBWORKMIN ) - 1 + WORK(1) = LWORKOPT +* + IF( LWORK .LT. LWORKMIN + $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN + INFO = -22 + ELSE IF( LRWORK .LT. LRWORKMIN + $ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN + INFO = -24 + ELSE + LORGQRWORK = LWORK - IORGQR + 1 + LORGLQWORK = LWORK - IORGLQ + 1 + LORBDBWORK = LWORK - IORBDB + 1 + LBBCSDWORK = LRWORK - IBBCSD + 1 + END IF + END IF +* +* Abort if any illegal arguments +* + IF( INFO .NE. 0 ) THEN + CALL XERBLA( 'CUNCSD', -INFO ) + RETURN + ELSE IF( LQUERY .OR. LRQUERY ) THEN + RETURN + END IF +* +* Transform to bidiagonal block form +* + CALL CUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, + $ LDX21, X22, LDX22, THETA, RWORK(IPHI), WORK(ITAUP1), + $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), + $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) +* +* Accumulate Householder reflectors +* + IF( COLMAJOR ) THEN + IF( WANTU1 .AND. P .GT. 0 ) THEN + CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) + CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + $ LORGQRWORK, INFO) + END IF + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL CLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) + CALL CUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), + $ WORK(IORGQR), LORGQRWORK, INFO ) + END IF + IF( WANTV1T .AND. Q .GT. 0 ) THEN + CALL CLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), + $ LDV1T ) + V1T(1, 1) = ONE + DO J = 2, Q + V1T(1,J) = ZERO + V1T(J,1) = ZERO + END DO + CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), + $ WORK(IORGLQ), LORGLQWORK, INFO ) + END IF + IF( WANTV2T .AND. M-Q .GT. 0 ) THEN + CALL CLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) + CALL CLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, + $ V2T(P+1,P+1), LDV2T ) + CALL CUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), + $ WORK(IORGLQ), LORGLQWORK, INFO ) + END IF + ELSE + IF( WANTU1 .AND. P .GT. 0 ) THEN + CALL CLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) + CALL CUNGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), + $ LORGLQWORK, INFO) + END IF + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL CLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) + CALL CUNGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), + $ WORK(IORGLQ), LORGLQWORK, INFO ) + END IF + IF( WANTV1T .AND. Q .GT. 0 ) THEN + CALL CLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), + $ LDV1T ) + V1T(1, 1) = ONE + DO J = 2, Q + V1T(1,J) = ZERO + V1T(J,1) = ZERO + END DO + CALL CUNGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), + $ WORK(IORGQR), LORGQRWORK, INFO ) + END IF + IF( WANTV2T .AND. M-Q .GT. 0 ) THEN + CALL CLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) + CALL CLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, + $ V2T(P+1,P+1), LDV2T ) + CALL CUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), + $ WORK(IORGQR), LORGQRWORK, INFO ) + END IF + END IF +* +* Compute the CSD of the matrix in bidiagonal-block form +* + CALL CBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, + $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, + $ LDV2T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), + $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), + $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), + $ LBBCSDWORK, INFO ) +* +* Permute rows and columns to place identity submatrices in top- +* left corner of (1,1)-block and/or bottom-right corner of (1,2)- +* block and/or bottom-right corner of (2,1)-block and/or top-left +* corner of (2,2)-block +* + IF( Q .GT. 0 .AND. WANTU2 ) THEN + DO I = 1, Q + IWORK(I) = M - P - Q + I + END DO + DO I = Q + 1, M - P + IWORK(I) = I - Q + END DO + IF( COLMAJOR ) THEN + CALL CLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) + ELSE + CALL CLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) + END IF + END IF + IF( M .GT. 0 .AND. WANTV2T ) THEN + DO I = 1, P + IWORK(I) = M - P - Q + I + END DO + DO I = P + 1, M - Q + IWORK(I) = I - P + END DO + IF( .NOT. COLMAJOR ) THEN + CALL CLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) + ELSE + CALL CLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) + END IF + END IF +* + RETURN +* +* End CUNCSD +* + END + -- cgit v1.2.3