! \brief DGGSVD3 computes the singular value decomposition (SVD) for OTHER matrices
!
! =========== DOCUMENTATION ===========
!
! Online html documentation available at
! http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DGGSVD3 + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \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 doubleGEsing
!
!> \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.7.0) --
! -- LAPACK is a software package provided by Univ. of Tennessee, --
! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! August 2015
program main
implicit none
integer, parameter :: dp = kind(0.d0)
INTEGER, PARAMETER :: M = 6, N = 4, P = 3, NN = 20, MM= 8
INTEGER, PARAMETER :: N1 = 6, N2 = 4, N3 = 3, N4 = 2
INTEGER, PARAMETER :: N5 = 2, M5 = 3, N6 = 3, M6 = 4
CHARACTER :: JOBQ = 'N', JOBU= 'N', JOBV= 'N'
! CHARACTER :: JOBQ = 'N', JOBU= 'A', JOBV= 'A'
! CHARACTER :: JOBQ = 'Q', JOBU = 'U', JOBV = 'V'
INTEGER :: INFO, K, L, KK, I, J
INTEGER :: LWORK = 40 * M
!LDA = M, LDB = P, LDU = M, LDV = P , LDQ = N,
INTEGER, ALLOCATABLE, DIMENSION ( : ) :: IWORK
REAL ( dp ), ALLOCATABLE, DIMENSION ( : ) :: ALPHA, BETA, WORK
Real (dp), ALLOCATABLE, DIMENSION ( : , : ) :: A, B
Real (dp), ALLOCATABLE, DIMENSION ( : , : ) :: U, V, Q
! Real (dp), DIMENSION ( M, N ) :: A
! Real (dp), DIMENSION ( P, N ) :: B
Real (dp), DIMENSION ( NN ) :: AA
Real (dp), DIMENSION ( MM ) :: BB
Real (dp), DIMENSION ( N1 ) :: A1
Real (dp), DIMENSION ( N2 ) :: A2
Real (dp), DIMENSION ( N3 ) :: A3
Real (dp), DIMENSION ( N4 ) :: A4
Real (dp), DIMENSION ( N5, M5 ) :: A5
Real (dp), DIMENSION ( N6, M6 ) :: A6
ALLOCATE ( IWORK ( N ) )
ALLOCATE ( ALPHA( N ), BETA( N ), WORK( LWORK ) )
! ALLOCATE ( Q( LDQ, N ), U( LDU, M ), V( LDV, P ) )
! ALLOCATE ( A( LDA, N ), B( LDB, N ) )
ALLOCATE ( Q( N , N ) )
ALLOCATE ( U( M , M ), V( P , P ) )
ALLOCATE ( A( M , N ) )
ALLOCATE ( B( P , N ) )
! INFO = 0
! IWORK = 0
! ALPHA = 0.0D00
! BETA = 0.0D00
! WORK = 0.0D00
! Q = 0.0D00
! U = 0.0D00
! V = 0.0D00
! A = 0.0D00
! B = 0.0D00
A = RESHAPE( (/ 1.0D00, 2.0D00 , 3.0D00, 4.0D00, 2.5D00, 3.6D00, 4.18D00,&
5.45D00 , 6.78D00 , 8.12D00 , 9.34D00 , 10.12D00 , 11.0D00, 12.0D00 , 34.0D00, -24.0D00, 42.5D00, 13.6D00, -41.18D00,&
46.45D00 , 16.12D00 , -14.2D00 , -.034D00 , 0.12D00 /) , (/ M , N /))
! A = RESHAPE( (/ 11.0D00, 2.0D00 , 3.0D00, 4.0D00, 2.5D00, 3.6D00, 4.18D00,&
! 5.45D00 , 6.78D00 , 8.12D00 , 9.34D00 , 10.12D00 , 11.0D00, 12.0D00 , 34.0D00, -24.0D00 /) , (/ LDA , N /))
B = RESHAPE( (/ 11.0D00, -3.0D00 , 0.0D00, 0.0D00, 0.0D00, -71.6D00, 4.98D00,&
0.00D00 , 0.0D00 , 0.0D00 , -9.34D00 , 10.12D00 /) , (/ P , N /))
write (*,*) "Hello world !! "
write (*,*) 1.0_dp/3
write (*,*) 1.0D00/3
write (*,*) 1.0D00/3.0D00
write (*,*) "INFO = ", INFO
write (*,*) " check 000 "
write (*,*) " dp = ", dp
write (*,*) " double = ", 1._8
! call dggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, LWORK, IWORK, INFO )
! call dggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, LDA, LDB, ALPHA, BETA, LDU, LDV, LDQ, WORK, LWORK, IWORK, INFO )
! call DGESVD( JOBU, JOBV, M, N, A, LDA, ALPHA, U, LDU, V, LDV, WORK, LWORK, INFO )
! call 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 )
call gen_svd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, B, ALPHA, BETA, U, V, Q, WORK, LWORK, IWORK, INFO )
write (*,*) "K , L = ", K, L
write (*,*) "------ E N D ----------"
A3 = 0.0_dp
A4 = 0.0_dp
A5 = 1.0_dp
! call test_3 (A1, A2, A3, A4, N1, N2, N3, N4, N5 )
! write (*,*) " N5 = ", N5
! call test_4 (A1, A2, A3, A4, A5, A6, N1, N2, N3, N4, KK, N5, M5, N6, M6)
! write (*,*) " KK = ", KK
! BB = 0.0_dp
! call test_2 (AA, BB, NN, MM)
write (*,*) "INFO = ", INFO
write (*,*) "K , L = ", K, L
WRITE (*,*) "---------------------------------------------------------------------"
DO I = 1 , N
write (*,*) " I, ALPHA (I), BETA(I) = ", I, ALPHA (I), BETA(I)
END DO
WRITE (*,*) "---------------------------------------------------------------------"
DO J = 1 , M
DO I = 1 , M
write (*,*) " J, I, U(J,I) = ", J, I, U(J,I)
END DO
END DO
WRITE (*,*) "---------------------------------------------------------------------"
DO J = 1 , P
DO I = 1 , P
write (*,*) " J, I, V(J,I) = ", J, I, V(J,I)
END DO
END DO
WRITE (*,*) "---------------------------------------------------------------------"
DO J = 1 , N
DO I = 1 , N
write (*,*) " J, I, Q(J,I) = ", J, I, Q(J,I)
END DO
END DO
write (*,*) "--------------- E N D ------------------------"
! DEALLOCATE ( A, B, U, V, Q )
! DEALLOCATE ( ALPHA, BETA, WORK, IWORK )
end program main