--- rpl/lapack/lapack/dtgsja.f 2011/07/22 07:38:12 1.8
+++ rpl/lapack/lapack/dtgsja.f 2011/11/21 20:43:06 1.9
@@ -1,11 +1,387 @@
+*> \brief \b DTGSJA
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DTGSJA + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE 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 )
+*
+* .. Scalar Arguments ..
+* CHARACTER JOBQ, JOBU, JOBV
+* INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
+* $ NCYCLE, P
+* DOUBLE PRECISION TOLA, TOLB
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ),
+* $ BETA( * ), Q( LDQ, * ), U( LDU, * ),
+* $ V( LDV, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DTGSJA computes the generalized singular value decomposition (GSVD)
+*> of two real upper triangular (or trapezoidal) matrices A and B.
+*>
+*> On entry, it is assumed that matrices A and B have the following
+*> forms, which may be obtained by the preprocessing subroutine DGGSVP
+*> from a general M-by-N matrix A and P-by-N matrix B:
+*>
+*> N-K-L K L
+*> A = K ( 0 A12 A13 ) if M-K-L >= 0;
+*> L ( 0 0 A23 )
+*> M-K-L ( 0 0 0 )
+*>
+*> N-K-L K L
+*> A = K ( 0 A12 A13 ) if M-K-L < 0;
+*> M-K ( 0 0 A23 )
+*>
+*> N-K-L K L
+*> B = 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.
+*>
+*> On exit,
+*>
+*> U**T *A*Q = D1*( 0 R ), V**T *B*Q = D2*( 0 R ),
+*>
+*> where U, V and Q are orthogonal matrices.
+*> R is a nonsingular upper triangular matrix, and D1 and D2 are
+*> ``diagonal'' matrices, which are of the following structures:
+*>
+*> 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 ) K
+*> L ( 0 0 R22 ) L
+*>
+*> 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.
+*>
+*> R = ( 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 computation of the orthogonal transformation matrices U, V or Q
+*> is optional. These matrices may either be formed explicitly, or they
+*> may be postmultiplied into input matrices U1, V1, or Q1.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] JOBU
+*> \verbatim
+*> JOBU is CHARACTER*1
+*> = 'U': U must contain an orthogonal matrix U1 on entry, and
+*> the product U1*U is returned;
+*> = 'I': U is initialized to the unit matrix, and the
+*> orthogonal matrix U is returned;
+*> = 'N': U is not computed.
+*> \endverbatim
+*>
+*> \param[in] JOBV
+*> \verbatim
+*> JOBV is CHARACTER*1
+*> = 'V': V must contain an orthogonal matrix V1 on entry, and
+*> the product V1*V is returned;
+*> = 'I': V is initialized to the unit matrix, and the
+*> orthogonal matrix V is returned;
+*> = 'N': V is not computed.
+*> \endverbatim
+*>
+*> \param[in] JOBQ
+*> \verbatim
+*> JOBQ is CHARACTER*1
+*> = 'Q': Q must contain an orthogonal matrix Q1 on entry, and
+*> the product Q1*Q is returned;
+*> = 'I': Q is initialized to the unit matrix, and the
+*> orthogonal matrix Q is returned;
+*> = '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] K
+*> \verbatim
+*> K is INTEGER
+*> \endverbatim
+*>
+*> \param[in] L
+*> \verbatim
+*> L is INTEGER
+*>
+*> K and L specify the subblocks in the input matrices A and B:
+*> A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N)
+*> of A and B, whose GSVD is going to be computed by DTGSJA.
+*> See Further Details.
+*> \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(N-K+1:N,1:MIN(K+L,M) ) 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, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
+*> a part of R. 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[in] TOLA
+*> \verbatim
+*> TOLA is DOUBLE PRECISION
+*> \endverbatim
+*>
+*> \param[in] TOLB
+*> \verbatim
+*> TOLB is DOUBLE PRECISION
+*>
+*> TOLA and TOLB are the convergence criteria for the Jacobi-
+*> Kogbetliantz iteration procedure. Generally, they are the
+*> same as used in the preprocessing step, say
+*> TOLA = max(M,N)*norm(A)*MAZHEPS,
+*> TOLB = max(P,N)*norm(B)*MAZHEPS.
+*> \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) = diag(C),
+*> BETA(K+1:K+L) = diag(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.
+*> Furthermore, if K+L < N,
+*> ALPHA(K+L+1:N) = 0 and
+*> BETA(K+L+1:N) = 0.
+*> \endverbatim
+*>
+*> \param[in,out] U
+*> \verbatim
+*> U is DOUBLE PRECISION array, dimension (LDU,M)
+*> On entry, if JOBU = 'U', U must contain a matrix U1 (usually
+*> the orthogonal matrix returned by DGGSVP).
+*> On exit,
+*> if JOBU = 'I', U contains the orthogonal matrix U;
+*> if JOBU = 'U', U contains the product U1*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[in,out] V
+*> \verbatim
+*> V is DOUBLE PRECISION array, dimension (LDV,P)
+*> On entry, if JOBV = 'V', V must contain a matrix V1 (usually
+*> the orthogonal matrix returned by DGGSVP).
+*> On exit,
+*> if JOBV = 'I', V contains the orthogonal matrix V;
+*> if JOBV = 'V', V contains the product V1*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[in,out] Q
+*> \verbatim
+*> Q is DOUBLE PRECISION array, dimension (LDQ,N)
+*> On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
+*> the orthogonal matrix returned by DGGSVP).
+*> On exit,
+*> if JOBQ = 'I', Q contains the orthogonal matrix Q;
+*> if JOBQ = 'Q', Q contains the product Q1*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 (2*N)
+*> \endverbatim
+*>
+*> \param[out] NCYCLE
+*> \verbatim
+*> NCYCLE is INTEGER
+*> The number of cycles required for convergence.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value.
+*> = 1: the procedure does not converge after MAXIT cycles.
+*> \endverbatim
+*>
+*> \verbatim
+*> Internal Parameters
+*> ===================
+*>
+*> MAXIT INTEGER
+*> MAXIT specifies the total loops that the iterative procedure
+*> may take. If after MAXIT cycles, the routine fails to
+*> converge, we return INFO = 1.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup doubleOTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> DTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
+*> min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
+*> matrix B13 to the form:
+*>
+*> U1**T *A13*Q1 = C1*R1; V1**T *B13*Q1 = S1*R1,
+*>
+*> where U1, V1 and Q1 are orthogonal matrix, and Z**T is the transpose
+*> of Z. C1 and S1 are diagonal matrices satisfying
+*>
+*> C1**2 + S1**2 = I,
+*>
+*> and R1 is an L-by-L nonsingular upper triangular matrix.
+*> \endverbatim
+*>
+* =====================================================================
SUBROUTINE 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 )
*
-* -- LAPACK routine (version 3.3.1) --
+* -- LAPACK computational routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* -- April 2009 --
+* November 2011
*
* .. Scalar Arguments ..
CHARACTER JOBQ, JOBU, JOBV
@@ -19,243 +395,6 @@
$ V( LDV, * ), WORK( * )
* ..
*
-* Purpose
-* =======
-*
-* DTGSJA computes the generalized singular value decomposition (GSVD)
-* of two real upper triangular (or trapezoidal) matrices A and B.
-*
-* On entry, it is assumed that matrices A and B have the following
-* forms, which may be obtained by the preprocessing subroutine DGGSVP
-* from a general M-by-N matrix A and P-by-N matrix B:
-*
-* N-K-L K L
-* A = K ( 0 A12 A13 ) if M-K-L >= 0;
-* L ( 0 0 A23 )
-* M-K-L ( 0 0 0 )
-*
-* N-K-L K L
-* A = K ( 0 A12 A13 ) if M-K-L < 0;
-* M-K ( 0 0 A23 )
-*
-* N-K-L K L
-* B = 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.
-*
-* On exit,
-*
-* U**T *A*Q = D1*( 0 R ), V**T *B*Q = D2*( 0 R ),
-*
-* where U, V and Q are orthogonal matrices.
-* R is a nonsingular upper triangular matrix, and D1 and D2 are
-* ``diagonal'' matrices, which are of the following structures:
-*
-* 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 ) K
-* L ( 0 0 R22 ) L
-*
-* 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.
-*
-* R = ( 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 computation of the orthogonal transformation matrices U, V or Q
-* is optional. These matrices may either be formed explicitly, or they
-* may be postmultiplied into input matrices U1, V1, or Q1.
-*
-* Arguments
-* =========
-*
-* JOBU (input) CHARACTER*1
-* = 'U': U must contain an orthogonal matrix U1 on entry, and
-* the product U1*U is returned;
-* = 'I': U is initialized to the unit matrix, and the
-* orthogonal matrix U is returned;
-* = 'N': U is not computed.
-*
-* JOBV (input) CHARACTER*1
-* = 'V': V must contain an orthogonal matrix V1 on entry, and
-* the product V1*V is returned;
-* = 'I': V is initialized to the unit matrix, and the
-* orthogonal matrix V is returned;
-* = 'N': V is not computed.
-*
-* JOBQ (input) CHARACTER*1
-* = 'Q': Q must contain an orthogonal matrix Q1 on entry, and
-* the product Q1*Q is returned;
-* = 'I': Q is initialized to the unit matrix, and the
-* orthogonal matrix Q is returned;
-* = 'N': Q is not computed.
-*
-* M (input) INTEGER
-* The number of rows of the matrix A. M >= 0.
-*
-* P (input) INTEGER
-* The number of rows of the matrix B. P >= 0.
-*
-* N (input) INTEGER
-* The number of columns of the matrices A and B. N >= 0.
-*
-* K (input) INTEGER
-* L (input) INTEGER
-* K and L specify the subblocks in the input matrices A and B:
-* A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N)
-* of A and B, whose GSVD is going to be computed by DTGSJA.
-* See Further Details.
-*
-* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-* On entry, the M-by-N matrix A.
-* On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
-* matrix R or part of R. See Purpose for details.
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,M).
-*
-* B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
-* On entry, the P-by-N matrix B.
-* On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
-* a part of R. See Purpose for details.
-*
-* LDB (input) INTEGER
-* The leading dimension of the array B. LDB >= max(1,P).
-*
-* TOLA (input) DOUBLE PRECISION
-* TOLB (input) DOUBLE PRECISION
-* TOLA and TOLB are the convergence criteria for the Jacobi-
-* Kogbetliantz iteration procedure. Generally, they are the
-* same as used in the preprocessing step, say
-* TOLA = max(M,N)*norm(A)*MAZHEPS,
-* TOLB = max(P,N)*norm(B)*MAZHEPS.
-*
-* ALPHA (output) DOUBLE PRECISION array, dimension (N)
-* BETA (output) 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) = diag(C),
-* BETA(K+1:K+L) = diag(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.
-* Furthermore, if K+L < N,
-* ALPHA(K+L+1:N) = 0 and
-* BETA(K+L+1:N) = 0.
-*
-* U (input/output) DOUBLE PRECISION array, dimension (LDU,M)
-* On entry, if JOBU = 'U', U must contain a matrix U1 (usually
-* the orthogonal matrix returned by DGGSVP).
-* On exit,
-* if JOBU = 'I', U contains the orthogonal matrix U;
-* if JOBU = 'U', U contains the product U1*U.
-* If JOBU = 'N', U is not referenced.
-*
-* LDU (input) INTEGER
-* The leading dimension of the array U. LDU >= max(1,M) if
-* JOBU = 'U'; LDU >= 1 otherwise.
-*
-* V (input/output) DOUBLE PRECISION array, dimension (LDV,P)
-* On entry, if JOBV = 'V', V must contain a matrix V1 (usually
-* the orthogonal matrix returned by DGGSVP).
-* On exit,
-* if JOBV = 'I', V contains the orthogonal matrix V;
-* if JOBV = 'V', V contains the product V1*V.
-* If JOBV = 'N', V is not referenced.
-*
-* LDV (input) INTEGER
-* The leading dimension of the array V. LDV >= max(1,P) if
-* JOBV = 'V'; LDV >= 1 otherwise.
-*
-* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
-* On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
-* the orthogonal matrix returned by DGGSVP).
-* On exit,
-* if JOBQ = 'I', Q contains the orthogonal matrix Q;
-* if JOBQ = 'Q', Q contains the product Q1*Q.
-* If JOBQ = 'N', Q is not referenced.
-*
-* LDQ (input) INTEGER
-* The leading dimension of the array Q. LDQ >= max(1,N) if
-* JOBQ = 'Q'; LDQ >= 1 otherwise.
-*
-* WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
-*
-* NCYCLE (output) INTEGER
-* The number of cycles required for convergence.
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: if INFO = -i, the i-th argument had an illegal value.
-* = 1: the procedure does not converge after MAXIT cycles.
-*
-* Internal Parameters
-* ===================
-*
-* MAXIT INTEGER
-* MAXIT specifies the total loops that the iterative procedure
-* may take. If after MAXIT cycles, the routine fails to
-* converge, we return INFO = 1.
-*
-* Further Details
-* ===============
-*
-* DTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
-* min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
-* matrix B13 to the form:
-*
-* U1**T *A13*Q1 = C1*R1; V1**T *B13*Q1 = S1*R1,
-*
-* where U1, V1 and Q1 are orthogonal matrix, and Z**T is the transpose
-* of Z. C1 and S1 are diagonal matrices satisfying
-*
-* C1**2 + S1**2 = I,
-*
-* and R1 is an L-by-L nonsingular upper triangular matrix.
-*
* =====================================================================
*
* .. Parameters ..