TODO: LAPACKE wrappers.
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
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
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
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
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 \
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 \
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 \
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 \
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
RETURN
END IF
*
-* Quick return if possible.
-*
- IF( MINMN.EQ.0 ) THEN
- RETURN
- END IF
-*
* Move initial columns up front.
*
NFXD = 1
*
END IF
*
- WORK( 1 ) = IWS
+ WORK( 1 ) = CMPLX( LWKOPT )
RETURN
*
* End of CGEQP3
*>
*> \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:
*>
--- /dev/null
+*> \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
*>
*> \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
--- /dev/null
+*> \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
RETURN
END IF
*
-* Quick return if possible.
-*
- IF( MINMN.EQ.0 ) THEN
- RETURN
- END IF
-*
* Move initial columns up front.
*
NFXD = 1
*>
*> \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:
*>
--- /dev/null
+*> \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
*>
*> \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
--- /dev/null
+*> \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
* ..
* .. Intrinsic Functions ..
INTRINSIC INT, MAX, MIN
-* ..
-* .. Executable Statements ..
+* Test input arguments
+* ====================
*
INFO = 0
LQUERY = ( LWORK.EQ.-1 )
RETURN
END IF
*
-* Quick return if possible.
-*
- IF( MINMN.EQ.0 ) THEN
- RETURN
- END IF
-*
* Move initial columns up front.
*
NFXD = 1
*>
*> \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:
*>
--- /dev/null
+*> \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
*>
*> \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
--- /dev/null
+*> \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
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
RETURN
END IF
*
-* Quick return if possible.
-*
- IF( MINMN.EQ.0 ) THEN
- RETURN
- END IF
-*
* Move initial columns up front.
*
NFXD = 1
*
END IF
*
- WORK( 1 ) = IWS
+ WORK( 1 ) = DCMPLX( LWKOPT )
RETURN
*
* End of ZGEQP3
*>
*> \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:
*>
--- /dev/null
+*> \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
*>
*> \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
--- /dev/null
+*> \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
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)
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
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)
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
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
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 \
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
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 \
*
* .. Parameters ..
INTEGER NTESTS
- PARAMETER ( NTESTS = 7 )
+ PARAMETER ( NTESTS = 12 )
INTEGER NTYPES
PARAMETER ( NTYPES = 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
CALL CGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
$ 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.
*>
*> 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:
* .. 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 ..
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
IFST = 1
ILST = 1
NT = 0
+ LWORK = 1
*
* Test error exits for the GG path.
*
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'
--- /dev/null
+*> \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
*
* .. Parameters ..
INTEGER NTESTS
- PARAMETER ( NTESTS = 7 )
+ PARAMETER ( NTESTS = 12 )
INTEGER NTYPES
PARAMETER ( NTYPES = 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
CALL DGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
$ 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.
*> \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:
* .. 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 ..
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
IFST = 1
ILST = 1
NT = 0
+ LWORK = 1
*
* Test error exits for the GG path.
*
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'
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'
--- /dev/null
+*> \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
*
* .. Parameters ..
INTEGER NTESTS
- PARAMETER ( NTESTS = 7 )
+ PARAMETER ( NTESTS = 12 )
INTEGER NTYPES
PARAMETER ( NTYPES = 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
CALL SGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
$ 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.
*>
*> 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:
* .. 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 ..
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
IFST = 1
ILST = 1
NT = 0
+ LWORK = 1
*
* Test error exits for the GG path.
*
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'
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'
--- /dev/null
+*> \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
*
* .. Parameters ..
INTEGER NTESTS
- PARAMETER ( NTESTS = 7 )
+ PARAMETER ( NTESTS = 12 )
INTEGER NTYPES
PARAMETER ( NTYPES = 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
CALL ZGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
$ 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.
*>
*> 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:
* .. 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 ..
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
IFST = 1
ILST = 1
NT = 0
+ LWORK = 1
*
* Test error exits for the GG path.
*
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'
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'
--- /dev/null
+*> \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
projects / platform / upstream / lapack.git / commitdiff