diff options
36 files changed, 6460 insertions, 71 deletions
diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index eaa37d60..143fee34 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -101,7 +101,7 @@ set(SLASRC sgetrs.f sggbak.f sggbal.f sgges.f sgges3.f sggesx.f sggev.f sggev3.f sggevx.f sggglm.f sgghrd.f sgghd3.f sgglse.f sggqrf.f - sggrqf.f sggsvd.f sggsvp.f sgtcon.f sgtrfs.f sgtsv.f + sggrqf.f sggsvd.f sggsvp.f sggsvd3.f sggsvp3.f sgtcon.f sgtrfs.f sgtsv.f sgtsvx.f sgttrf.f sgttrs.f sgtts2.f shgeqz.f shsein.f shseqr.f slabrd.f slacon.f slacn2.f slaein.f slaexc.f slag2.f slags2.f slagtm.f slagv2.f slahqr.f @@ -175,7 +175,7 @@ set(CLASRC cggbak.f cggbal.f cgges.f cgges3.f cggesx.f cggev.f cggev3.f cggevx.f cggglm.f cgghrd.f cgghd3.f cgglse.f cggqrf.f cggrqf.f - cggsvd.f cggsvp.f + cggsvd.f cggsvp.f cggsvd3.f cggsvp3.f cgtcon.f cgtrfs.f cgtsv.f cgtsvx.f cgttrf.f cgttrs.f cgtts2.f chbev.f chbevd.f chbevx.f chbgst.f chbgv.f chbgvd.f chbgvx.f chbtrd.f checon.f cheev.f cheevd.f cheevr.f cheevx.f chegs2.f chegst.f @@ -259,7 +259,7 @@ set(DLASRC dgetrs.f dggbak.f dggbal.f dgges.f dgges3.f dggesx.f dggev.f dggev3.f dggevx.f dggglm.f dgghrd.f dgghd3.f dgglse.f dggqrf.f - dggrqf.f dggsvd.f dggsvp.f dgtcon.f dgtrfs.f dgtsv.f + dggrqf.f dggsvd.f dggsvp.f dggsvd3.f dggsvp3.f dgtcon.f dgtrfs.f dgtsv.f dgtsvx.f dgttrf.f dgttrs.f dgtts2.f dhgeqz.f dhsein.f dhseqr.f dlabrd.f dlacon.f dlacn2.f dlaein.f dlaexc.f dlag2.f dlags2.f dlagtm.f dlagv2.f dlahqr.f @@ -332,7 +332,7 @@ set(ZLASRC zggbak.f zggbal.f zgges.f zgges3.f zggesx.f zggev.f zggev3.f zggevx.f zggglm.f zgghrd.f zgghd3.f zgglse.f zggqrf.f zggrqf.f - zggsvd.f zggsvp.f + zggsvd.f zggsvp.f zggsvd3.f zggsvp3.f zgtcon.f zgtrfs.f zgtsv.f zgtsvx.f zgttrf.f zgttrs.f zgtts2.f zhbev.f zhbevd.f zhbevx.f zhbgst.f zhbgv.f zhbgvd.f zhbgvx.f zhbtrd.f zhecon.f zheev.f zheevd.f zheevr.f zheevx.f zhegs2.f zhegst.f diff --git a/SRC/Makefile b/SRC/Makefile index c31611e5..cb417dae 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -108,7 +108,7 @@ SLASRC = \ sggbak.o sggbal.o sgges.o sgges3.o sggesx.o \ sggev.o sggev3.o sggevx.o \ sggglm.o sgghrd.o sgghd3.o sgglse.o sggqrf.o \ - sggrqf.o sggsvd.o sggsvp.o sgtcon.o sgtrfs.o sgtsv.o \ + sggrqf.o sggsvd.o sggsvp.o sggsvd3.o sggsvp3.o sgtcon.o sgtrfs.o sgtsv.o \ sgtsvx.o sgttrf.o sgttrs.o sgtts2.o shgeqz.o \ shsein.o shseqr.o slabrd.o slacon.o slacn2.o \ slaein.o slaexc.o slag2.o slags2.o slagtm.o slagv2.o slahqr.o \ @@ -184,7 +184,7 @@ CLASRC = \ cggbak.o cggbal.o cgges.o cgges3.o cggesx.o \ cggev.o cggev3.o cggevx.o cggglm.o\ cgghrd.o cgghd3.o cgglse.o cggqrf.o cggrqf.o \ - cggsvd.o cggsvp.o \ + cggsvd.o cggsvp.o cggsvd3.o cggsvp3.o \ cgtcon.o cgtrfs.o cgtsv.o cgtsvx.o cgttrf.o cgttrs.o cgtts2.o chbev.o \ chbevd.o chbevx.o chbgst.o chbgv.o chbgvd.o chbgvx.o chbtrd.o \ checon.o cheev.o cheevd.o cheevr.o cheevx.o chegs2.o chegst.o \ @@ -270,7 +270,7 @@ DLASRC = \ dgetrs.o dggbak.o dggbal.o dgges.o dgges3.o dggesx.o \ dggev.o dggev3.o dggevx.o \ dggglm.o dgghrd.o dgghd3.o dgglse.o dggqrf.o \ - dggrqf.o dggsvd.o dggsvp.o dgtcon.o dgtrfs.o dgtsv.o \ + dggrqf.o dggsvd.o dggsvp.o dggsvd3.o dggsvp3.o dgtcon.o dgtrfs.o dgtsv.o \ dgtsvx.o dgttrf.o dgttrs.o dgtts2.o dhgeqz.o \ dhsein.o dhseqr.o dlabrd.o dlacon.o dlacn2.o \ dlaein.o dlaexc.o dlag2.o dlags2.o dlagtm.o dlagv2.o dlahqr.o \ @@ -345,7 +345,7 @@ ZLASRC = \ zggbak.o zggbal.o zgges.o zgges3.o zggesx.o \ zggev.o zggev3.o zggevx.o zggglm.o \ zgghrd.o zgghd3.o zgglse.o zggqrf.o zggrqf.o \ - zggsvd.o zggsvp.o \ + zggsvd.o zggsvp.o zggsvd3.o zggsvp3.o \ zgtcon.o zgtrfs.o zgtsv.o zgtsvx.o zgttrf.o zgttrs.o zgtts2.o zhbev.o \ zhbevd.o zhbevx.o zhbgst.o zhbgv.o zhbgvd.o zhbgvx.o zhbtrd.o \ zhecon.o zheev.o zheevd.o zheevr.o zheevx.o zhegs2.o zhegst.o \ diff --git a/SRC/cgeqp3.f b/SRC/cgeqp3.f index cc78b0b4..69f79132 100644 --- a/SRC/cgeqp3.f +++ b/SRC/cgeqp3.f @@ -220,7 +220,7 @@ NB = ILAENV( INB, 'CGEQRF', ' ', M, N, -1, -1 ) LWKOPT = ( N + 1 )*NB END IF - WORK( 1 ) = LWKOPT + WORK( 1 ) = CMPLX( LWKOPT ) * IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN INFO = -8 @@ -234,12 +234,6 @@ RETURN END IF * -* Quick return if possible. -* - IF( MINMN.EQ.0 ) THEN - RETURN - END IF -* * Move initial columns up front. * NFXD = 1 @@ -370,7 +364,7 @@ * END IF * - WORK( 1 ) = IWS + WORK( 1 ) = CMPLX( LWKOPT ) RETURN * * End of CGEQP3 diff --git a/SRC/cggsvd.f b/SRC/cggsvd.f index e7ee4087..080ef0ac 100644 --- a/SRC/cggsvd.f +++ b/SRC/cggsvd.f @@ -39,6 +39,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine CGGSVD3. +*> *> CGGSVD computes the generalized singular value decomposition (GSVD) *> of an M-by-N complex matrix A and P-by-N complex matrix B: *> diff --git a/SRC/cggsvd3.f b/SRC/cggsvd3.f new file mode 100644 index 00000000..5cef7d47 --- /dev/null +++ b/SRC/cggsvd3.f @@ -0,0 +1,506 @@ +*> \brief <b> CGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CGGSVD3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvd3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvd3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvd3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, +* LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, +* LWORK, RWORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL ALPHA( * ), BETA( * ), RWORK( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), +* $ U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CGGSVD3 computes the generalized singular value decomposition (GSVD) +*> of an M-by-N complex matrix A and P-by-N complex matrix B: +*> +*> U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R ) +*> +*> where U, V and Q are unitary matrices. +*> Let K+L = the effective numerical rank of the +*> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper +*> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" +*> matrices and of the following structures, respectively: +*> +*> If M-K-L >= 0, +*> +*> K L +*> D1 = K ( I 0 ) +*> L ( 0 C ) +*> M-K-L ( 0 0 ) +*> +*> K L +*> D2 = L ( 0 S ) +*> P-L ( 0 0 ) +*> +*> N-K-L K L +*> ( 0 R ) = K ( 0 R11 R12 ) +*> L ( 0 0 R22 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), +*> S = diag( BETA(K+1), ... , BETA(K+L) ), +*> C**2 + S**2 = I. +*> +*> R is stored in A(1:K+L,N-K-L+1:N) on exit. +*> +*> If M-K-L < 0, +*> +*> K M-K K+L-M +*> D1 = K ( I 0 0 ) +*> M-K ( 0 C 0 ) +*> +*> K M-K K+L-M +*> D2 = M-K ( 0 S 0 ) +*> K+L-M ( 0 0 I ) +*> P-L ( 0 0 0 ) +*> +*> N-K-L K M-K K+L-M +*> ( 0 R ) = K ( 0 R11 R12 R13 ) +*> M-K ( 0 0 R22 R23 ) +*> K+L-M ( 0 0 0 R33 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(M) ), +*> S = diag( BETA(K+1), ... , BETA(M) ), +*> C**2 + S**2 = I. +*> +*> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored +*> ( 0 R22 R23 ) +*> in B(M-K+1:L,N+M-K-L+1:N) on exit. +*> +*> The routine computes C, S, R, and optionally the unitary +*> transformation matrices U, V and Q. +*> +*> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of +*> A and B implicitly gives the SVD of A*inv(B): +*> A*inv(B) = U*(D1*inv(D2))*V**H. +*> If ( A**H,B**H)**H has orthonormal columns, then the GSVD of A and B is also +*> equal to the CS decomposition of A and B. Furthermore, the GSVD can +*> be used to derive the solution of the eigenvalue problem: +*> A**H*A x = lambda* B**H*B x. +*> In some literature, the GSVD of A and B is presented in the form +*> U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 ) +*> where U and V are orthogonal and X is nonsingular, and D1 and D2 are +*> ``diagonal''. The former GSVD form can be converted to the latter +*> form by taking the nonsingular matrix X as +*> +*> X = Q*( I 0 ) +*> ( 0 inv(R) ) +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Unitary matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Unitary matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Unitary matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose. +*> K + L = effective numerical rank of (A**H,B**H)**H. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular matrix R, or part of R. +*> See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains part of the triangular matrix R if +*> M-K-L < 0. See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is REAL array, dimension (N) +*> +*> On exit, ALPHA and BETA contain the generalized singular +*> value pairs of A and B; +*> ALPHA(1:K) = 1, +*> BETA(1:K) = 0, +*> and if M-K-L >= 0, +*> ALPHA(K+1:K+L) = C, +*> BETA(K+1:K+L) = S, +*> or if M-K-L < 0, +*> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 +*> BETA(K+1:M) =S, BETA(M+1:K+L) =1 +*> and +*> ALPHA(K+L+1:N) = 0 +*> BETA(K+L+1:N) = 0 +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX array, dimension (LDU,M) +*> If JOBU = 'U', U contains the M-by-M unitary matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX array, dimension (LDV,P) +*> If JOBV = 'V', V contains the P-by-P unitary matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is COMPLEX array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> On exit, IWORK stores the sorting information. More +*> precisely, the following loop will sort ALPHA +*> for I = K+1, min(M,K+L) +*> swap ALPHA(I) and ALPHA(IWORK(I)) +*> endfor +*> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: if INFO = 1, the Jacobi-type procedure failed to +*> converge. For further details, see subroutine CTGSJA. +*> \endverbatim +* +*> \par Internal Parameters: +* ========================= +*> +*> \verbatim +*> TOLA REAL +*> TOLB REAL +*> TOLA and TOLB are the thresholds to determine the effective +*> rank of (A**H,B**H)**H. Generally, they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup complexOTHERsing +* +*> \par Contributors: +* ================== +*> +*> Ming Gu and Huan Ren, Computer Science Division, University of +*> California at Berkeley, USA +*> +* +*> \par Further Details: +* ===================== +*> +*> CGGSVD3 replaces the deprecated subroutine CGGSVD. +*> +* ===================================================================== + SUBROUTINE CGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, + $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, LWORK, RWORK, IWORK, INFO ) +* +* -- LAPACK driver routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL ALPHA( * ), BETA( * ), RWORK( * ) + COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL WANTQ, WANTU, WANTV, LQUERY + INTEGER I, IBND, ISUB, J, NCYCLE, LWKOPT + REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL +* .. +* .. External Functions .. + LOGICAL LSAME + REAL CLANGE, SLAMCH + EXTERNAL LSAME, CLANGE, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CGGSVP3, CTGSJA, SCOPY, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Decode and test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( P.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -10 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -12 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, + $ WORK, WORK, -1, INFO ) + LWKOPT = N + INT( WORK( 1 ) ) + LWKOPT = MAX( 2*N, LWKOPT ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = CMPLX( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGGSVD3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* Compute the Frobenius norm of matrices A and B +* + ANORM = CLANGE( '1', M, N, A, LDA, RWORK ) + BNORM = CLANGE( '1', P, N, B, LDB, RWORK ) +* +* Get machine precision and set up threshold for determining +* the effective numerical rank of the matrices A and B. +* + ULP = SLAMCH( 'Precision' ) + UNFL = SLAMCH( 'Safe Minimum' ) + TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP + TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP +* + CALL CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, + $ WORK, WORK( N+1 ), LWORK-N, INFO ) +* +* Compute the GSVD of two upper "triangular" matrices +* + CALL CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, + $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, NCYCLE, INFO ) +* +* Sort the singular values and store the pivot indices in IWORK +* Copy ALPHA to RWORK, then sort ALPHA in RWORK +* + CALL SCOPY( N, ALPHA, 1, RWORK, 1 ) + IBND = MIN( L, M-K ) + DO 20 I = 1, IBND +* +* Scan for largest ALPHA(K+I) +* + ISUB = I + SMAX = RWORK( K+I ) + DO 10 J = I + 1, IBND + TEMP = RWORK( K+J ) + IF( TEMP.GT.SMAX ) THEN + ISUB = J + SMAX = TEMP + END IF + 10 CONTINUE + IF( ISUB.NE.I ) THEN + RWORK( K+ISUB ) = RWORK( K+I ) + RWORK( K+I ) = SMAX + IWORK( K+I ) = K + ISUB + ELSE + IWORK( K+I ) = K + I + END IF + 20 CONTINUE +* + WORK( 1 ) = CMPLX( LWKOPT ) + RETURN +* +* End of CGGSVD3 +* + END diff --git a/SRC/cggsvp.f b/SRC/cggsvp.f index 4a53b688..daf67eb9 100644 --- a/SRC/cggsvp.f +++ b/SRC/cggsvp.f @@ -40,6 +40,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine CGGSVP3. +*> *> CGGSVP computes unitary matrices U, V and Q such that *> *> N-K-L K L diff --git a/SRC/cggsvp3.f b/SRC/cggsvp3.f new file mode 100644 index 00000000..36fe9913 --- /dev/null +++ b/SRC/cggsvp3.f @@ -0,0 +1,580 @@ +*> \brief \b CGGSVP3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CGGSVP3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggsvp3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggsvp3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggsvp3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, +* TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, +* IWORK, RWORK, TAU, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* REAL TOLA, TOLB +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL RWORK( * ) +* COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), +* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CGGSVP3 computes unitary matrices U, V and Q such that +*> +*> N-K-L K L +*> U**H*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; +*> L ( 0 0 A23 ) +*> M-K-L ( 0 0 0 ) +*> +*> N-K-L K L +*> = K ( 0 A12 A13 ) if M-K-L < 0; +*> M-K ( 0 0 A23 ) +*> +*> N-K-L K L +*> V**H*B*Q = L ( 0 0 B13 ) +*> P-L ( 0 0 0 ) +*> +*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular +*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, +*> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective +*> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H. +*> +*> This decomposition is the preprocessing step for computing the +*> Generalized Singular Value Decomposition (GSVD), see subroutine +*> CGGSVD3. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Unitary matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Unitary matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Unitary matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular (or trapezoidal) matrix +*> described in the Purpose section. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains the triangular matrix described in +*> the Purpose section. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[in] TOLA +*> \verbatim +*> TOLA is REAL +*> \endverbatim +*> +*> \param[in] TOLB +*> \verbatim +*> TOLB is REAL +*> +*> TOLA and TOLB are the thresholds to determine the effective +*> numerical rank of matrix B and a subblock of A. Generally, +*> they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose section. +*> K + L = effective numerical rank of (A**H,B**H)**H. +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX array, dimension (LDU,M) +*> If JOBU = 'U', U contains the unitary matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX array, dimension (LDV,P) +*> If JOBV = 'V', V contains the unitary matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is COMPLEX array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the unitary matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup complexOTHERcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The subroutine uses LAPACK subroutine CGEQP3 for the QR factorization +*> with column pivoting to detect the effective numerical rank of the +*> a matrix. It may be replaced by a better rank determination strategy. +*> +*> CGGSVP3 replaces the deprecated subroutine CGGSVP. +*> +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE CGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, + $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, + $ IWORK, RWORK, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK + REAL TOLA, TOLB +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL RWORK( * ) + COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), + $ CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY + INTEGER I, J, LWKOPT + COMPLEX T +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CGEQP3, CGEQR2, CGERQ2, CLACPY, CLAPMT, + $ CLASET, CUNG2R, CUNM2R, CUNMR2, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, MAX, MIN, REAL +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + FORWRD = .TRUE. + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( P.LT.0 ) THEN + INFO = -5 + ELSE IF( N.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -8 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -10 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL CGEQP3( P, N, B, LDB, IWORK, TAU, WORK, -1, RWORK, INFO ) + LWKOPT = INT( WORK ( 1 ) ) + IF( WANTV ) THEN + LWKOPT = MAX( LWKOPT, P ) + END IF + LWKOPT = MAX( LWKOPT, MIN( N, P ) ) + LWKOPT = MAX( LWKOPT, M ) + IF( WANTQ ) THEN + LWKOPT = MAX( LWKOPT, N ) + END IF + CALL CGEQP3( M, N, A, LDA, IWORK, TAU, WORK, -1, RWORK, INFO ) + LWKOPT = MAX( LWKOPT, INT( WORK ( 1 ) ) ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = CMPLX( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CGGSVP3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* QR with column pivoting of B: B*P = V*( S11 S12 ) +* ( 0 0 ) +* + DO 10 I = 1, N + IWORK( I ) = 0 + 10 CONTINUE + CALL CGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, RWORK, INFO ) +* +* Update A := A*P +* + CALL CLAPMT( FORWRD, M, N, A, LDA, IWORK ) +* +* Determine the effective rank of matrix B. +* + L = 0 + DO 20 I = 1, MIN( P, N ) + IF( ABS( B( I, I ) ).GT.TOLB ) + $ L = L + 1 + 20 CONTINUE +* + IF( WANTV ) THEN +* +* Copy the details of V, and form V. +* + CALL CLASET( 'Full', P, P, CZERO, CZERO, V, LDV ) + IF( P.GT.1 ) + $ CALL CLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), + $ LDV ) + CALL CUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) + END IF +* +* Clean up B +* + DO 40 J = 1, L - 1 + DO 30 I = J + 1, L + B( I, J ) = CZERO + 30 CONTINUE + 40 CONTINUE + IF( P.GT.L ) + $ CALL CLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB ) +* + IF( WANTQ ) THEN +* +* Set Q = I and Update Q := Q*P +* + CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) + CALL CLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) + END IF +* + IF( P.GE.L .AND. N.NE.L ) THEN +* +* RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z +* + CALL CGERQ2( L, N, B, LDB, TAU, WORK, INFO ) +* +* Update A := A*Z**H +* + CALL CUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB, + $ TAU, A, LDA, WORK, INFO ) + IF( WANTQ ) THEN +* +* Update Q := Q*Z**H +* + CALL CUNMR2( 'Right', 'Conjugate transpose', N, N, L, B, + $ LDB, TAU, Q, LDQ, WORK, INFO ) + END IF +* +* Clean up B +* + CALL CLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB ) + DO 60 J = N - L + 1, N + DO 50 I = J - N + L + 1, L + B( I, J ) = CZERO + 50 CONTINUE + 60 CONTINUE +* + END IF +* +* Let N-L L +* A = ( A11 A12 ) M, +* +* then the following does the complete QR decomposition of A11: +* +* A11 = U*( 0 T12 )*P1**H +* ( 0 0 ) +* + DO 70 I = 1, N - L + IWORK( I ) = 0 + 70 CONTINUE + CALL CGEQP3( M, N-L, A, LDA, IWORK, TAU, WORK, LWORK, RWORK, + $ INFO ) +* +* Determine the effective rank of A11 +* + K = 0 + DO 80 I = 1, MIN( M, N-L ) + IF( ABS( A( I, I ) ).GT.TOLA ) + $ K = K + 1 + 80 CONTINUE +* +* Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N ) +* + CALL CUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ), + $ A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Copy the details of U, and form U +* + CALL CLASET( 'Full', M, M, CZERO, CZERO, U, LDU ) + IF( M.GT.1 ) + $ CALL CLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), + $ LDU ) + CALL CUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) + END IF +* + IF( WANTQ ) THEN +* +* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 +* + CALL CLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) + END IF +* +* Clean up A: set the strictly lower triangular part of +* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. +* + DO 100 J = 1, K - 1 + DO 90 I = J + 1, K + A( I, J ) = CZERO + 90 CONTINUE + 100 CONTINUE + IF( M.GT.K ) + $ CALL CLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA ) +* + IF( N-L.GT.K ) THEN +* +* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 +* + CALL CGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) +* + IF( WANTQ ) THEN +* +* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H +* + CALL CUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A, + $ LDA, TAU, Q, LDQ, WORK, INFO ) + END IF +* +* Clean up A +* + CALL CLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA ) + DO 120 J = N - L - K + 1, N - L + DO 110 I = J - N + L + K + 1, K + A( I, J ) = CZERO + 110 CONTINUE + 120 CONTINUE +* + END IF +* + IF( M.GT.K ) THEN +* +* QR factorization of A( K+1:M,N-L+1:N ) +* + CALL CGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Update U(:,K+1:M) := U(:,K+1:M)*U1 +* + CALL CUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), + $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, + $ WORK, INFO ) + END IF +* +* Clean up +* + DO 140 J = N - L + 1, N + DO 130 I = J - N + K + L + 1, M + A( I, J ) = CZERO + 130 CONTINUE + 140 CONTINUE +* + END IF +* + WORK( 1 ) = CMPLX( LWKOPT ) + RETURN +* +* End of CGGSVP3 +* + END diff --git a/SRC/dgeqp3.f b/SRC/dgeqp3.f index 62a172f7..f4cb1792 100644 --- a/SRC/dgeqp3.f +++ b/SRC/dgeqp3.f @@ -225,12 +225,6 @@ RETURN END IF * -* Quick return if possible. -* - IF( MINMN.EQ.0 ) THEN - RETURN - END IF -* * Move initial columns up front. * NFXD = 1 diff --git a/SRC/dggsvd.f b/SRC/dggsvd.f index 17dac941..6d7ace44 100644 --- a/SRC/dggsvd.f +++ b/SRC/dggsvd.f @@ -39,6 +39,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine DGGSVD3. +*> *> DGGSVD computes the generalized singular value decomposition (GSVD) *> of an M-by-N real matrix A and P-by-N real matrix B: *> diff --git a/SRC/dggsvd3.f b/SRC/dggsvd3.f new file mode 100644 index 00000000..baad1532 --- /dev/null +++ b/SRC/dggsvd3.f @@ -0,0 +1,503 @@ +*> \brief <b> DGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGGSVD3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvd3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvd3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvd3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, +* LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, +* LWORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ), +* $ BETA( * ), Q( LDQ, * ), U( LDU, * ), +* $ V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGGSVD3 computes the generalized singular value decomposition (GSVD) +*> of an M-by-N real matrix A and P-by-N real matrix B: +*> +*> U**T*A*Q = D1*( 0 R ), V**T*B*Q = D2*( 0 R ) +*> +*> where U, V and Q are orthogonal matrices. +*> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T, +*> then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and +*> D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the +*> following structures, respectively: +*> +*> If M-K-L >= 0, +*> +*> K L +*> D1 = K ( I 0 ) +*> L ( 0 C ) +*> M-K-L ( 0 0 ) +*> +*> K L +*> D2 = L ( 0 S ) +*> P-L ( 0 0 ) +*> +*> N-K-L K L +*> ( 0 R ) = K ( 0 R11 R12 ) +*> L ( 0 0 R22 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), +*> S = diag( BETA(K+1), ... , BETA(K+L) ), +*> C**2 + S**2 = I. +*> +*> R is stored in A(1:K+L,N-K-L+1:N) on exit. +*> +*> If M-K-L < 0, +*> +*> K M-K K+L-M +*> D1 = K ( I 0 0 ) +*> M-K ( 0 C 0 ) +*> +*> K M-K K+L-M +*> D2 = M-K ( 0 S 0 ) +*> K+L-M ( 0 0 I ) +*> P-L ( 0 0 0 ) +*> +*> N-K-L K M-K K+L-M +*> ( 0 R ) = K ( 0 R11 R12 R13 ) +*> M-K ( 0 0 R22 R23 ) +*> K+L-M ( 0 0 0 R33 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(M) ), +*> S = diag( BETA(K+1), ... , BETA(M) ), +*> C**2 + S**2 = I. +*> +*> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored +*> ( 0 R22 R23 ) +*> in B(M-K+1:L,N+M-K-L+1:N) on exit. +*> +*> The routine computes C, S, R, and optionally the orthogonal +*> transformation matrices U, V and Q. +*> +*> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of +*> A and B implicitly gives the SVD of A*inv(B): +*> A*inv(B) = U*(D1*inv(D2))*V**T. +*> If ( A**T,B**T)**T has orthonormal columns, then the GSVD of A and B is +*> also equal to the CS decomposition of A and B. Furthermore, the GSVD +*> can be used to derive the solution of the eigenvalue problem: +*> A**T*A x = lambda* B**T*B x. +*> In some literature, the GSVD of A and B is presented in the form +*> U**T*A*X = ( 0 D1 ), V**T*B*X = ( 0 D2 ) +*> where U and V are orthogonal and X is nonsingular, D1 and D2 are +*> ``diagonal''. The former GSVD form can be converted to the latter +*> form by taking the nonsingular matrix X as +*> +*> X = Q*( I 0 ) +*> ( 0 inv(R) ). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Orthogonal matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Orthogonal matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Orthogonal matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose. +*> K + L = effective numerical rank of (A**T,B**T)**T. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular matrix R, or part of R. +*> See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains the triangular matrix R if M-K-L < 0. +*> See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION array, dimension (N) +*> +*> On exit, ALPHA and BETA contain the generalized singular +*> value pairs of A and B; +*> ALPHA(1:K) = 1, +*> BETA(1:K) = 0, +*> and if M-K-L >= 0, +*> ALPHA(K+1:K+L) = C, +*> BETA(K+1:K+L) = S, +*> or if M-K-L < 0, +*> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 +*> BETA(K+1:M) =S, BETA(M+1:K+L) =1 +*> and +*> ALPHA(K+L+1:N) = 0 +*> BETA(K+L+1:N) = 0 +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is DOUBLE PRECISION array, dimension (LDU,M) +*> If JOBU = 'U', U contains the M-by-M orthogonal matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension (LDV,P) +*> If JOBV = 'V', V contains the P-by-P orthogonal matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is DOUBLE PRECISION array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> On exit, IWORK stores the sorting information. More +*> precisely, the following loop will sort ALPHA +*> for I = K+1, min(M,K+L) +*> swap ALPHA(I) and ALPHA(IWORK(I)) +*> endfor +*> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: if INFO = 1, the Jacobi-type procedure failed to +*> converge. For further details, see subroutine DTGSJA. +*> \endverbatim +* +*> \par Internal Parameters: +* ========================= +*> +*> \verbatim +*> TOLA DOUBLE PRECISION +*> TOLB DOUBLE PRECISION +*> TOLA and TOLB are the thresholds to determine the effective +*> rank of (A**T,B**T)**T. Generally, they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup doubleOTHERsing +* +*> \par Contributors: +* ================== +*> +*> Ming Gu and Huan Ren, Computer Science Division, University of +*> California at Berkeley, USA +*> +* +*> \par Further Details: +* ===================== +*> +*> DGGSVD3 replaces the deprecated subroutine DGGSVD. +*> +* ===================================================================== + SUBROUTINE DGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, + $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, LWORK, IWORK, INFO ) +* +* -- LAPACK driver routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ), + $ BETA( * ), Q( LDQ, * ), U( LDU, * ), + $ V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL WANTQ, WANTU, WANTV, LQUERY + INTEGER I, IBND, ISUB, J, NCYCLE, LWKOPT + DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, DLANGE + EXTERNAL LSAME, DLAMCH, DLANGE +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DGGSVP3, DTGSJA, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Decode and test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( P.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -10 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -12 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK, + $ WORK, -1, INFO ) + LWKOPT = N + INT( WORK( 1 ) ) + LWKOPT = MAX( 2*N, LWKOPT ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = DBLE( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGGSVD3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* Compute the Frobenius norm of matrices A and B +* + ANORM = DLANGE( '1', M, N, A, LDA, WORK ) + BNORM = DLANGE( '1', P, N, B, LDB, WORK ) +* +* Get machine precision and set up threshold for determining +* the effective numerical rank of the matrices A and B. +* + ULP = DLAMCH( 'Precision' ) + UNFL = DLAMCH( 'Safe Minimum' ) + TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP + TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP +* +* Preprocessing +* + CALL DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK, + $ WORK( N+1 ), LWORK-N, INFO ) +* +* Compute the GSVD of two upper "triangular" matrices +* + CALL DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, + $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, NCYCLE, INFO ) +* +* Sort the singular values and store the pivot indices in IWORK +* Copy ALPHA to WORK, then sort ALPHA in WORK +* + CALL DCOPY( N, ALPHA, 1, WORK, 1 ) + IBND = MIN( L, M-K ) + DO 20 I = 1, IBND +* +* Scan for largest ALPHA(K+I) +* + ISUB = I + SMAX = WORK( K+I ) + DO 10 J = I + 1, IBND + TEMP = WORK( K+J ) + IF( TEMP.GT.SMAX ) THEN + ISUB = J + SMAX = TEMP + END IF + 10 CONTINUE + IF( ISUB.NE.I ) THEN + WORK( K+ISUB ) = WORK( K+I ) + WORK( K+I ) = SMAX + IWORK( K+I ) = K + ISUB + ELSE + IWORK( K+I ) = K + I + END IF + 20 CONTINUE +* + WORK( 1 ) = DBLE( LWKOPT ) + RETURN +* +* End of DGGSVD3 +* + END diff --git a/SRC/dggsvp.f b/SRC/dggsvp.f index 567133d6..7e195b0b 100644 --- a/SRC/dggsvp.f +++ b/SRC/dggsvp.f @@ -39,6 +39,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine DGGSVP3. +*> *> DGGSVP computes orthogonal matrices U, V and Q such that *> *> N-K-L K L diff --git a/SRC/dggsvp3.f b/SRC/dggsvp3.f new file mode 100644 index 00000000..6041e4ea --- /dev/null +++ b/SRC/dggsvp3.f @@ -0,0 +1,571 @@ +*> \brief \b DGGSVP3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DGGSVP3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvp3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvp3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvp3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, +* TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, +* IWORK, TAU, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* DOUBLE PRECISION TOLA, TOLB +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), +* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGGSVP3 computes orthogonal matrices U, V and Q such that +*> +*> N-K-L K L +*> U**T*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; +*> L ( 0 0 A23 ) +*> M-K-L ( 0 0 0 ) +*> +*> N-K-L K L +*> = K ( 0 A12 A13 ) if M-K-L < 0; +*> M-K ( 0 0 A23 ) +*> +*> N-K-L K L +*> V**T*B*Q = L ( 0 0 B13 ) +*> P-L ( 0 0 0 ) +*> +*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular +*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, +*> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective +*> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T. +*> +*> This decomposition is the preprocessing step for computing the +*> Generalized Singular Value Decomposition (GSVD), see subroutine +*> DGGSVD3. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Orthogonal matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Orthogonal matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Orthogonal matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular (or trapezoidal) matrix +*> described in the Purpose section. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains the triangular matrix described in +*> the Purpose section. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[in] TOLA +*> \verbatim +*> TOLA is DOUBLE PRECISION +*> \endverbatim +*> +*> \param[in] TOLB +*> \verbatim +*> TOLB is DOUBLE PRECISION +*> +*> TOLA and TOLB are the thresholds to determine the effective +*> numerical rank of matrix B and a subblock of A. Generally, +*> they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose section. +*> K + L = effective numerical rank of (A**T,B**T)**T. +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is DOUBLE PRECISION array, dimension (LDU,M) +*> If JOBU = 'U', U contains the orthogonal matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension (LDV,P) +*> If JOBV = 'V', V contains the orthogonal matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is DOUBLE PRECISION array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the orthogonal matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup doubleOTHERcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The subroutine uses LAPACK subroutine DGEQP3 for the QR factorization +*> with column pivoting to detect the effective numerical rank of the +*> a matrix. It may be replaced by a better rank determination strategy. +*> +*> DGGSVP3 replaces the deprecated subroutine DGGSVP. +*> +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, + $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, + $ IWORK, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK + DOUBLE PRECISION TOLA, TOLB +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY + INTEGER I, J, LWKOPT +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DGEQP3, DGEQR2, DGERQ2, DLACPY, DLAPMT, + $ DLASET, DORG2R, DORM2R, DORMR2, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + FORWRD = .TRUE. + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( P.LT.0 ) THEN + INFO = -5 + ELSE IF( N.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -8 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -10 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL DGEQP3( P, N, B, LDB, IWORK, TAU, WORK, -1, INFO ) + LWKOPT = INT( WORK ( 1 ) ) + IF( WANTV ) THEN + LWKOPT = MAX( LWKOPT, P ) + END IF + LWKOPT = MAX( LWKOPT, MIN( N, P ) ) + LWKOPT = MAX( LWKOPT, M ) + IF( WANTQ ) THEN + LWKOPT = MAX( LWKOPT, N ) + END IF + CALL DGEQP3( M, N, A, LDA, IWORK, TAU, WORK, -1, INFO ) + LWKOPT = MAX( LWKOPT, INT( WORK ( 1 ) ) ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = DBLE( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DGGSVP3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* QR with column pivoting of B: B*P = V*( S11 S12 ) +* ( 0 0 ) +* + DO 10 I = 1, N + IWORK( I ) = 0 + 10 CONTINUE + CALL DGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, INFO ) +* +* Update A := A*P +* + CALL DLAPMT( FORWRD, M, N, A, LDA, IWORK ) +* +* Determine the effective rank of matrix B. +* + L = 0 + DO 20 I = 1, MIN( P, N ) + IF( ABS( B( I, I ) ).GT.TOLB ) + $ L = L + 1 + 20 CONTINUE +* + IF( WANTV ) THEN +* +* Copy the details of V, and form V. +* + CALL DLASET( 'Full', P, P, ZERO, ZERO, V, LDV ) + IF( P.GT.1 ) + $ CALL DLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), + $ LDV ) + CALL DORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) + END IF +* +* Clean up B +* + DO 40 J = 1, L - 1 + DO 30 I = J + 1, L + B( I, J ) = ZERO + 30 CONTINUE + 40 CONTINUE + IF( P.GT.L ) + $ CALL DLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB ) +* + IF( WANTQ ) THEN +* +* Set Q = I and Update Q := Q*P +* + CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) + CALL DLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) + END IF +* + IF( P.GE.L .AND. N.NE.L ) THEN +* +* RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z +* + CALL DGERQ2( L, N, B, LDB, TAU, WORK, INFO ) +* +* Update A := A*Z**T +* + CALL DORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A, + $ LDA, WORK, INFO ) +* + IF( WANTQ ) THEN +* +* Update Q := Q*Z**T +* + CALL DORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q, + $ LDQ, WORK, INFO ) + END IF +* +* Clean up B +* + CALL DLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB ) + DO 60 J = N - L + 1, N + DO 50 I = J - N + L + 1, L + B( I, J ) = ZERO + 50 CONTINUE + 60 CONTINUE +* + END IF +* +* Let N-L L +* A = ( A11 A12 ) M, +* +* then the following does the complete QR decomposition of A11: +* +* A11 = U*( 0 T12 )*P1**T +* ( 0 0 ) +* + DO 70 I = 1, N - L + IWORK( I ) = 0 + 70 CONTINUE + CALL DGEQP3( M, N-L, A, LDA, IWORK, TAU, WORK, LWORK, INFO ) +* +* Determine the effective rank of A11 +* + K = 0 + DO 80 I = 1, MIN( M, N-L ) + IF( ABS( A( I, I ) ).GT.TOLA ) + $ K = K + 1 + 80 CONTINUE +* +* Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N ) +* + CALL DORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA, + $ TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Copy the details of U, and form U +* + CALL DLASET( 'Full', M, M, ZERO, ZERO, U, LDU ) + IF( M.GT.1 ) + $ CALL DLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), + $ LDU ) + CALL DORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) + END IF +* + IF( WANTQ ) THEN +* +* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 +* + CALL DLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) + END IF +* +* Clean up A: set the strictly lower triangular part of +* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. +* + DO 100 J = 1, K - 1 + DO 90 I = J + 1, K + A( I, J ) = ZERO + 90 CONTINUE + 100 CONTINUE + IF( M.GT.K ) + $ CALL DLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA ) +* + IF( N-L.GT.K ) THEN +* +* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 +* + CALL DGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) +* + IF( WANTQ ) THEN +* +* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T +* + CALL DORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU, + $ Q, LDQ, WORK, INFO ) + END IF +* +* Clean up A +* + CALL DLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA ) + DO 120 J = N - L - K + 1, N - L + DO 110 I = J - N + L + K + 1, K + A( I, J ) = ZERO + 110 CONTINUE + 120 CONTINUE +* + END IF +* + IF( M.GT.K ) THEN +* +* QR factorization of A( K+1:M,N-L+1:N ) +* + CALL DGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Update U(:,K+1:M) := U(:,K+1:M)*U1 +* + CALL DORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), + $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, + $ WORK, INFO ) + END IF +* +* Clean up +* + DO 140 J = N - L + 1, N + DO 130 I = J - N + K + L + 1, M + A( I, J ) = ZERO + 130 CONTINUE + 140 CONTINUE +* + END IF +* + WORK( 1 ) = DBLE( LWKOPT ) + RETURN +* +* End of DGGSVP3 +* + END diff --git a/SRC/sgeqp3.f b/SRC/sgeqp3.f index 23f3af25..abc61ac5 100644 --- a/SRC/sgeqp3.f +++ b/SRC/sgeqp3.f @@ -185,8 +185,8 @@ * .. * .. Intrinsic Functions .. INTRINSIC INT, MAX, MIN -* .. -* .. Executable Statements .. +* Test input arguments +* ==================== * INFO = 0 LQUERY = ( LWORK.EQ.-1 ) @@ -222,12 +222,6 @@ RETURN END IF * -* Quick return if possible. -* - IF( MINMN.EQ.0 ) THEN - RETURN - END IF -* * Move initial columns up front. * NFXD = 1 diff --git a/SRC/sggsvd.f b/SRC/sggsvd.f index 63201976..0bf38806 100644 --- a/SRC/sggsvd.f +++ b/SRC/sggsvd.f @@ -39,6 +39,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine SGGSVD3. +*> *> SGGSVD computes the generalized singular value decomposition (GSVD) *> of an M-by-N real matrix A and P-by-N real matrix B: *> diff --git a/SRC/sggsvd3.f b/SRC/sggsvd3.f new file mode 100644 index 00000000..3e7b114e --- /dev/null +++ b/SRC/sggsvd3.f @@ -0,0 +1,503 @@ +*> \brief <b> SGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SGGSVD3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggsvd3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggsvd3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggsvd3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, +* LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, +* LWORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL A( LDA, * ), ALPHA( * ), B( LDB, * ), +* $ BETA( * ), Q( LDQ, * ), U( LDU, * ), +* $ V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SGGSVD3 computes the generalized singular value decomposition (GSVD) +*> of an M-by-N real matrix A and P-by-N real matrix B: +*> +*> U**T*A*Q = D1*( 0 R ), V**T*B*Q = D2*( 0 R ) +*> +*> where U, V and Q are orthogonal matrices. +*> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T, +*> then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and +*> D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the +*> following structures, respectively: +*> +*> If M-K-L >= 0, +*> +*> K L +*> D1 = K ( I 0 ) +*> L ( 0 C ) +*> M-K-L ( 0 0 ) +*> +*> K L +*> D2 = L ( 0 S ) +*> P-L ( 0 0 ) +*> +*> N-K-L K L +*> ( 0 R ) = K ( 0 R11 R12 ) +*> L ( 0 0 R22 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), +*> S = diag( BETA(K+1), ... , BETA(K+L) ), +*> C**2 + S**2 = I. +*> +*> R is stored in A(1:K+L,N-K-L+1:N) on exit. +*> +*> If M-K-L < 0, +*> +*> K M-K K+L-M +*> D1 = K ( I 0 0 ) +*> M-K ( 0 C 0 ) +*> +*> K M-K K+L-M +*> D2 = M-K ( 0 S 0 ) +*> K+L-M ( 0 0 I ) +*> P-L ( 0 0 0 ) +*> +*> N-K-L K M-K K+L-M +*> ( 0 R ) = K ( 0 R11 R12 R13 ) +*> M-K ( 0 0 R22 R23 ) +*> K+L-M ( 0 0 0 R33 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(M) ), +*> S = diag( BETA(K+1), ... , BETA(M) ), +*> C**2 + S**2 = I. +*> +*> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored +*> ( 0 R22 R23 ) +*> in B(M-K+1:L,N+M-K-L+1:N) on exit. +*> +*> The routine computes C, S, R, and optionally the orthogonal +*> transformation matrices U, V and Q. +*> +*> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of +*> A and B implicitly gives the SVD of A*inv(B): +*> A*inv(B) = U*(D1*inv(D2))*V**T. +*> If ( A**T,B**T)**T has orthonormal columns, then the GSVD of A and B is +*> also equal to the CS decomposition of A and B. Furthermore, the GSVD +*> can be used to derive the solution of the eigenvalue problem: +*> A**T*A x = lambda* B**T*B x. +*> In some literature, the GSVD of A and B is presented in the form +*> U**T*A*X = ( 0 D1 ), V**T*B*X = ( 0 D2 ) +*> where U and V are orthogonal and X is nonsingular, D1 and D2 are +*> ``diagonal''. The former GSVD form can be converted to the latter +*> form by taking the nonsingular matrix X as +*> +*> X = Q*( I 0 ) +*> ( 0 inv(R) ). +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Orthogonal matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Orthogonal matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Orthogonal matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose. +*> K + L = effective numerical rank of (A**T,B**T)**T. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular matrix R, or part of R. +*> See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is REAL array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains the triangular matrix R if M-K-L < 0. +*> See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is REAL array, dimension (N) +*> +*> On exit, ALPHA and BETA contain the generalized singular +*> value pairs of A and B; +*> ALPHA(1:K) = 1, +*> BETA(1:K) = 0, +*> and if M-K-L >= 0, +*> ALPHA(K+1:K+L) = C, +*> BETA(K+1:K+L) = S, +*> or if M-K-L < 0, +*> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 +*> BETA(K+1:M) =S, BETA(M+1:K+L) =1 +*> and +*> ALPHA(K+L+1:N) = 0 +*> BETA(K+L+1:N) = 0 +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is REAL array, dimension (LDU,M) +*> If JOBU = 'U', U contains the M-by-M orthogonal matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is REAL array, dimension (LDV,P) +*> If JOBV = 'V', V contains the P-by-P orthogonal matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is REAL array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> On exit, IWORK stores the sorting information. More +*> precisely, the following loop will sort ALPHA +*> for I = K+1, min(M,K+L) +*> swap ALPHA(I) and ALPHA(IWORK(I)) +*> endfor +*> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: if INFO = 1, the Jacobi-type procedure failed to +*> converge. For further details, see subroutine STGSJA. +*> \endverbatim +* +*> \par Internal Parameters: +* ========================= +*> +*> \verbatim +*> TOLA REAL +*> TOLB REAL +*> TOLA and TOLB are the thresholds to determine the effective +*> rank of (A**T,B**T)**T. Generally, they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup realOTHERsing +* +*> \par Contributors: +* ================== +*> +*> Ming Gu and Huan Ren, Computer Science Division, University of +*> California at Berkeley, USA +*> +* +*> \par Further Details: +* ===================== +*> +*> SGGSVD3 replaces the deprecated subroutine SGGSVD. +*> +* ===================================================================== + SUBROUTINE SGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, + $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, LWORK, IWORK, INFO ) +* +* -- LAPACK driver routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL A( LDA, * ), ALPHA( * ), B( LDB, * ), + $ BETA( * ), Q( LDQ, * ), U( LDU, * ), + $ V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL WANTQ, WANTU, WANTV, LQUERY + INTEGER I, IBND, ISUB, J, NCYCLE, LWKOPT + REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL +* .. +* .. External Functions .. + LOGICAL LSAME + REAL SLAMCH, SLANGE + EXTERNAL LSAME, SLAMCH, SLANGE +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SGGSVP3, STGSJA, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Decode and test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( P.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -10 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -12 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL SGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK, + $ WORK, -1, INFO ) + LWKOPT = N + INT( WORK( 1 ) ) + LWKOPT = MAX( 2*N, LWKOPT ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = REAL( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SGGSVD3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* Compute the Frobenius norm of matrices A and B +* + ANORM = SLANGE( '1', M, N, A, LDA, WORK ) + BNORM = SLANGE( '1', P, N, B, LDB, WORK ) +* +* Get machine precision and set up threshold for determining +* the effective numerical rank of the matrices A and B. +* + ULP = SLAMCH( 'Precision' ) + UNFL = SLAMCH( 'Safe Minimum' ) + TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP + TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP +* +* Preprocessing +* + CALL SGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK, + $ WORK( N+1 ), LWORK-N, INFO ) +* +* Compute the GSVD of two upper "triangular" matrices +* + CALL STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, + $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, NCYCLE, INFO ) +* +* Sort the singular values and store the pivot indices in IWORK +* Copy ALPHA to WORK, then sort ALPHA in WORK +* + CALL SCOPY( N, ALPHA, 1, WORK, 1 ) + IBND = MIN( L, M-K ) + DO 20 I = 1, IBND +* +* Scan for largest ALPHA(K+I) +* + ISUB = I + SMAX = WORK( K+I ) + DO 10 J = I + 1, IBND + TEMP = WORK( K+J ) + IF( TEMP.GT.SMAX ) THEN + ISUB = J + SMAX = TEMP + END IF + 10 CONTINUE + IF( ISUB.NE.I ) THEN + WORK( K+ISUB ) = WORK( K+I ) + WORK( K+I ) = SMAX + IWORK( K+I ) = K + ISUB + ELSE + IWORK( K+I ) = K + I + END IF + 20 CONTINUE +* + WORK( 1 ) = REAL( LWKOPT ) + RETURN +* +* End of SGGSVD3 +* + END diff --git a/SRC/sggsvp.f b/SRC/sggsvp.f index f89c352b..0bbb30b9 100644 --- a/SRC/sggsvp.f +++ b/SRC/sggsvp.f @@ -39,6 +39,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine SGGSVP3. +*> *> SGGSVP computes orthogonal matrices U, V and Q such that *> *> N-K-L K L diff --git a/SRC/sggsvp3.f b/SRC/sggsvp3.f new file mode 100644 index 00000000..f5496246 --- /dev/null +++ b/SRC/sggsvp3.f @@ -0,0 +1,571 @@ +*> \brief \b SGGSVP3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SGGSVP3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggsvp3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggsvp3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggsvp3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, +* TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, +* IWORK, TAU, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* REAL TOLA, TOLB +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), +* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SGGSVP3 computes orthogonal matrices U, V and Q such that +*> +*> N-K-L K L +*> U**T*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; +*> L ( 0 0 A23 ) +*> M-K-L ( 0 0 0 ) +*> +*> N-K-L K L +*> = K ( 0 A12 A13 ) if M-K-L < 0; +*> M-K ( 0 0 A23 ) +*> +*> N-K-L K L +*> V**T*B*Q = L ( 0 0 B13 ) +*> P-L ( 0 0 0 ) +*> +*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular +*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, +*> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective +*> numerical rank of the (M+P)-by-N matrix (A**T,B**T)**T. +*> +*> This decomposition is the preprocessing step for computing the +*> Generalized Singular Value Decomposition (GSVD), see subroutine +*> SGGSVD3. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Orthogonal matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Orthogonal matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Orthogonal matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular (or trapezoidal) matrix +*> described in the Purpose section. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is REAL array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains the triangular matrix described in +*> the Purpose section. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[in] TOLA +*> \verbatim +*> TOLA is REAL +*> \endverbatim +*> +*> \param[in] TOLB +*> \verbatim +*> TOLB is REAL +*> +*> TOLA and TOLB are the thresholds to determine the effective +*> numerical rank of matrix B and a subblock of A. Generally, +*> they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose section. +*> K + L = effective numerical rank of (A**T,B**T)**T. +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is REAL array, dimension (LDU,M) +*> If JOBU = 'U', U contains the orthogonal matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is REAL array, dimension (LDV,P) +*> If JOBV = 'V', V contains the orthogonal matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is REAL array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the orthogonal matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup realOTHERcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The subroutine uses LAPACK subroutine SGEQP3 for the QR factorization +*> with column pivoting to detect the effective numerical rank of the +*> a matrix. It may be replaced by a better rank determination strategy. +*> +*> SGGSVP3 replaces the deprecated subroutine SGGSVP. +*> +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE SGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, + $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, + $ IWORK, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK + REAL TOLA, TOLB +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY + INTEGER I, J, LWKOPT +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL SGEQP3, SGEQR2, SGERQ2, SLACPY, SLAPMT, + $ SLASET, SORG2R, SORM2R, SORMR2, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + FORWRD = .TRUE. + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( P.LT.0 ) THEN + INFO = -5 + ELSE IF( N.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -8 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -10 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL SGEQP3( P, N, B, LDB, IWORK, TAU, WORK, -1, INFO ) + LWKOPT = INT( WORK ( 1 ) ) + IF( WANTV ) THEN + LWKOPT = MAX( LWKOPT, P ) + END IF + LWKOPT = MAX( LWKOPT, MIN( N, P ) ) + LWKOPT = MAX( LWKOPT, M ) + IF( WANTQ ) THEN + LWKOPT = MAX( LWKOPT, N ) + END IF + CALL SGEQP3( M, N, A, LDA, IWORK, TAU, WORK, -1, INFO ) + LWKOPT = MAX( LWKOPT, INT( WORK ( 1 ) ) ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = REAL( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SGGSVP3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* QR with column pivoting of B: B*P = V*( S11 S12 ) +* ( 0 0 ) +* + DO 10 I = 1, N + IWORK( I ) = 0 + 10 CONTINUE + CALL SGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, INFO ) +* +* Update A := A*P +* + CALL SLAPMT( FORWRD, M, N, A, LDA, IWORK ) +* +* Determine the effective rank of matrix B. +* + L = 0 + DO 20 I = 1, MIN( P, N ) + IF( ABS( B( I, I ) ).GT.TOLB ) + $ L = L + 1 + 20 CONTINUE +* + IF( WANTV ) THEN +* +* Copy the details of V, and form V. +* + CALL SLASET( 'Full', P, P, ZERO, ZERO, V, LDV ) + IF( P.GT.1 ) + $ CALL SLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), + $ LDV ) + CALL SORG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) + END IF +* +* Clean up B +* + DO 40 J = 1, L - 1 + DO 30 I = J + 1, L + B( I, J ) = ZERO + 30 CONTINUE + 40 CONTINUE + IF( P.GT.L ) + $ CALL SLASET( 'Full', P-L, N, ZERO, ZERO, B( L+1, 1 ), LDB ) +* + IF( WANTQ ) THEN +* +* Set Q = I and Update Q := Q*P +* + CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) + CALL SLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) + END IF +* + IF( P.GE.L .AND. N.NE.L ) THEN +* +* RQ factorization of (S11 S12): ( S11 S12 ) = ( 0 S12 )*Z +* + CALL SGERQ2( L, N, B, LDB, TAU, WORK, INFO ) +* +* Update A := A*Z**T +* + CALL SORMR2( 'Right', 'Transpose', M, N, L, B, LDB, TAU, A, + $ LDA, WORK, INFO ) +* + IF( WANTQ ) THEN +* +* Update Q := Q*Z**T +* + CALL SORMR2( 'Right', 'Transpose', N, N, L, B, LDB, TAU, Q, + $ LDQ, WORK, INFO ) + END IF +* +* Clean up B +* + CALL SLASET( 'Full', L, N-L, ZERO, ZERO, B, LDB ) + DO 60 J = N - L + 1, N + DO 50 I = J - N + L + 1, L + B( I, J ) = ZERO + 50 CONTINUE + 60 CONTINUE +* + END IF +* +* Let N-L L +* A = ( A11 A12 ) M, +* +* then the following does the complete QR decomposition of A11: +* +* A11 = U*( 0 T12 )*P1**T +* ( 0 0 ) +* + DO 70 I = 1, N - L + IWORK( I ) = 0 + 70 CONTINUE + CALL SGEQP3( M, N-L, A, LDA, IWORK, TAU, WORK, LWORK, INFO ) +* +* Determine the effective rank of A11 +* + K = 0 + DO 80 I = 1, MIN( M, N-L ) + IF( ABS( A( I, I ) ).GT.TOLA ) + $ K = K + 1 + 80 CONTINUE +* +* Update A12 := U**T*A12, where A12 = A( 1:M, N-L+1:N ) +* + CALL SORM2R( 'Left', 'Transpose', M, L, MIN( M, N-L ), A, LDA, + $ TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Copy the details of U, and form U +* + CALL SLASET( 'Full', M, M, ZERO, ZERO, U, LDU ) + IF( M.GT.1 ) + $ CALL SLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), + $ LDU ) + CALL SORG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) + END IF +* + IF( WANTQ ) THEN +* +* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 +* + CALL SLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) + END IF +* +* Clean up A: set the strictly lower triangular part of +* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. +* + DO 100 J = 1, K - 1 + DO 90 I = J + 1, K + A( I, J ) = ZERO + 90 CONTINUE + 100 CONTINUE + IF( M.GT.K ) + $ CALL SLASET( 'Full', M-K, N-L, ZERO, ZERO, A( K+1, 1 ), LDA ) +* + IF( N-L.GT.K ) THEN +* +* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 +* + CALL SGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) +* + IF( WANTQ ) THEN +* +* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**T +* + CALL SORMR2( 'Right', 'Transpose', N, N-L, K, A, LDA, TAU, + $ Q, LDQ, WORK, INFO ) + END IF +* +* Clean up A +* + CALL SLASET( 'Full', K, N-L-K, ZERO, ZERO, A, LDA ) + DO 120 J = N - L - K + 1, N - L + DO 110 I = J - N + L + K + 1, K + A( I, J ) = ZERO + 110 CONTINUE + 120 CONTINUE +* + END IF +* + IF( M.GT.K ) THEN +* +* QR factorization of A( K+1:M,N-L+1:N ) +* + CALL SGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Update U(:,K+1:M) := U(:,K+1:M)*U1 +* + CALL SORM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), + $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, + $ WORK, INFO ) + END IF +* +* Clean up +* + DO 140 J = N - L + 1, N + DO 130 I = J - N + K + L + 1, M + A( I, J ) = ZERO + 130 CONTINUE + 140 CONTINUE +* + END IF +* + WORK( 1 ) = REAL( LWKOPT ) + RETURN +* +* End of SGGSVP3 +* + END diff --git a/SRC/zgeqp3.f b/SRC/zgeqp3.f index 8b4b620a..fc8d5933 100644 --- a/SRC/zgeqp3.f +++ b/SRC/zgeqp3.f @@ -220,7 +220,7 @@ NB = ILAENV( INB, 'ZGEQRF', ' ', M, N, -1, -1 ) LWKOPT = ( N + 1 )*NB END IF - WORK( 1 ) = LWKOPT + WORK( 1 ) = DCMPLX( LWKOPT ) * IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN INFO = -8 @@ -234,12 +234,6 @@ RETURN END IF * -* Quick return if possible. -* - IF( MINMN.EQ.0 ) THEN - RETURN - END IF -* * Move initial columns up front. * NFXD = 1 @@ -370,7 +364,7 @@ * END IF * - WORK( 1 ) = IWS + WORK( 1 ) = DCMPLX( LWKOPT ) RETURN * * End of ZGEQP3 diff --git a/SRC/zggsvd.f b/SRC/zggsvd.f index b27c7fbd..db82910e 100644 --- a/SRC/zggsvd.f +++ b/SRC/zggsvd.f @@ -39,6 +39,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine ZGGSVD3. +*> *> ZGGSVD computes the generalized singular value decomposition (GSVD) *> of an M-by-N complex matrix A and P-by-N complex matrix B: *> diff --git a/SRC/zggsvd3.f b/SRC/zggsvd3.f new file mode 100644 index 00000000..d478d292 --- /dev/null +++ b/SRC/zggsvd3.f @@ -0,0 +1,505 @@ +*> \brief <b> ZGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices</b> +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGGSVD3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggsvd3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggsvd3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggsvd3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, +* LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, +* LWORK, RWORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), +* $ U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZGGSVD3 computes the generalized singular value decomposition (GSVD) +*> of an M-by-N complex matrix A and P-by-N complex matrix B: +*> +*> U**H*A*Q = D1*( 0 R ), V**H*B*Q = D2*( 0 R ) +*> +*> where U, V and Q are unitary matrices. +*> Let K+L = the effective numerical rank of the +*> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper +*> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) "diagonal" +*> matrices and of the following structures, respectively: +*> +*> If M-K-L >= 0, +*> +*> K L +*> D1 = K ( I 0 ) +*> L ( 0 C ) +*> M-K-L ( 0 0 ) +*> +*> K L +*> D2 = L ( 0 S ) +*> P-L ( 0 0 ) +*> +*> N-K-L K L +*> ( 0 R ) = K ( 0 R11 R12 ) +*> L ( 0 0 R22 ) +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), +*> S = diag( BETA(K+1), ... , BETA(K+L) ), +*> C**2 + S**2 = I. +*> +*> R is stored in A(1:K+L,N-K-L+1:N) on exit. +*> +*> If M-K-L < 0, +*> +*> K M-K K+L-M +*> D1 = K ( I 0 0 ) +*> M-K ( 0 C 0 ) +*> +*> K M-K K+L-M +*> D2 = M-K ( 0 S 0 ) +*> K+L-M ( 0 0 I ) +*> P-L ( 0 0 0 ) +*> +*> N-K-L K M-K K+L-M +*> ( 0 R ) = K ( 0 R11 R12 R13 ) +*> M-K ( 0 0 R22 R23 ) +*> K+L-M ( 0 0 0 R33 ) +*> +*> where +*> +*> C = diag( ALPHA(K+1), ... , ALPHA(M) ), +*> S = diag( BETA(K+1), ... , BETA(M) ), +*> C**2 + S**2 = I. +*> +*> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored +*> ( 0 R22 R23 ) +*> in B(M-K+1:L,N+M-K-L+1:N) on exit. +*> +*> The routine computes C, S, R, and optionally the unitary +*> transformation matrices U, V and Q. +*> +*> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of +*> A and B implicitly gives the SVD of A*inv(B): +*> A*inv(B) = U*(D1*inv(D2))*V**H. +*> If ( A**H,B**H)**H has orthonormal columns, then the GSVD of A and B is also +*> equal to the CS decomposition of A and B. Furthermore, the GSVD can +*> be used to derive the solution of the eigenvalue problem: +*> A**H*A x = lambda* B**H*B x. +*> In some literature, the GSVD of A and B is presented in the form +*> U**H*A*X = ( 0 D1 ), V**H*B*X = ( 0 D2 ) +*> where U and V are orthogonal and X is nonsingular, and D1 and D2 are +*> ``diagonal''. The former GSVD form can be converted to the latter +*> form by taking the nonsingular matrix X as +*> +*> X = Q*( I 0 ) +*> ( 0 inv(R) ) +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Unitary matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Unitary matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Unitary matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose. +*> K + L = effective numerical rank of (A**H,B**H)**H. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular matrix R, or part of R. +*> See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains part of the triangular matrix R if +*> M-K-L < 0. See Purpose for details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION array, dimension (N) +*> +*> On exit, ALPHA and BETA contain the generalized singular +*> value pairs of A and B; +*> ALPHA(1:K) = 1, +*> BETA(1:K) = 0, +*> and if M-K-L >= 0, +*> ALPHA(K+1:K+L) = C, +*> BETA(K+1:K+L) = S, +*> or if M-K-L < 0, +*> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 +*> BETA(K+1:M) =S, BETA(M+1:K+L) =1 +*> and +*> ALPHA(K+L+1:N) = 0 +*> BETA(K+L+1:N) = 0 +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX*16 array, dimension (LDU,M) +*> If JOBU = 'U', U contains the M-by-M unitary matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX*16 array, dimension (LDV,P) +*> If JOBV = 'V', V contains the P-by-P unitary matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is COMPLEX*16 array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> On exit, IWORK stores the sorting information. More +*> precisely, the following loop will sort ALPHA +*> for I = K+1, min(M,K+L) +*> swap ALPHA(I) and ALPHA(IWORK(I)) +*> endfor +*> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: if INFO = 1, the Jacobi-type procedure failed to +*> converge. For further details, see subroutine ZTGSJA. +*> \endverbatim +* +*> \par Internal Parameters: +* ========================= +*> +*> \verbatim +*> TOLA DOUBLE PRECISION +*> TOLB DOUBLE PRECISION +*> TOLA and TOLB are the thresholds to determine the effective +*> rank of (A**H,B**H)**H. Generally, they are set to +*> TOLA = MAX(M,N)*norm(A)*MACHEPS, +*> TOLB = MAX(P,N)*norm(B)*MACHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup complex16OTHERsing +* +*> \par Contributors: +* ================== +*> +*> Ming Gu and Huan Ren, Computer Science Division, University of +*> California at Berkeley, USA +*> +* +*> \par Further Details: +* ===================== +*> +*> ZGGSVD3 replaces the deprecated subroutine ZGGSVD. +*> +* ===================================================================== + SUBROUTINE ZGGSVD3( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, + $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, LWORK, RWORK, IWORK, INFO ) +* +* -- LAPACK driver routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * ) + COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Local Scalars .. + LOGICAL WANTQ, WANTU, WANTV, LQUERY + INTEGER I, IBND, ISUB, J, NCYCLE, LWKOPT + DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL +* .. +* .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH, ZLANGE + EXTERNAL LSAME, DLAMCH, ZLANGE +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, XERBLA, ZGGSVP, ZTGSJA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. Executable Statements .. +* +* Decode and test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( P.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -10 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -12 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL ZGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, + $ WORK, WORK, -1, INFO ) + LWKOPT = N + INT( WORK( 1 ) ) + LWKOPT = MAX( 2*N, LWKOPT ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = DCMPLX( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGGSVD3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* Compute the Frobenius norm of matrices A and B +* + ANORM = ZLANGE( '1', M, N, A, LDA, RWORK ) + BNORM = ZLANGE( '1', P, N, B, LDB, RWORK ) +* +* Get machine precision and set up threshold for determining +* the effective numerical rank of the matrices A and B. +* + ULP = DLAMCH( 'Precision' ) + UNFL = DLAMCH( 'Safe Minimum' ) + TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP + TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP +* + CALL ZGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, + $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, RWORK, + $ WORK, WORK( N+1 ), LWORK-N, INFO ) +* +* Compute the GSVD of two upper "triangular" matrices +* + CALL ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, + $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, + $ WORK, NCYCLE, INFO ) +* +* Sort the singular values and store the pivot indices in IWORK +* Copy ALPHA to RWORK, then sort ALPHA in RWORK +* + CALL DCOPY( N, ALPHA, 1, RWORK, 1 ) + IBND = MIN( L, M-K ) + DO 20 I = 1, IBND +* +* Scan for largest ALPHA(K+I) +* + ISUB = I + SMAX = RWORK( K+I ) + DO 10 J = I + 1, IBND + TEMP = RWORK( K+J ) + IF( TEMP.GT.SMAX ) THEN + ISUB = J + SMAX = TEMP + END IF + 10 CONTINUE + IF( ISUB.NE.I ) THEN + RWORK( K+ISUB ) = RWORK( K+I ) + RWORK( K+I ) = SMAX + IWORK( K+I ) = K + ISUB + ELSE + IWORK( K+I ) = K + I + END IF + 20 CONTINUE +* + WORK( 1 ) = DCMPLX( LWKOPT ) + RETURN +* +* End of ZGGSVD3 +* + END diff --git a/SRC/zggsvp.f b/SRC/zggsvp.f index 40d8c5b0..aff6c665 100644 --- a/SRC/zggsvp.f +++ b/SRC/zggsvp.f @@ -40,6 +40,8 @@ *> *> \verbatim *> +*> This routine is deprecated and has been replaced by routine ZGGSVP3. +*> *> ZGGSVP computes unitary matrices U, V and Q such that *> *> N-K-L K L diff --git a/SRC/zggsvp3.f b/SRC/zggsvp3.f new file mode 100644 index 00000000..b397651c --- /dev/null +++ b/SRC/zggsvp3.f @@ -0,0 +1,580 @@ +*> \brief \b ZGGSVP3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZGGSVP3 + dependencies +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggsvp3.f"> +*> [TGZ]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggsvp3.f"> +*> [ZIP]</a> +*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggsvp3.f"> +*> [TXT]</a> +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, +* TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, +* IWORK, RWORK, TAU, WORK, LWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBQ, JOBU, JOBV +* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, LWORK +* DOUBLE PRECISION TOLA, TOLB +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION RWORK( * ) +* COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), +* $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZGGSVP3 computes unitary matrices U, V and Q such that +*> +*> N-K-L K L +*> U**H*A*Q = K ( 0 A12 A13 ) if M-K-L >= 0; +*> L ( 0 0 A23 ) +*> M-K-L ( 0 0 0 ) +*> +*> N-K-L K L +*> = K ( 0 A12 A13 ) if M-K-L < 0; +*> M-K ( 0 0 A23 ) +*> +*> N-K-L K L +*> V**H*B*Q = L ( 0 0 B13 ) +*> P-L ( 0 0 0 ) +*> +*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular +*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, +*> otherwise A23 is (M-K)-by-L upper trapezoidal. K+L = the effective +*> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H. +*> +*> This decomposition is the preprocessing step for computing the +*> Generalized Singular Value Decomposition (GSVD), see subroutine +*> ZGGSVD3. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] JOBU +*> \verbatim +*> JOBU is CHARACTER*1 +*> = 'U': Unitary matrix U is computed; +*> = 'N': U is not computed. +*> \endverbatim +*> +*> \param[in] JOBV +*> \verbatim +*> JOBV is CHARACTER*1 +*> = 'V': Unitary matrix V is computed; +*> = 'N': V is not computed. +*> \endverbatim +*> +*> \param[in] JOBQ +*> \verbatim +*> JOBQ is CHARACTER*1 +*> = 'Q': Unitary matrix Q is computed; +*> = 'N': Q is not computed. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the M-by-N matrix A. +*> On exit, A contains the triangular (or trapezoidal) matrix +*> described in the Purpose section. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,M). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,N) +*> On entry, the P-by-N matrix B. +*> On exit, B contains the triangular matrix described in +*> the Purpose section. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,P). +*> \endverbatim +*> +*> \param[in] TOLA +*> \verbatim +*> TOLA is DOUBLE PRECISION +*> \endverbatim +*> +*> \param[in] TOLB +*> \verbatim +*> TOLB is DOUBLE PRECISION +*> +*> TOLA and TOLB are the thresholds to determine the effective +*> numerical rank of matrix B and a subblock of A. Generally, +*> they are set to +*> TOLA = MAX(M,N)*norm(A)*MAZHEPS, +*> TOLB = MAX(P,N)*norm(B)*MAZHEPS. +*> The size of TOLA and TOLB may affect the size of backward +*> errors of the decomposition. +*> \endverbatim +*> +*> \param[out] K +*> \verbatim +*> K is INTEGER +*> \endverbatim +*> +*> \param[out] L +*> \verbatim +*> L is INTEGER +*> +*> On exit, K and L specify the dimension of the subblocks +*> described in Purpose section. +*> K + L = effective numerical rank of (A**H,B**H)**H. +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX*16 array, dimension (LDU,M) +*> If JOBU = 'U', U contains the unitary matrix U. +*> If JOBU = 'N', U is not referenced. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M) if +*> JOBU = 'U'; LDU >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX*16 array, dimension (LDV,P) +*> If JOBV = 'V', V contains the unitary matrix V. +*> If JOBV = 'N', V is not referenced. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P) if +*> JOBV = 'V'; LDV >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is COMPLEX*16 array, dimension (LDQ,N) +*> If JOBQ = 'Q', Q contains the unitary matrix Q. +*> If JOBQ = 'N', Q is not referenced. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N) if +*> JOBQ = 'Q'; LDQ >= 1 otherwise. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX*16 array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) +*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is 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. +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup complex16OTHERcomputational +* +*> \par Further Details: +* ===================== +* +*> \verbatim +*> +*> The subroutine uses LAPACK subroutine ZGEQP3 for the QR factorization +*> with column pivoting to detect the effective numerical rank of the +*> a matrix. It may be replaced by a better rank determination strategy. +*> +*> ZGGSVP3 replaces the deprecated subroutine ZGGSVP. +*> +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE ZGGSVP3( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, + $ TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, + $ IWORK, RWORK, TAU, WORK, LWORK, INFO ) +* +* -- LAPACK computational routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* + IMPLICIT NONE +* +* .. Scalar Arguments .. + CHARACTER JOBQ, JOBU, JOBV + INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P, + $ LWORK + DOUBLE PRECISION TOLA, TOLB +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION RWORK( * ) + COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), + $ TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), + $ CONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL FORWRD, WANTQ, WANTU, WANTV, LQUERY + INTEGER I, J, LWKOPT + COMPLEX*16 T +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZGEQP3, ZGEQR2, ZGERQ2, ZLACPY, ZLAPMT, + $ ZLASET, ZUNG2R, ZUNM2R, ZUNMR2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DIMAG, MAX, MIN +* .. +* .. Executable Statements .. +* +* Test the input parameters +* + WANTU = LSAME( JOBU, 'U' ) + WANTV = LSAME( JOBV, 'V' ) + WANTQ = LSAME( JOBQ, 'Q' ) + FORWRD = .TRUE. + LQUERY = ( LWORK.EQ.-1 ) + LWKOPT = 1 +* +* Test the input arguments +* + INFO = 0 + IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN + INFO = -1 + ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN + INFO = -2 + ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( P.LT.0 ) THEN + INFO = -5 + ELSE IF( N.LT.0 ) THEN + INFO = -6 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -8 + ELSE IF( LDB.LT.MAX( 1, P ) ) THEN + INFO = -10 + ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN + INFO = -16 + ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN + INFO = -18 + ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN + INFO = -20 + ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN + INFO = -24 + END IF +* +* Compute workspace +* + IF( INFO.EQ.0 ) THEN + CALL ZGEQP3( P, N, B, LDB, IWORK, TAU, WORK, -1, RWORK, INFO ) + LWKOPT = INT( WORK ( 1 ) ) + IF( WANTV ) THEN + LWKOPT = MAX( LWKOPT, P ) + END IF + LWKOPT = MAX( LWKOPT, MIN( N, P ) ) + LWKOPT = MAX( LWKOPT, M ) + IF( WANTQ ) THEN + LWKOPT = MAX( LWKOPT, N ) + END IF + CALL ZGEQP3( M, N, A, LDA, IWORK, TAU, WORK, -1, RWORK, INFO ) + LWKOPT = MAX( LWKOPT, INT( WORK ( 1 ) ) ) + LWKOPT = MAX( 1, LWKOPT ) + WORK( 1 ) = DCMPLX( LWKOPT ) + END IF +* + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZGGSVP3', -INFO ) + RETURN + END IF + IF( LQUERY ) THEN + RETURN + ENDIF +* +* QR with column pivoting of B: B*P = V*( S11 S12 ) +* ( 0 0 ) +* + DO 10 I = 1, N + IWORK( I ) = 0 + 10 CONTINUE + CALL ZGEQP3( P, N, B, LDB, IWORK, TAU, WORK, LWORK, RWORK, INFO ) +* +* Update A := A*P +* + CALL ZLAPMT( FORWRD, M, N, A, LDA, IWORK ) +* +* Determine the effective rank of matrix B. +* + L = 0 + DO 20 I = 1, MIN( P, N ) + IF( ABS( B( I, I ) ).GT.TOLB ) + $ L = L + 1 + 20 CONTINUE +* + IF( WANTV ) THEN +* +* Copy the details of V, and form V. +* + CALL ZLASET( 'Full', P, P, CZERO, CZERO, V, LDV ) + IF( P.GT.1 ) + $ CALL ZLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ), + $ LDV ) + CALL ZUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO ) + END IF +* +* Clean up B +* + DO 40 J = 1, L - 1 + DO 30 I = J + 1, L + B( I, J ) = CZERO + 30 CONTINUE + 40 CONTINUE + IF( P.GT.L ) + $ CALL ZLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB ) +* + IF( WANTQ ) THEN +* +* Set Q = I and Update Q := Q*P +* + CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) + CALL ZLAPMT( FORWRD, N, N, Q, LDQ, IWORK ) + END IF +* + IF( P.GE.L .AND. N.NE.L ) THEN +* +* RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z +* + CALL ZGERQ2( L, N, B, LDB, TAU, WORK, INFO ) +* +* Update A := A*Z**H +* + CALL ZUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB, + $ TAU, A, LDA, WORK, INFO ) + IF( WANTQ ) THEN +* +* Update Q := Q*Z**H +* + CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N, L, B, + $ LDB, TAU, Q, LDQ, WORK, INFO ) + END IF +* +* Clean up B +* + CALL ZLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB ) + DO 60 J = N - L + 1, N + DO 50 I = J - N + L + 1, L + B( I, J ) = CZERO + 50 CONTINUE + 60 CONTINUE +* + END IF +* +* Let N-L L +* A = ( A11 A12 ) M, +* +* then the following does the complete QR decomposition of A11: +* +* A11 = U*( 0 T12 )*P1**H +* ( 0 0 ) +* + DO 70 I = 1, N - L + IWORK( I ) = 0 + 70 CONTINUE + CALL ZGEQP3( M, N-L, A, LDA, IWORK, TAU, WORK, LWORK, RWORK, + $ INFO ) +* +* Determine the effective rank of A11 +* + K = 0 + DO 80 I = 1, MIN( M, N-L ) + IF( ABS( A( I, I ) ).GT.TOLA ) + $ K = K + 1 + 80 CONTINUE +* +* Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N ) +* + CALL ZUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ), + $ A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Copy the details of U, and form U +* + CALL ZLASET( 'Full', M, M, CZERO, CZERO, U, LDU ) + IF( M.GT.1 ) + $ CALL ZLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ), + $ LDU ) + CALL ZUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO ) + END IF +* + IF( WANTQ ) THEN +* +* Update Q( 1:N, 1:N-L ) = Q( 1:N, 1:N-L )*P1 +* + CALL ZLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK ) + END IF +* +* Clean up A: set the strictly lower triangular part of +* A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0. +* + DO 100 J = 1, K - 1 + DO 90 I = J + 1, K + A( I, J ) = CZERO + 90 CONTINUE + 100 CONTINUE + IF( M.GT.K ) + $ CALL ZLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA ) +* + IF( N-L.GT.K ) THEN +* +* RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1 +* + CALL ZGERQ2( K, N-L, A, LDA, TAU, WORK, INFO ) +* + IF( WANTQ ) THEN +* +* Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H +* + CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A, + $ LDA, TAU, Q, LDQ, WORK, INFO ) + END IF +* +* Clean up A +* + CALL ZLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA ) + DO 120 J = N - L - K + 1, N - L + DO 110 I = J - N + L + K + 1, K + A( I, J ) = CZERO + 110 CONTINUE + 120 CONTINUE +* + END IF +* + IF( M.GT.K ) THEN +* +* QR factorization of A( K+1:M,N-L+1:N ) +* + CALL ZGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO ) +* + IF( WANTU ) THEN +* +* Update U(:,K+1:M) := U(:,K+1:M)*U1 +* + CALL ZUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ), + $ A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU, + $ WORK, INFO ) + END IF +* +* Clean up +* + DO 140 J = N - L + 1, N + DO 130 I = J - N + K + L + 1, M + A( I, J ) = CZERO + 130 CONTINUE + 140 CONTINUE +* + END IF +* + WORK( 1 ) = DCMPLX( LWKOPT ) + RETURN +* +* End of ZGGSVP3 +* + END diff --git a/TESTING/EIG/CMakeLists.txt b/TESTING/EIG/CMakeLists.txt index 574cf70b..1fcc73bd 100644 --- a/TESTING/EIG/CMakeLists.txt +++ b/TESTING/EIG/CMakeLists.txt @@ -58,7 +58,7 @@ set(SEIGTST schkee.f sget02.f sget10.f sget22.f sget23.f sget24.f sget31.f sget32.f sget33.f sget34.f sget35.f sget36.f sget37.f sget38.f sget39.f sget51.f sget52.f sget53.f - sget54.f sglmts.f sgqrts.f sgrqts.f sgsvts.f + sget54.f sglmts.f sgqrts.f sgrqts.f sgsvts.f sgsvts3.f shst01.f slarfy.f slarhs.f slatm4.f slctes.f slctsx.f slsets.f sort01.f sort03.f ssbt21.f ssgt01.f sslect.f sspt21.f sstt21.f sstt22.f ssyt21.f ssyt22.f) @@ -74,7 +74,7 @@ set(CEIGTST cchkee.f cerrbd.f cerrec.f cerred.f cerrgg.f cerrhs.f cerrst.f cget02.f cget10.f cget22.f cget23.f cget24.f cget35.f cget36.f cget37.f cget38.f cget51.f cget52.f - cget54.f cglmts.f cgqrts.f cgrqts.f cgsvts.f + cget54.f cglmts.f cgqrts.f cgrqts.f cgsvts.f cgsvts3.f chbt21.f chet21.f chet22.f chpt21.f chst01.f clarfy.f clarhs.f clatm4.f clctes.f clctsx.f clsets.f csbmv.f csgt01.f cslect.f @@ -95,7 +95,7 @@ set(DEIGTST dchkee.f dget02.f dget10.f dget22.f dget23.f dget24.f dget31.f dget32.f dget33.f dget34.f dget35.f dget36.f dget37.f dget38.f dget39.f dget51.f dget52.f dget53.f - dget54.f dglmts.f dgqrts.f dgrqts.f dgsvts.f + dget54.f dglmts.f dgqrts.f dgrqts.f dgsvts.f dgsvts3.f dhst01.f dlarfy.f dlarhs.f dlatm4.f dlctes.f dlctsx.f dlsets.f dort01.f dort03.f dsbt21.f dsgt01.f dslect.f dspt21.f dstt21.f dstt22.f dsyt21.f dsyt22.f) @@ -111,7 +111,7 @@ set(ZEIGTST zchkee.f zerrbd.f zerrec.f zerred.f zerrgg.f zerrhs.f zerrst.f zget02.f zget10.f zget22.f zget23.f zget24.f zget35.f zget36.f zget37.f zget38.f zget51.f zget52.f - zget54.f zglmts.f zgqrts.f zgrqts.f zgsvts.f + zget54.f zglmts.f zgqrts.f zgrqts.f zgsvts.f zgsvts3.f zhbt21.f zhet21.f zhet22.f zhpt21.f zhst01.f zlarfy.f zlarhs.f zlatm4.f zlctes.f zlctsx.f zlsets.f zsbmv.f zsgt01.f zslect.f diff --git a/TESTING/EIG/Makefile b/TESTING/EIG/Makefile index 41a14611..06543c2f 100644 --- a/TESTING/EIG/Makefile +++ b/TESTING/EIG/Makefile @@ -60,7 +60,7 @@ SEIGTST = schkee.o \ sget02.o sget10.o sget22.o sget23.o sget24.o sget31.o \ sget32.o sget33.o sget34.o sget35.o sget36.o \ sget37.o sget38.o sget39.o sget51.o sget52.o sget53.o \ - sget54.o sglmts.o sgqrts.o sgrqts.o sgsvts.o \ + sget54.o sglmts.o sgqrts.o sgrqts.o sgsvts.o sgsvts3.o \ shst01.o slarfy.o slarhs.o slatm4.o slctes.o slctsx.o slsets.o sort01.o \ sort03.o ssbt21.o ssgt01.o sslect.o sspt21.o sstt21.o \ sstt22.o ssyt21.o ssyt22.o @@ -76,7 +76,7 @@ CEIGTST = cchkee.o \ cerrbd.o cerrec.o cerred.o cerrgg.o cerrhs.o cerrst.o \ cget02.o cget10.o cget22.o cget23.o cget24.o \ cget35.o cget36.o cget37.o cget38.o cget51.o cget52.o \ - cget54.o cglmts.o cgqrts.o cgrqts.o cgsvts.o \ + cget54.o cglmts.o cgqrts.o cgrqts.o cgsvts.o cgsvts3.o \ chbt21.o chet21.o chet22.o chpt21.o chst01.o \ clarfy.o clarhs.o clatm4.o clctes.o clctsx.o clsets.o csbmv.o \ csgt01.o cslect.o \ @@ -97,7 +97,7 @@ DEIGTST = dchkee.o \ dget02.o dget10.o dget22.o dget23.o dget24.o dget31.o \ dget32.o dget33.o dget34.o dget35.o dget36.o \ dget37.o dget38.o dget39.o dget51.o dget52.o dget53.o \ - dget54.o dglmts.o dgqrts.o dgrqts.o dgsvts.o \ + dget54.o dglmts.o dgqrts.o dgrqts.o dgsvts.o dgsvts3.o \ dhst01.o dlarfy.o dlarhs.o dlatm4.o dlctes.o dlctsx.o dlsets.o dort01.o \ dort03.o dsbt21.o dsgt01.o dslect.o dspt21.o dstt21.o \ dstt22.o dsyt21.o dsyt22.o @@ -113,7 +113,7 @@ ZEIGTST = zchkee.o \ zerrbd.o zerrec.o zerred.o zerrgg.o zerrhs.o zerrst.o \ zget02.o zget10.o zget22.o zget23.o zget24.o \ zget35.o zget36.o zget37.o zget38.o zget51.o zget52.o \ - zget54.o zglmts.o zgqrts.o zgrqts.o zgsvts.o \ + zget54.o zglmts.o zgqrts.o zgrqts.o zgsvts.o zgsvts3.o \ zhbt21.o zhet21.o zhet22.o zhpt21.o zhst01.o \ zlarfy.o zlarhs.o zlatm4.o zlctes.o zlctsx.o zlsets.o zsbmv.o \ zsgt01.o zslect.o \ diff --git a/TESTING/EIG/cckgsv.f b/TESTING/EIG/cckgsv.f index fa189bfd..f11e6c75 100644 --- a/TESTING/EIG/cckgsv.f +++ b/TESTING/EIG/cckgsv.f @@ -219,7 +219,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 12 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) * .. @@ -237,7 +237,8 @@ REAL RESULT( NTESTS ) * .. * .. External Subroutines .. - EXTERNAL ALAHDG, ALAREQ, ALASUM, CGSVTS, CLATMS, SLATB9 + EXTERNAL ALAHDG, ALAREQ, ALASUM, CLATMS, SLATB9, CGSVTS, + $ CGSVTS3 * .. * .. Intrinsic Functions .. INTRINSIC ABS @@ -309,6 +310,12 @@ $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, $ LWORK, RWORK, RESULT ) * + CALL CGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT( NT+1 ) ) +* + NT = NT + 6 +* * Print information about the tests that did not * pass the threshold. * diff --git a/TESTING/EIG/cerrgg.f b/TESTING/EIG/cerrgg.f index 6f7e0508..70669b63 100644 --- a/TESTING/EIG/cerrgg.f +++ b/TESTING/EIG/cerrgg.f @@ -23,8 +23,8 @@ *> *> CERRGG tests the error exits for CGGES, CGGESX, CGGEV, CGGEVX, *> CGGES3, CGGEV3, CGGGLM, CGGHRD, CGGLSE, CGGQRF, CGGRQF, CGGSVD, -*> CGGSVP, CHGEQZ, CTGEVC, CTGEXC, CTGSEN, CTGSJA, CTGSNA, CTGSYL, -*> and CUNCSD. +*> CGGSVD3, CGGSVP, CGGSVP3, CHGEQZ, CTGEVC, CTGEXC, CTGSEN, CTGSJA, +*> CTGSNA, CTGSYL, and CUNCSD. *> \endverbatim * * Arguments: @@ -78,7 +78,7 @@ * .. Local Scalars .. CHARACTER*2 C2 INTEGER DUMMYK, DUMMYL, I, IFST, IHI, ILO, ILST, INFO, - $ J, M, NCYCLE, NT, SDIM + $ J, M, NCYCLE, NT, SDIM, LWORK REAL ANRM, BNRM, DIF, SCALE, TOLA, TOLB * .. * .. Local Arrays .. @@ -99,7 +99,8 @@ EXTERNAL CGGES, CGGESX, CGGEV, CGGEVX, CGGGLM, CGGHRD, $ CGGLSE, CGGQRF, CGGRQF, CGGSVD, CGGSVP, CHGEQZ, $ CHKXER, CTGEVC, CTGEXC, CTGSEN, CTGSJA, CTGSNA, - $ CTGSYL, CUNCSD, CGGES3, CGGEV3, CGGHD3 + $ CTGSYL, CUNCSD, CGGES3, CGGEV3, CGGHD3, + $ CGGSVD3, CGGSVP3 * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -135,6 +136,7 @@ IFST = 1 ILST = 1 NT = 0 + LWORK = 1 * * Test error exits for the GG path. * @@ -348,6 +350,126 @@ CALL CHKXER( 'CGGSVD', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* CGGSVD3 +* + SRNAMT = 'CGGSVD3' + INFOT = 1 + CALL CGGSVD3( '/', 'N', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CGGSVD3( 'N', '/', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL CGGSVD3( 'N', 'N', '/', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL CGGSVD3( 'N', 'N', 'N', -1, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL CGGSVD3( 'N', 'N', 'N', 0, -1, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL CGGSVD3( 'N', 'N', 'N', 0, 0, -1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL CGGSVD3( 'N', 'N', 'N', 2, 1, 1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 12 + CALL CGGSVD3( 'N', 'N', 'N', 1, 1, 2, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL CGGSVD3( 'U', 'N', 'N', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL CGGSVD3( 'N', 'V', 'N', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 2, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL CGGSVD3( 'N', 'N', 'Q', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 2, V, 2, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'CGGSVD3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* +* CGGSVP3 +* + SRNAMT = 'CGGSVP3' + INFOT = 1 + CALL CGGSVP3( '/', 'N', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL CGGSVP3( 'N', '/', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL CGGSVP3( 'N', 'N', '/', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL CGGSVP3( 'N', 'N', 'N', -1, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL CGGSVP3( 'N', 'N', 'N', 0, -1, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL CGGSVP3( 'N', 'N', 'N', 0, 0, -1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL CGGSVP3( 'N', 'N', 'N', 2, 1, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL CGGSVP3( 'N', 'N', 'N', 1, 2, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL CGGSVP3( 'U', 'N', 'N', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL CGGSVP3( 'N', 'V', 'N', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 2, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL CGGSVP3( 'N', 'N', 'Q', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 2, V, 2, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'CGGSVP3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * CGGSVP * SRNAMT = 'CGGSVP' diff --git a/TESTING/EIG/cgsvts3.f b/TESTING/EIG/cgsvts3.f new file mode 100644 index 00000000..db57a0bc --- /dev/null +++ b/TESTING/EIG/cgsvts3.f @@ -0,0 +1,396 @@ +*> \brief \b CGSVTS3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE CGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, +* LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, +* LWORK, RWORK, RESULT ) +* +* .. Scalar Arguments .. +* INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) +* COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ), +* $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), +* $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CGSVTS3 tests CGGSVD3, which computes the GSVD of an M-by-N matrix A +*> and a P-by-N matrix B: +*> U'*A*Q = D1*R and V'*B*Q = D2*R. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,M) +*> The M-by-N matrix A. +*> \endverbatim +*> +*> \param[out] AF +*> \verbatim +*> AF is COMPLEX array, dimension (LDA,N) +*> Details of the GSVD of A and B, as returned by CGGSVD3, +*> see CGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the arrays A and AF. +*> LDA >= max( 1,M ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is COMPLEX array, dimension (LDB,P) +*> On entry, the P-by-N matrix B. +*> \endverbatim +*> +*> \param[out] BF +*> \verbatim +*> BF is COMPLEX array, dimension (LDB,N) +*> Details of the GSVD of A and B, as returned by CGGSVD3, +*> see CGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the arrays B and BF. +*> LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX array, dimension(LDU,M) +*> The M by M unitary matrix U. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M). +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX array, dimension(LDV,M) +*> The P by P unitary matrix V. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P). +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is COMPLEX array, dimension(LDQ,N) +*> The N by N unitary matrix Q. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is REAL array, dimension (N) +*> +*> The generalized singular value pairs of A and B, the +*> ``diagonal'' matrices D1 and D2 are constructed from +*> ALPHA and BETA, see subroutine CGGSVD3 for details. +*> \endverbatim +*> +*> \param[out] R +*> \verbatim +*> R is COMPLEX array, dimension(LDQ,N) +*> The upper triangular matrix R. +*> \endverbatim +*> +*> \param[in] LDR +*> \verbatim +*> LDR is INTEGER +*> The leading dimension of the array R. LDR >= max(1,N). +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (LWORK) +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK, +*> LWORK >= max(M,P,N)*max(M,P,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (max(M,P,N)) +*> \endverbatim +*> +*> \param[out] RESULT +*> \verbatim +*> RESULT is REAL array, dimension (6) +*> The test ratios: +*> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) +*> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) +*> RESULT(3) = norm( I - U'*U ) / ( M*ULP ) +*> RESULT(4) = norm( I - V'*V ) / ( P*ULP ) +*> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) +*> RESULT(6) = 0 if ALPHA is in decreasing order; +*> = ULPINV otherwise. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup complex_eig +* +* ===================================================================== + SUBROUTINE CGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) + COMPLEX A( LDA, * ), AF( LDA, * ), B( LDB, * ), + $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), + $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + COMPLEX CZERO, CONE + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), + $ CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, INFO, J, K, L + REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL +* .. +* .. External Functions .. + REAL CLANGE, CLANHE, SLAMCH + EXTERNAL CLANGE, CLANHE, SLAMCH +* .. +* .. External Subroutines .. + EXTERNAL CGEMM, CGGSVD3, CHERK, CLACPY, CLASET, SCOPY +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN, REAL +* .. +* .. Executable Statements .. +* + ULP = SLAMCH( 'Precision' ) + ULPINV = ONE / ULP + UNFL = SLAMCH( 'Safe minimum' ) +* +* Copy the matrix A to the array AF. +* + CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB ) +* + ANORM = MAX( CLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) + BNORM = MAX( CLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) +* +* Factorize the matrices A and B in the arrays AF and BF. +* + CALL CGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, + $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, + $ RWORK, IWORK, INFO ) +* +* Copy R +* + DO 20 I = 1, MIN( K+L, M ) + DO 10 J = I, K + L + R( I, J ) = AF( I, N-K-L+J ) + 10 CONTINUE + 20 CONTINUE +* + IF( M-K-L.LT.0 ) THEN + DO 40 I = M + 1, K + L + DO 30 J = I, K + L + R( I, J ) = BF( I-K, N-K-L+J ) + 30 CONTINUE + 40 CONTINUE + END IF +* +* Compute A:= U'*A*Q - D1*R +* + CALL CGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA, + $ Q, LDQ, CZERO, WORK, LDA ) +* + CALL CGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE, + $ U, LDU, WORK, LDA, CZERO, A, LDA ) +* + DO 60 I = 1, K + DO 50 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) + 50 CONTINUE + 60 CONTINUE +* + DO 80 I = K + 1, MIN( K+L, M ) + DO 70 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) + 70 CONTINUE + 80 CONTINUE +* +* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . +* + RESID = CLANGE( '1', M, N, A, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) / + $ ULP + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute B := V'*B*Q - D2*R +* + CALL CGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB, + $ Q, LDQ, CZERO, WORK, LDB ) +* + CALL CGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE, + $ V, LDV, WORK, LDB, CZERO, B, LDB ) +* + DO 100 I = 1, L + DO 90 J = I, L + B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) + 90 CONTINUE + 100 CONTINUE +* +* Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . +* + RESID = CLANGE( '1', P, N, B, LDB, RWORK ) + IF( BNORM.GT.ZERO ) THEN + RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) / + $ ULP + ELSE + RESULT( 2 ) = ZERO + END IF +* +* Compute I - U'*U +* + CALL CLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ ) + CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU, + $ ONE, WORK, LDU ) +* +* Compute norm( I - U'*U ) / ( M * ULP ) . +* + RESID = CLANHE( '1', 'Upper', M, WORK, LDU, RWORK ) + RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP +* +* Compute I - V'*V +* + CALL CLASET( 'Full', P, P, CZERO, CONE, WORK, LDV ) + CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV, + $ ONE, WORK, LDV ) +* +* Compute norm( I - V'*V ) / ( P * ULP ) . +* + RESID = CLANHE( '1', 'Upper', P, WORK, LDV, RWORK ) + RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP +* +* Compute I - Q'*Q +* + CALL CLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ ) + CALL CHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ, + $ ONE, WORK, LDQ ) +* +* Compute norm( I - Q'*Q ) / ( N * ULP ) . +* + RESID = CLANHE( '1', 'Upper', N, WORK, LDQ, RWORK ) + RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP +* +* Check sorting +* + CALL SCOPY( N, ALPHA, 1, RWORK, 1 ) + DO 110 I = K + 1, MIN( K+L, M ) + J = IWORK( I ) + IF( I.NE.J ) THEN + TEMP = RWORK( I ) + RWORK( I ) = RWORK( J ) + RWORK( J ) = TEMP + END IF + 110 CONTINUE +* + RESULT( 6 ) = ZERO + DO 120 I = K + 1, MIN( K+L, M ) - 1 + IF( RWORK( I ).LT.RWORK( I+1 ) ) + $ RESULT( 6 ) = ULPINV + 120 CONTINUE +* + RETURN +* +* End of CGSVTS3 +* + END diff --git a/TESTING/EIG/dckgsv.f b/TESTING/EIG/dckgsv.f index cd9e5ae0..5dc7fcee 100644 --- a/TESTING/EIG/dckgsv.f +++ b/TESTING/EIG/dckgsv.f @@ -219,7 +219,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 12 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) * .. @@ -237,7 +237,8 @@ DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Subroutines .. - EXTERNAL ALAHDG, ALAREQ, ALASUM, DGSVTS, DLATB9, DLATMS + EXTERNAL ALAHDG, ALAREQ, ALASUM, DGSVTS, DLATB9, DLATMS, + $ DGSVTS3 * .. * .. Intrinsic Functions .. INTRINSIC ABS @@ -307,6 +308,12 @@ $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, $ LWORK, RWORK, RESULT ) * + CALL DGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT( NT+1 ) ) +* + NT = NT + 6 +* * Print information about the tests that did not * pass the threshold. * diff --git a/TESTING/EIG/derrgg.f b/TESTING/EIG/derrgg.f index e43ce273..b6492ff2 100644 --- a/TESTING/EIG/derrgg.f +++ b/TESTING/EIG/derrgg.f @@ -22,9 +22,9 @@ *> \verbatim *> *> DERRGG tests the error exits for DGGES, DGGESX, DGGEV, DGGEVX, -*> DGGGLM, DGGHRD, DGGLSE, DGGQRF, DGGRQF, DGGSVD, DGGSVP, DHGEQZ, -*> DORCSD, DTGEVC, DTGEXC, DTGSEN, DTGSJA, DTGSNA, DGGES3, DGGEV3, -*> and DTGSYL. +*> DGGGLM, DGGHRD, DGGLSE, DGGQRF, DGGRQF, DGGSVD, DGGSVD3, DGGSVP, +*> DGGSVP3, DHGEQZ, DORCSD, DTGEVC, DTGEXC, DTGSEN, DTGSJA, DTGSNA, +*> DGGES3, DGGEV3, and DTGSYL. *> \endverbatim * * Arguments: @@ -78,7 +78,7 @@ * .. Local Scalars .. CHARACTER*2 C2 INTEGER DUMMYK, DUMMYL, I, IFST, ILO, IHI, ILST, INFO, - $ J, M, NCYCLE, NT, SDIM + $ J, M, NCYCLE, NT, SDIM, LWORK DOUBLE PRECISION ANRM, BNRM, DIF, SCALE, TOLA, TOLB * .. * .. Local Arrays .. @@ -98,7 +98,8 @@ EXTERNAL CHKXER, DGGES, DGGESX, DGGEV, DGGEVX, DGGGLM, $ DGGHRD, DGGLSE, DGGQRF, DGGRQF, DGGSVD, DGGSVP, $ DHGEQZ, DORCSD, DTGEVC, DTGEXC, DTGSEN, DTGSJA, - $ DTGSNA, DTGSYL, DGGHD3, DGGES3, DGGEV3 + $ DTGSNA, DTGSYL, DGGHD3, DGGES3, DGGEV3, + $ DGGSVD3, DGGSVP3 * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -134,6 +135,7 @@ IFST = 1 ILST = 1 NT = 0 + LWORK = 1 * * Test error exits for the GG path. * @@ -347,6 +349,55 @@ CALL CHKXER( 'DGGSVD', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* DGGSVD3 +* + SRNAMT = 'DGGSVD3' + INFOT = 1 + CALL DGGSVD3( '/', 'N', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DGGSVD3( 'N', '/', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL DGGSVD3( 'N', 'N', '/', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL DGGSVD3( 'N', 'N', 'N', -1, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL DGGSVD3( 'N', 'N', 'N', 0, -1, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL DGGSVD3( 'N', 'N', 'N', 0, 0, -1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL DGGSVD3( 'N', 'N', 'N', 2, 1, 1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 12 + CALL DGGSVD3( 'N', 'N', 'N', 1, 1, 2, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL DGGSVD3( 'U', 'N', 'N', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL DGGSVD3( 'N', 'V', 'N', 1, 1, 2, DUMMYK, DUMMYL, A, 1, B, + $ 2, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL DGGSVD3( 'N', 'N', 'Q', 1, 2, 1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'DGGSVD3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * DGGSVP * SRNAMT = 'DGGSVP' @@ -407,6 +458,66 @@ CALL CHKXER( 'DGGSVP', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* DGGSVP3 +* + SRNAMT = 'DGGSVP3' + INFOT = 1 + CALL DGGSVP3( '/', 'N', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL DGGSVP3( 'N', '/', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL DGGSVP3( 'N', 'N', '/', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL DGGSVP3( 'N', 'N', 'N', -1, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL DGGSVP3( 'N', 'N', 'N', 0, -1, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL DGGSVP3( 'N', 'N', 'N', 0, 0, -1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL DGGSVP3( 'N', 'N', 'N', 2, 1, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL DGGSVP3( 'N', 'N', 'N', 1, 2, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL DGGSVP3( 'U', 'N', 'N', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL DGGSVP3( 'N', 'V', 'N', 1, 2, 1, A, 1, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL DGGSVP3( 'N', 'N', 'Q', 1, 1, 2, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'DGGSVP3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * DTGSJA * SRNAMT = 'DTGSJA' diff --git a/TESTING/EIG/dgsvts3.f b/TESTING/EIG/dgsvts3.f new file mode 100644 index 00000000..d7ceeb65 --- /dev/null +++ b/TESTING/EIG/dgsvts3.f @@ -0,0 +1,396 @@ +*> \brief \b DGSVTS3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE DGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, +* LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, +* LWORK, RWORK, RESULT ) +* +* .. Scalar Arguments .. +* INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), ALPHA( * ), +* $ B( LDB, * ), BETA( * ), BF( LDB, * ), +* $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), +* $ RWORK( * ), U( LDU, * ), V( LDV, * ), +* $ WORK( LWORK ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DGSVTS3 tests DGGSVD3, which computes the GSVD of an M-by-N matrix A +*> and a P-by-N matrix B: +*> U'*A*Q = D1*R and V'*B*Q = D2*R. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,M) +*> The M-by-N matrix A. +*> \endverbatim +*> +*> \param[out] AF +*> \verbatim +*> AF is DOUBLE PRECISION array, dimension (LDA,N) +*> Details of the GSVD of A and B, as returned by DGGSVD3, +*> see DGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the arrays A and AF. +*> LDA >= max( 1,M ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (LDB,P) +*> On entry, the P-by-N matrix B. +*> \endverbatim +*> +*> \param[out] BF +*> \verbatim +*> BF is DOUBLE PRECISION array, dimension (LDB,N) +*> Details of the GSVD of A and B, as returned by DGGSVD3, +*> see DGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the arrays B and BF. +*> LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is DOUBLE PRECISION array, dimension(LDU,M) +*> The M by M orthogonal matrix U. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M). +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension(LDV,M) +*> The P by P orthogonal matrix V. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P). +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is DOUBLE PRECISION array, dimension(LDQ,N) +*> The N by N orthogonal matrix Q. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION array, dimension (N) +*> +*> The generalized singular value pairs of A and B, the +*> ``diagonal'' matrices D1 and D2 are constructed from +*> ALPHA and BETA, see subroutine DGGSVD3 for details. +*> \endverbatim +*> +*> \param[out] R +*> \verbatim +*> R is DOUBLE PRECISION array, dimension(LDQ,N) +*> The upper triangular matrix R. +*> \endverbatim +*> +*> \param[in] LDR +*> \verbatim +*> LDR is INTEGER +*> The leading dimension of the array R. LDR >= max(1,N). +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (LWORK) +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK, +*> LWORK >= max(M,P,N)*max(M,P,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (max(M,P,N)) +*> \endverbatim +*> +*> \param[out] RESULT +*> \verbatim +*> RESULT is DOUBLE PRECISION array, dimension (6) +*> The test ratios: +*> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) +*> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) +*> RESULT(3) = norm( I - U'*U ) / ( M*ULP ) +*> RESULT(4) = norm( I - V'*V ) / ( P*ULP ) +*> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) +*> RESULT(6) = 0 if ALPHA is in decreasing order; +*> = ULPINV otherwise. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup double_eig +* +* ===================================================================== + SUBROUTINE DGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION A( LDA, * ), AF( LDA, * ), ALPHA( * ), + $ B( LDB, * ), BETA( * ), BF( LDB, * ), + $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), + $ RWORK( * ), U( LDU, * ), V( LDV, * ), + $ WORK( LWORK ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, INFO, J, K, L + DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLANGE, DLANSY + EXTERNAL DLAMCH, DLANGE, DLANSY +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DGEMM, DGGSVD3, DLACPY, DLASET, DSYRK +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, MAX, MIN +* .. +* .. Executable Statements .. +* + ULP = DLAMCH( 'Precision' ) + ULPINV = ONE / ULP + UNFL = DLAMCH( 'Safe minimum' ) +* +* Copy the matrix A to the array AF. +* + CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL DLACPY( 'Full', P, N, B, LDB, BF, LDB ) +* + ANORM = MAX( DLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) + BNORM = MAX( DLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) +* +* Factorize the matrices A and B in the arrays AF and BF. +* + CALL DGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, + $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, + $ IWORK, INFO ) +* +* Copy R +* + DO 20 I = 1, MIN( K+L, M ) + DO 10 J = I, K + L + R( I, J ) = AF( I, N-K-L+J ) + 10 CONTINUE + 20 CONTINUE +* + IF( M-K-L.LT.0 ) THEN + DO 40 I = M + 1, K + L + DO 30 J = I, K + L + R( I, J ) = BF( I-K, N-K-L+J ) + 30 CONTINUE + 40 CONTINUE + END IF +* +* Compute A:= U'*A*Q - D1*R +* + CALL DGEMM( 'No transpose', 'No transpose', M, N, N, ONE, A, LDA, + $ Q, LDQ, ZERO, WORK, LDA ) +* + CALL DGEMM( 'Transpose', 'No transpose', M, N, M, ONE, U, LDU, + $ WORK, LDA, ZERO, A, LDA ) +* + DO 60 I = 1, K + DO 50 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) + 50 CONTINUE + 60 CONTINUE +* + DO 80 I = K + 1, MIN( K+L, M ) + DO 70 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) + 70 CONTINUE + 80 CONTINUE +* +* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . +* + RESID = DLANGE( '1', M, N, A, LDA, RWORK ) +* + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / + $ ULP + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute B := V'*B*Q - D2*R +* + CALL DGEMM( 'No transpose', 'No transpose', P, N, N, ONE, B, LDB, + $ Q, LDQ, ZERO, WORK, LDB ) +* + CALL DGEMM( 'Transpose', 'No transpose', P, N, P, ONE, V, LDV, + $ WORK, LDB, ZERO, B, LDB ) +* + DO 100 I = 1, L + DO 90 J = I, L + B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) + 90 CONTINUE + 100 CONTINUE +* +* Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . +* + RESID = DLANGE( '1', P, N, B, LDB, RWORK ) + IF( BNORM.GT.ZERO ) THEN + RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, N ) ) ) / BNORM ) / + $ ULP + ELSE + RESULT( 2 ) = ZERO + END IF +* +* Compute I - U'*U +* + CALL DLASET( 'Full', M, M, ZERO, ONE, WORK, LDQ ) + CALL DSYRK( 'Upper', 'Transpose', M, M, -ONE, U, LDU, ONE, WORK, + $ LDU ) +* +* Compute norm( I - U'*U ) / ( M * ULP ) . +* + RESID = DLANSY( '1', 'Upper', M, WORK, LDU, RWORK ) + RESULT( 3 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / ULP +* +* Compute I - V'*V +* + CALL DLASET( 'Full', P, P, ZERO, ONE, WORK, LDV ) + CALL DSYRK( 'Upper', 'Transpose', P, P, -ONE, V, LDV, ONE, WORK, + $ LDV ) +* +* Compute norm( I - V'*V ) / ( P * ULP ) . +* + RESID = DLANSY( '1', 'Upper', P, WORK, LDV, RWORK ) + RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP +* +* Compute I - Q'*Q +* + CALL DLASET( 'Full', N, N, ZERO, ONE, WORK, LDQ ) + CALL DSYRK( 'Upper', 'Transpose', N, N, -ONE, Q, LDQ, ONE, WORK, + $ LDQ ) +* +* Compute norm( I - Q'*Q ) / ( N * ULP ) . +* + RESID = DLANSY( '1', 'Upper', N, WORK, LDQ, RWORK ) + RESULT( 5 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP +* +* Check sorting +* + CALL DCOPY( N, ALPHA, 1, WORK, 1 ) + DO 110 I = K + 1, MIN( K+L, M ) + J = IWORK( I ) + IF( I.NE.J ) THEN + TEMP = WORK( I ) + WORK( I ) = WORK( J ) + WORK( J ) = TEMP + END IF + 110 CONTINUE +* + RESULT( 6 ) = ZERO + DO 120 I = K + 1, MIN( K+L, M ) - 1 + IF( WORK( I ).LT.WORK( I+1 ) ) + $ RESULT( 6 ) = ULPINV + 120 CONTINUE +* + RETURN +* +* End of DGSVTS3 +* + END diff --git a/TESTING/EIG/sckgsv.f b/TESTING/EIG/sckgsv.f index ccddbc79..ecd18fa0 100644 --- a/TESTING/EIG/sckgsv.f +++ b/TESTING/EIG/sckgsv.f @@ -219,7 +219,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 12 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) * .. @@ -237,7 +237,8 @@ REAL RESULT( NTESTS ) * .. * .. External Subroutines .. - EXTERNAL ALAHDG, ALAREQ, ALASUM, SGSVTS, SLATB9, SLATMS + EXTERNAL ALAHDG, ALAREQ, ALASUM, SGSVTS, SLATB9, SLATMS, + $ SGSVTS3 * .. * .. Intrinsic Functions .. INTRINSIC ABS @@ -307,6 +308,12 @@ $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, $ LWORK, RWORK, RESULT ) * + CALL SGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT( NT+1 ) ) +* + NT = NT + 6 +* * Print information about the tests that did not * pass the threshold. * diff --git a/TESTING/EIG/serrgg.f b/TESTING/EIG/serrgg.f index 89610320..2973f171 100644 --- a/TESTING/EIG/serrgg.f +++ b/TESTING/EIG/serrgg.f @@ -23,8 +23,8 @@ *> *> SERRGG tests the error exits for SGGES, SGGESX, SGGEV, SGGEVX, *> SGGES3, SGGEV3, SGGGLM, SGGHRD, SGGLSE, SGGQRF, SGGRQF, SGGSVD, -*> SGGSVP, SHGEQZ, SORCSD, STGEVC, STGEXC, STGSEN, STGSJA, STGSNA, -*> and STGSYL. +*> SGGSVD3, SGGSVP, SGGSVP3, SHGEQZ, SORCSD, STGEVC, STGEXC, STGSEN, +*> STGSJA, STGSNA, and STGSYL. *> \endverbatim * * Arguments: @@ -78,7 +78,7 @@ * .. Local Scalars .. CHARACTER*2 C2 INTEGER DUMMYK, DUMMYL, I, IFST, ILO, IHI, ILST, INFO, - $ J, M, NCYCLE, NT, SDIM + $ J, M, NCYCLE, NT, SDIM, LWORK REAL ANRM, BNRM, DIF, SCALE, TOLA, TOLB * .. * .. Local Arrays .. @@ -98,7 +98,8 @@ EXTERNAL CHKXER, SGGES, SGGESX, SGGEV, SGGEVX, SGGGLM, $ SGGHRD, SGGLSE, SGGQRF, SGGRQF, SGGSVD, SGGSVP, $ SHGEQZ, SORCSD, STGEVC, STGEXC, STGSEN, STGSJA, - $ STGSNA, STGSYL, SGGES3, SGGEV3, SGGHD3 + $ STGSNA, STGSYL, SGGES3, SGGEV3, SGGHD3, + $ SGGSVD3, SGGSVP3 * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -134,6 +135,7 @@ IFST = 1 ILST = 1 NT = 0 + LWORK = 1 * * Test error exits for the GG path. * @@ -347,6 +349,55 @@ CALL CHKXER( 'SGGSVD', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* SGGSVD3 +* + SRNAMT = 'SGGSVD3' + INFOT = 1 + CALL SGGSVD3( '/', 'N', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL SGGSVD3( 'N', '/', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL SGGSVD3( 'N', 'N', '/', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL SGGSVD3( 'N', 'N', 'N', -1, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL SGGSVD3( 'N', 'N', 'N', 0, -1, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL SGGSVD3( 'N', 'N', 'N', 0, 0, -1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL SGGSVD3( 'N', 'N', 'N', 2, 1, 1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 12 + CALL SGGSVD3( 'N', 'N', 'N', 1, 1, 2, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL SGGSVD3( 'U', 'N', 'N', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL SGGSVD3( 'N', 'V', 'N', 1, 1, 2, DUMMYK, DUMMYL, A, 1, B, + $ 2, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL SGGSVD3( 'N', 'N', 'Q', 1, 2, 1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, IW, LWORK, INFO ) + CALL CHKXER( 'SGGSVD3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * SGGSVP * SRNAMT = 'SGGSVP' @@ -407,6 +458,66 @@ CALL CHKXER( 'SGGSVP', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* SGGSVP3 +* + SRNAMT = 'SGGSVP3' + INFOT = 1 + CALL SGGSVP3( '/', 'N', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL SGGSVP3( 'N', '/', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL SGGSVP3( 'N', 'N', '/', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL SGGSVP3( 'N', 'N', 'N', -1, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL SGGSVP3( 'N', 'N', 'N', 0, -1, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL SGGSVP3( 'N', 'N', 'N', 0, 0, -1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL SGGSVP3( 'N', 'N', 'N', 2, 1, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL SGGSVP3( 'N', 'N', 'N', 1, 2, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL SGGSVP3( 'U', 'N', 'N', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL SGGSVP3( 'N', 'V', 'N', 1, 2, 1, A, 1, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL SGGSVP3( 'N', 'N', 'Q', 1, 1, 2, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'SGGSVP3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * STGSJA * SRNAMT = 'STGSJA' diff --git a/TESTING/EIG/sgsvts3.f b/TESTING/EIG/sgsvts3.f new file mode 100644 index 00000000..23aa62c7 --- /dev/null +++ b/TESTING/EIG/sgsvts3.f @@ -0,0 +1,396 @@ +*> \brief \b SGSVTS3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE SGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, +* LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, +* LWORK, RWORK, RESULT ) +* +* .. Scalar Arguments .. +* INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* REAL A( LDA, * ), AF( LDA, * ), ALPHA( * ), +* $ B( LDB, * ), BETA( * ), BF( LDB, * ), +* $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), +* $ RWORK( * ), U( LDU, * ), V( LDV, * ), +* $ WORK( LWORK ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SGSVTS3 tests SGGSVD3, which computes the GSVD of an M-by-N matrix A +*> and a P-by-N matrix B: +*> U'*A*Q = D1*R and V'*B*Q = D2*R. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is REAL array, dimension (LDA,M) +*> The M-by-N matrix A. +*> \endverbatim +*> +*> \param[out] AF +*> \verbatim +*> AF is REAL array, dimension (LDA,N) +*> Details of the GSVD of A and B, as returned by SGGSVD3, +*> see SGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the arrays A and AF. +*> LDA >= max( 1,M ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is REAL array, dimension (LDB,P) +*> On entry, the P-by-N matrix B. +*> \endverbatim +*> +*> \param[out] BF +*> \verbatim +*> BF is REAL array, dimension (LDB,N) +*> Details of the GSVD of A and B, as returned by SGGSVD3, +*> see SGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the arrays B and BF. +*> LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is REAL array, dimension(LDU,M) +*> The M by M orthogonal matrix U. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M). +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is REAL array, dimension(LDV,M) +*> The P by P orthogonal matrix V. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P). +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is REAL array, dimension(LDQ,N) +*> The N by N orthogonal matrix Q. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is REAL array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is REAL array, dimension (N) +*> +*> The generalized singular value pairs of A and B, the +*> ``diagonal'' matrices D1 and D2 are constructed from +*> ALPHA and BETA, see subroutine SGGSVD3 for details. +*> \endverbatim +*> +*> \param[out] R +*> \verbatim +*> R is REAL array, dimension(LDQ,N) +*> The upper triangular matrix R. +*> \endverbatim +*> +*> \param[in] LDR +*> \verbatim +*> LDR is INTEGER +*> The leading dimension of the array R. LDR >= max(1,N). +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (LWORK) +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK, +*> LWORK >= max(M,P,N)*max(M,P,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is REAL array, dimension (max(M,P,N)) +*> \endverbatim +*> +*> \param[out] RESULT +*> \verbatim +*> RESULT is REAL array, dimension (6) +*> The test ratios: +*> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) +*> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) +*> RESULT(3) = norm( I - U'*U ) / ( M*ULP ) +*> RESULT(4) = norm( I - V'*V ) / ( P*ULP ) +*> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) +*> RESULT(6) = 0 if ALPHA is in decreasing order; +*> = ULPINV otherwise. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup single_eig +* +* ===================================================================== + SUBROUTINE SGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + REAL A( LDA, * ), AF( LDA, * ), ALPHA( * ), + $ B( LDB, * ), BETA( * ), BF( LDB, * ), + $ Q( LDQ, * ), R( LDR, * ), RESULT( 6 ), + $ RWORK( * ), U( LDU, * ), V( LDV, * ), + $ WORK( LWORK ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, INFO, J, K, L + REAL ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL +* .. +* .. External Functions .. + REAL SLAMCH, SLANGE, SLANSY + EXTERNAL SLAMCH, SLANGE, SLANSY +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SGEMM, SGGSVD3, SLACPY, SLASET, SSYRK +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN, REAL +* .. +* .. Executable Statements .. +* + ULP = SLAMCH( 'Precision' ) + ULPINV = ONE / ULP + UNFL = SLAMCH( 'Safe minimum' ) +* +* Copy the matrix A to the array AF. +* + CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB ) +* + ANORM = MAX( SLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) + BNORM = MAX( SLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) +* +* Factorize the matrices A and B in the arrays AF and BF. +* + CALL SGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, + $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, + $ IWORK, INFO ) +* +* Copy R +* + DO 20 I = 1, MIN( K+L, M ) + DO 10 J = I, K + L + R( I, J ) = AF( I, N-K-L+J ) + 10 CONTINUE + 20 CONTINUE +* + IF( M-K-L.LT.0 ) THEN + DO 40 I = M + 1, K + L + DO 30 J = I, K + L + R( I, J ) = BF( I-K, N-K-L+J ) + 30 CONTINUE + 40 CONTINUE + END IF +* +* Compute A:= U'*A*Q - D1*R +* + CALL SGEMM( 'No transpose', 'No transpose', M, N, N, ONE, A, LDA, + $ Q, LDQ, ZERO, WORK, LDA ) +* + CALL SGEMM( 'Transpose', 'No transpose', M, N, M, ONE, U, LDU, + $ WORK, LDA, ZERO, A, LDA ) +* + DO 60 I = 1, K + DO 50 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) + 50 CONTINUE + 60 CONTINUE +* + DO 80 I = K + 1, MIN( K+L, M ) + DO 70 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) + 70 CONTINUE + 80 CONTINUE +* +* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . +* + RESID = SLANGE( '1', M, N, A, LDA, RWORK ) +* + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / REAL( MAX( 1, M, N ) ) ) / ANORM ) / + $ ULP + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute B := V'*B*Q - D2*R +* + CALL SGEMM( 'No transpose', 'No transpose', P, N, N, ONE, B, LDB, + $ Q, LDQ, ZERO, WORK, LDB ) +* + CALL SGEMM( 'Transpose', 'No transpose', P, N, P, ONE, V, LDV, + $ WORK, LDB, ZERO, B, LDB ) +* + DO 100 I = 1, L + DO 90 J = I, L + B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) + 90 CONTINUE + 100 CONTINUE +* +* Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . +* + RESID = SLANGE( '1', P, N, B, LDB, RWORK ) + IF( BNORM.GT.ZERO ) THEN + RESULT( 2 ) = ( ( RESID / REAL( MAX( 1, P, N ) ) ) / BNORM ) / + $ ULP + ELSE + RESULT( 2 ) = ZERO + END IF +* +* Compute I - U'*U +* + CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDQ ) + CALL SSYRK( 'Upper', 'Transpose', M, M, -ONE, U, LDU, ONE, WORK, + $ LDU ) +* +* Compute norm( I - U'*U ) / ( M * ULP ) . +* + RESID = SLANSY( '1', 'Upper', M, WORK, LDU, RWORK ) + RESULT( 3 ) = ( RESID / REAL( MAX( 1, M ) ) ) / ULP +* +* Compute I - V'*V +* + CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDV ) + CALL SSYRK( 'Upper', 'Transpose', P, P, -ONE, V, LDV, ONE, WORK, + $ LDV ) +* +* Compute norm( I - V'*V ) / ( P * ULP ) . +* + RESID = SLANSY( '1', 'Upper', P, WORK, LDV, RWORK ) + RESULT( 4 ) = ( RESID / REAL( MAX( 1, P ) ) ) / ULP +* +* Compute I - Q'*Q +* + CALL SLASET( 'Full', N, N, ZERO, ONE, WORK, LDQ ) + CALL SSYRK( 'Upper', 'Transpose', N, N, -ONE, Q, LDQ, ONE, WORK, + $ LDQ ) +* +* Compute norm( I - Q'*Q ) / ( N * ULP ) . +* + RESID = SLANSY( '1', 'Upper', N, WORK, LDQ, RWORK ) + RESULT( 5 ) = ( RESID / REAL( MAX( 1, N ) ) ) / ULP +* +* Check sorting +* + CALL SCOPY( N, ALPHA, 1, WORK, 1 ) + DO 110 I = K + 1, MIN( K+L, M ) + J = IWORK( I ) + IF( I.NE.J ) THEN + TEMP = WORK( I ) + WORK( I ) = WORK( J ) + WORK( J ) = TEMP + END IF + 110 CONTINUE +* + RESULT( 6 ) = ZERO + DO 120 I = K + 1, MIN( K+L, M ) - 1 + IF( WORK( I ).LT.WORK( I+1 ) ) + $ RESULT( 6 ) = ULPINV + 120 CONTINUE +* + RETURN +* +* End of SGSVTS3 +* + END diff --git a/TESTING/EIG/zckgsv.f b/TESTING/EIG/zckgsv.f index c36373e3..e6fa7365 100644 --- a/TESTING/EIG/zckgsv.f +++ b/TESTING/EIG/zckgsv.f @@ -219,7 +219,7 @@ * * .. Parameters .. INTEGER NTESTS - PARAMETER ( NTESTS = 7 ) + PARAMETER ( NTESTS = 12 ) INTEGER NTYPES PARAMETER ( NTYPES = 8 ) * .. @@ -237,7 +237,8 @@ DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Subroutines .. - EXTERNAL ALAHDG, ALAREQ, ALASUM, DLATB9, ZGSVTS, ZLATMS + EXTERNAL ALAHDG, ALAREQ, ALASUM, DLATB9, ZGSVTS, ZLATMS, + $ ZGSVTS3 * .. * .. Intrinsic Functions .. INTRINSIC ABS @@ -309,6 +310,12 @@ $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, $ LWORK, RWORK, RESULT ) * + CALL ZGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT( NT+1 ) ) +* + NT = NT + 6 +* * Print information about the tests that did not * pass the threshold. * diff --git a/TESTING/EIG/zerrgg.f b/TESTING/EIG/zerrgg.f index 5ed7ee61..420b6241 100644 --- a/TESTING/EIG/zerrgg.f +++ b/TESTING/EIG/zerrgg.f @@ -23,8 +23,8 @@ *> *> ZERRGG tests the error exits for ZGGES, ZGGESX, ZGGEV, ZGGEVX, *> ZGGES3, ZGGEV3, ZGGGLM, ZGGHRD, ZGGLSE, ZGGQRF, ZGGRQF, ZGGSVD, -*> ZGGSVP, ZHGEQZ, ZTGEVC, ZTGEXC, ZTGSEN, ZTGSJA, ZTGSNA, ZTGSYL, -*> and ZUNCSD. +*> ZGGSVD3, ZGGSVP, ZGGSVP3, ZHGEQZ, ZTGEVC, ZTGEXC, ZTGSEN, ZTGSJA, +*> ZTGSNA, ZTGSYL, and ZUNCSD. *> \endverbatim * * Arguments: @@ -78,7 +78,7 @@ * .. Local Scalars .. CHARACTER*2 C2 INTEGER DUMMYK, DUMMYL, I, IFST, IHI, ILO, ILST, INFO, - $ J, M, NCYCLE, NT, SDIM + $ J, M, NCYCLE, NT, SDIM, LWORK DOUBLE PRECISION ANRM, BNRM, DIF, SCALE, TOLA, TOLB * .. * .. Local Arrays .. @@ -99,7 +99,8 @@ EXTERNAL CHKXER, ZGGES, ZGGESX, ZGGEV, ZGGEVX, ZGGGLM, $ ZGGHRD, ZGGLSE, ZGGQRF, ZGGRQF, ZGGSVD, ZGGSVP, $ ZHGEQZ, ZTGEVC, ZTGEXC, ZTGSEN, ZTGSJA, ZTGSNA, - $ ZTGSYL, ZUNCSD, ZGGES3, ZGGEV3, ZGGHD3 + $ ZTGSYL, ZUNCSD, ZGGES3, ZGGEV3, ZGGHD3, + $ ZGGSVD3, ZGGSVP3 * .. * .. Scalars in Common .. LOGICAL LERR, OK @@ -135,6 +136,7 @@ IFST = 1 ILST = 1 NT = 0 + LWORK = 1 * * Test error exits for the GG path. * @@ -348,6 +350,66 @@ CALL CHKXER( 'ZGGSVD', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* ZGGSVD3 +* + SRNAMT = 'ZGGSVD3' + INFOT = 1 + CALL ZGGSVD3( '/', 'N', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZGGSVD3( 'N', '/', 'N', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL ZGGSVD3( 'N', 'N', '/', 0, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL ZGGSVD3( 'N', 'N', 'N', -1, 0, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL ZGGSVD3( 'N', 'N', 'N', 0, -1, 0, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL ZGGSVD3( 'N', 'N', 'N', 0, 0, -1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL ZGGSVD3( 'N', 'N', 'N', 2, 1, 1, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 12 + CALL ZGGSVD3( 'N', 'N', 'N', 1, 1, 2, DUMMYK, DUMMYL, A, 1, B, + $ 1, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL ZGGSVD3( 'U', 'N', 'N', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 1, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL ZGGSVD3( 'N', 'V', 'N', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 2, V, 1, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL ZGGSVD3( 'N', 'N', 'Q', 2, 2, 2, DUMMYK, DUMMYL, A, 2, B, + $ 2, R1, R2, U, 2, V, 2, Q, 1, W, RW, IW, LWORK, + $ INFO ) + CALL CHKXER( 'ZGGSVD3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * ZGGSVP * SRNAMT = 'ZGGSVP' @@ -408,6 +470,66 @@ CALL CHKXER( 'ZGGSVP', INFOT, NOUT, LERR, OK ) NT = NT + 11 * +* ZGGSVP3 +* + SRNAMT = 'ZGGSVP3' + INFOT = 1 + CALL ZGGSVP3( '/', 'N', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 2 + CALL ZGGSVP3( 'N', '/', 'N', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 3 + CALL ZGGSVP3( 'N', 'N', '/', 0, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 4 + CALL ZGGSVP3( 'N', 'N', 'N', -1, 0, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 5 + CALL ZGGSVP3( 'N', 'N', 'N', 0, -1, 0, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 6 + CALL ZGGSVP3( 'N', 'N', 'N', 0, 0, -1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 8 + CALL ZGGSVP3( 'N', 'N', 'N', 2, 1, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 10 + CALL ZGGSVP3( 'N', 'N', 'N', 1, 2, 1, A, 1, B, 1, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 16 + CALL ZGGSVP3( 'U', 'N', 'N', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 1, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 18 + CALL ZGGSVP3( 'N', 'V', 'N', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 2, V, 1, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + INFOT = 20 + CALL ZGGSVP3( 'N', 'N', 'Q', 2, 2, 2, A, 2, B, 2, TOLA, TOLB, + $ DUMMYK, DUMMYL, U, 2, V, 2, Q, 1, IW, RW, TAU, W, + $ LWORK, INFO ) + CALL CHKXER( 'ZGGSVP3', INFOT, NOUT, LERR, OK ) + NT = NT + 11 +* * ZTGSJA * SRNAMT = 'ZTGSJA' diff --git a/TESTING/EIG/zgsvts3.f b/TESTING/EIG/zgsvts3.f new file mode 100644 index 00000000..8a5c2d4c --- /dev/null +++ b/TESTING/EIG/zgsvts3.f @@ -0,0 +1,396 @@ +*> \brief \b ZGSVTS3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +* Definition: +* =========== +* +* SUBROUTINE ZGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, +* LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, +* LWORK, RWORK, RESULT ) +* +* .. Scalar Arguments .. +* INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) +* COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), +* $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), +* $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZGSVTS3 tests ZGGSVD3, which computes the GSVD of an M-by-N matrix A +*> and a P-by-N matrix B: +*> U'*A*Q = D1*R and V'*B*Q = D2*R. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix A. M >= 0. +*> \endverbatim +*> +*> \param[in] P +*> \verbatim +*> P is INTEGER +*> The number of rows of the matrix B. P >= 0. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,M) +*> The M-by-N matrix A. +*> \endverbatim +*> +*> \param[out] AF +*> \verbatim +*> AF is COMPLEX*16 array, dimension (LDA,N) +*> Details of the GSVD of A and B, as returned by ZGGSVD3, +*> see ZGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the arrays A and AF. +*> LDA >= max( 1,M ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is COMPLEX*16 array, dimension (LDB,P) +*> On entry, the P-by-N matrix B. +*> \endverbatim +*> +*> \param[out] BF +*> \verbatim +*> BF is COMPLEX*16 array, dimension (LDB,N) +*> Details of the GSVD of A and B, as returned by ZGGSVD3, +*> see ZGGSVD3 for further details. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the arrays B and BF. +*> LDB >= max(1,P). +*> \endverbatim +*> +*> \param[out] U +*> \verbatim +*> U is COMPLEX*16 array, dimension(LDU,M) +*> The M by M unitary matrix U. +*> \endverbatim +*> +*> \param[in] LDU +*> \verbatim +*> LDU is INTEGER +*> The leading dimension of the array U. LDU >= max(1,M). +*> \endverbatim +*> +*> \param[out] V +*> \verbatim +*> V is COMPLEX*16 array, dimension(LDV,M) +*> The P by P unitary matrix V. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. LDV >= max(1,P). +*> \endverbatim +*> +*> \param[out] Q +*> \verbatim +*> Q is COMPLEX*16 array, dimension(LDQ,N) +*> The N by N unitary matrix Q. +*> \endverbatim +*> +*> \param[in] LDQ +*> \verbatim +*> LDQ is INTEGER +*> The leading dimension of the array Q. LDQ >= max(1,N). +*> \endverbatim +*> +*> \param[out] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION array, dimension (N) +*> \endverbatim +*> +*> \param[out] BETA +*> \verbatim +*> BETA is DOUBLE PRECISION array, dimension (N) +*> +*> The generalized singular value pairs of A and B, the +*> ``diagonal'' matrices D1 and D2 are constructed from +*> ALPHA and BETA, see subroutine ZGGSVD3 for details. +*> \endverbatim +*> +*> \param[out] R +*> \verbatim +*> R is COMPLEX*16 array, dimension(LDQ,N) +*> The upper triangular matrix R. +*> \endverbatim +*> +*> \param[in] LDR +*> \verbatim +*> LDR is INTEGER +*> The leading dimension of the array R. LDR >= max(1,N). +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (N) +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (LWORK) +*> \endverbatim +*> +*> \param[in] LWORK +*> \verbatim +*> LWORK is INTEGER +*> The dimension of the array WORK, +*> LWORK >= max(M,P,N)*max(M,P,N). +*> \endverbatim +*> +*> \param[out] RWORK +*> \verbatim +*> RWORK is DOUBLE PRECISION array, dimension (max(M,P,N)) +*> \endverbatim +*> +*> \param[out] RESULT +*> \verbatim +*> RESULT is DOUBLE PRECISION array, dimension (6) +*> The test ratios: +*> RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) +*> RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) +*> RESULT(3) = norm( I - U'*U ) / ( M*ULP ) +*> RESULT(4) = norm( I - V'*V ) / ( P*ULP ) +*> RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) +*> RESULT(6) = 0 if ALPHA is in decreasing order; +*> = ULPINV otherwise. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date August 2015 +* +*> \ingroup complex16_eig +* +* ===================================================================== + SUBROUTINE ZGSVTS3( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V, + $ LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK, + $ LWORK, RWORK, RESULT ) +* +* -- LAPACK test routine (version 3.6.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* August 2015 +* +* .. Scalar Arguments .. + INTEGER LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P +* .. +* .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION ALPHA( * ), BETA( * ), RESULT( 6 ), RWORK( * ) + COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ), + $ BF( LDB, * ), Q( LDQ, * ), R( LDR, * ), + $ U( LDU, * ), V( LDV, * ), WORK( LWORK ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + COMPLEX*16 CZERO, CONE + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), + $ CONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, INFO, J, K, L + DOUBLE PRECISION ANORM, BNORM, RESID, TEMP, ULP, ULPINV, UNFL +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE + EXTERNAL DLAMCH, ZLANGE, ZLANHE +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, ZGEMM, ZGGSVD3, ZHERK, ZLACPY, ZLASET +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, MAX, MIN +* .. +* .. Executable Statements .. +* + ULP = DLAMCH( 'Precision' ) + ULPINV = ONE / ULP + UNFL = DLAMCH( 'Safe minimum' ) +* +* Copy the matrix A to the array AF. +* + CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL ZLACPY( 'Full', P, N, B, LDB, BF, LDB ) +* + ANORM = MAX( ZLANGE( '1', M, N, A, LDA, RWORK ), UNFL ) + BNORM = MAX( ZLANGE( '1', P, N, B, LDB, RWORK ), UNFL ) +* +* Factorize the matrices A and B in the arrays AF and BF. +* + CALL ZGGSVD3( 'U', 'V', 'Q', M, N, P, K, L, AF, LDA, BF, LDB, + $ ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, + $ RWORK, IWORK, INFO ) +* +* Copy R +* + DO 20 I = 1, MIN( K+L, M ) + DO 10 J = I, K + L + R( I, J ) = AF( I, N-K-L+J ) + 10 CONTINUE + 20 CONTINUE +* + IF( M-K-L.LT.0 ) THEN + DO 40 I = M + 1, K + L + DO 30 J = I, K + L + R( I, J ) = BF( I-K, N-K-L+J ) + 30 CONTINUE + 40 CONTINUE + END IF +* +* Compute A:= U'*A*Q - D1*R +* + CALL ZGEMM( 'No transpose', 'No transpose', M, N, N, CONE, A, LDA, + $ Q, LDQ, CZERO, WORK, LDA ) +* + CALL ZGEMM( 'Conjugate transpose', 'No transpose', M, N, M, CONE, + $ U, LDU, WORK, LDA, CZERO, A, LDA ) +* + DO 60 I = 1, K + DO 50 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - R( I, J ) + 50 CONTINUE + 60 CONTINUE +* + DO 80 I = K + 1, MIN( K+L, M ) + DO 70 J = I, K + L + A( I, N-K-L+J ) = A( I, N-K-L+J ) - ALPHA( I )*R( I, J ) + 70 CONTINUE + 80 CONTINUE +* +* Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . +* + RESID = ZLANGE( '1', M, N, A, LDA, RWORK ) + IF( ANORM.GT.ZERO ) THEN + RESULT( 1 ) = ( ( RESID / DBLE( MAX( 1, M, N ) ) ) / ANORM ) / + $ ULP + ELSE + RESULT( 1 ) = ZERO + END IF +* +* Compute B := V'*B*Q - D2*R +* + CALL ZGEMM( 'No transpose', 'No transpose', P, N, N, CONE, B, LDB, + $ Q, LDQ, CZERO, WORK, LDB ) +* + CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, N, P, CONE, + $ V, LDV, WORK, LDB, CZERO, B, LDB ) +* + DO 100 I = 1, L + DO 90 J = I, L + B( I, N-L+J ) = B( I, N-L+J ) - BETA( K+I )*R( K+I, K+J ) + 90 CONTINUE + 100 CONTINUE +* +* Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . +* + RESID = ZLANGE( '1', P, N, B, LDB, RWORK ) + IF( BNORM.GT.ZERO ) THEN + RESULT( 2 ) = ( ( RESID / DBLE( MAX( 1, P, N ) ) ) / BNORM ) / + $ ULP + ELSE + RESULT( 2 ) = ZERO + END IF +* +* Compute I - U'*U +* + CALL ZLASET( 'Full', M, M, CZERO, CONE, WORK, LDQ ) + CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, U, LDU, + $ ONE, WORK, LDU ) +* +* Compute norm( I - U'*U ) / ( M * ULP ) . +* + RESID = ZLANHE( '1', 'Upper', M, WORK, LDU, RWORK ) + RESULT( 3 ) = ( RESID / DBLE( MAX( 1, M ) ) ) / ULP +* +* Compute I - V'*V +* + CALL ZLASET( 'Full', P, P, CZERO, CONE, WORK, LDV ) + CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, V, LDV, + $ ONE, WORK, LDV ) +* +* Compute norm( I - V'*V ) / ( P * ULP ) . +* + RESID = ZLANHE( '1', 'Upper', P, WORK, LDV, RWORK ) + RESULT( 4 ) = ( RESID / DBLE( MAX( 1, P ) ) ) / ULP +* +* Compute I - Q'*Q +* + CALL ZLASET( 'Full', N, N, CZERO, CONE, WORK, LDQ ) + CALL ZHERK( 'Upper', 'Conjugate transpose', N, N, -ONE, Q, LDQ, + $ ONE, WORK, LDQ ) +* +* Compute norm( I - Q'*Q ) / ( N * ULP ) . +* + RESID = ZLANHE( '1', 'Upper', N, WORK, LDQ, RWORK ) + RESULT( 5 ) = ( RESID / DBLE( MAX( 1, N ) ) ) / ULP +* +* Check sorting +* + CALL DCOPY( N, ALPHA, 1, RWORK, 1 ) + DO 110 I = K + 1, MIN( K+L, M ) + J = IWORK( I ) + IF( I.NE.J ) THEN + TEMP = RWORK( I ) + RWORK( I ) = RWORK( J ) + RWORK( J ) = TEMP + END IF + 110 CONTINUE +* + RESULT( 6 ) = ZERO + DO 120 I = K + 1, MIN( K+L, M ) - 1 + IF( RWORK( I ).LT.RWORK( I+1 ) ) + $ RESULT( 6 ) = ULPINV + 120 CONTINUE +* + RETURN +* +* End of ZGSVTS3 +* + END |