Modules
module	la_lapack_svd_comp
	SVD components: bidiagonal reduction and its orthogonal factors, Jacobi sweeps, generalized SVD.

Functions/Subroutines
pure subroutine, public	la_lapack_svd_comp::la_sgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, info)
	SGBBRD: reduces a real general m-by-n band matrix A to upper bidiagonal form B by an orthogonal transformation: Q*T A * P = B. The routine computes B, and optionally forms Q or PT, or computes QT*C for a given matrix C.

pure subroutine, public	la_lapack_svd_comp::la_dgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, info)
	DGBBRD: reduces a real general m-by-n band matrix A to upper bidiagonal form B by an orthogonal transformation: Q*T A * P = B. The routine computes B, and optionally forms Q or PT, or computes QT*C for a given matrix C.

pure subroutine, public	la_lapack_svd_comp::la_qgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, info)
	QGBBRD: reduces a real general m-by-n band matrix A to upper bidiagonal form B by an orthogonal transformation: Q*T A * P = B. The routine computes B, and optionally forms Q or PT, or computes QT*C for a given matrix C.

pure subroutine, public	la_lapack_svd_comp::la_sgebd2 (m, n, a, lda, d, e, tauq, taup, work, info)
	SGEBD2: reduces a real general m by n matrix A to upper or lower bidiagonal form B by an orthogonal transformation: Q*T A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_dgebd2 (m, n, a, lda, d, e, tauq, taup, work, info)
	DGEBD2: reduces a real general m by n matrix A to upper or lower bidiagonal form B by an orthogonal transformation: Q*T A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_qgebd2 (m, n, a, lda, d, e, tauq, taup, work, info)
	QGEBD2: reduces a real general m by n matrix A to upper or lower bidiagonal form B by an orthogonal transformation: Q*T A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_sgsvj0 (jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	SGSVJ0: is called from SGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as SGESVJ does, but it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer.

pure subroutine, public	la_lapack_svd_comp::la_dgsvj0 (jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	DGSVJ0: is called from DGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as DGESVJ does, but it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer.

pure subroutine, public	la_lapack_svd_comp::la_qgsvj0 (jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	QGSVJ0: is called from QGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as QGESVJ does, but it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer.

pure subroutine, public	la_lapack_svd_comp::la_sgsvj1 (jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	SGSVJ1: is called from SGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as SGESVJ does, but it targets only particular pivots and it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer. Further Details.

pure subroutine, public	la_lapack_svd_comp::la_dgsvj1 (jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	DGSVJ1: is called from DGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as DGESVJ does, but it targets only particular pivots and it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer. Further Details.

pure subroutine, public	la_lapack_svd_comp::la_qgsvj1 (jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	QGSVJ1: is called from QGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as QGESVJ does, but it targets only particular pivots and it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer. Further Details.

pure subroutine, public	la_lapack_svd_comp::la_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)
	STGSJA: 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 SGGSVP 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 AQ = D1( 0 R ), V*T BQ = 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) ), C2 + S2 = 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) ), C2 + S2 = 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.

pure subroutine, public	la_lapack_svd_comp::la_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)
	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 AQ = D1( 0 R ), V*T BQ = 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) ), C2 + S2 = 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) ), C2 + S2 = 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.

pure subroutine, public	la_lapack_svd_comp::la_qtgsja (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)
	QTGSJA: 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 AQ = D1( 0 R ), V*T BQ = 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) ), C2 + S2 = 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) ), C2 + S2 = 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.

pure subroutine, public	la_lapack_svd_comp::la_sgebrd (m, n, a, lda, d, e, tauq, taup, work, lwork, info)
	SGEBRD: reduces a general real M-by-N matrix A to upper or lower bidiagonal form B by an orthogonal transformation: Q*T A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_dgebrd (m, n, a, lda, d, e, tauq, taup, work, lwork, info)
	DGEBRD: reduces a general real M-by-N matrix A to upper or lower bidiagonal form B by an orthogonal transformation: Q*T A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_qgebrd (m, n, a, lda, d, e, tauq, taup, work, lwork, info)
	QGEBRD: reduces a general real M-by-N matrix A to upper or lower bidiagonal form B by an orthogonal transformation: Q*T A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_sorgbr (vect, m, n, k, a, lda, tau, work, lwork, info)
	SORGBR: generates one of the real orthogonal matrices Q or P*T determined by SGEBRD when reducing a real matrix A to bidiagonal form: A = Q B * PT. Q and PT are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and SORGBR returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and SORGBR returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and PT is of order N: if k < n, PT = G(k) . . . G(2) G(1) and SORGBR returns the first m rows of PT, where n >= m >= k; if k >= n, PT = G(n-1) . . . G(2) G(1) and SORGBR returns P**T as an N-by-N matrix.

pure subroutine, public	la_lapack_svd_comp::la_dorgbr (vect, m, n, k, a, lda, tau, work, lwork, info)
	DORGBR: generates one of the real orthogonal matrices Q or P*T determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q B * PT. Q and PT are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and PT is of order N: if k < n, PT = G(k) . . . G(2) G(1) and DORGBR returns the first m rows of PT, where n >= m >= k; if k >= n, PT = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as an N-by-N matrix.

pure subroutine, public	la_lapack_svd_comp::la_qorgbr (vect, m, n, k, a, lda, tau, work, lwork, info)
	QORGBR: generates one of the real orthogonal matrices Q or P*T determined by QGEBRD when reducing a real matrix A to bidiagonal form: A = Q B * PT. Q and PT are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and QORGBR returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and QORGBR returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and PT is of order N: if k < n, PT = G(k) . . . G(2) G(1) and QORGBR returns the first m rows of PT, where n >= m >= k; if k >= n, PT = G(n-1) . . . G(2) G(1) and QORGBR returns P**T as an N-by-N matrix.

pure subroutine, public	la_lapack_svd_comp::la_sormbr (vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
	If VECT = 'Q', SORMBR: overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q*T C C * Q*T If VECT = 'P', SORMBR overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P C C * P TRANS = 'T': P*T C C * PT Here Q and PT are the orthogonal matrices determined by SGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * PT. Q and PT are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the orthogonal matrix Q or P**T that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1).

pure subroutine, public	la_lapack_svd_comp::la_dormbr (vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
	If VECT = 'Q', DORMBR: overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q*T C C * Q*T If VECT = 'P', DORMBR overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P C C * P TRANS = 'T': P*T C C * PT Here Q and PT are the orthogonal matrices determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * PT. Q and PT are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the orthogonal matrix Q or P**T that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1).

pure subroutine, public	la_lapack_svd_comp::la_qormbr (vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
	If VECT = 'Q', QORMBR: overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'T': Q*T C C * Q*T If VECT = 'P', QORMBR overwrites the general real M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P C C * P TRANS = 'T': P*T C C * PT Here Q and PT are the orthogonal matrices determined by QGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * PT. Q and PT are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the orthogonal matrix Q or P**T that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1).

pure subroutine, public	la_lapack_svd_comp::la_cgebd2 (m, n, a, lda, d, e, tauq, taup, work, info)
	CGEBD2: reduces a complex general m by n matrix A to upper or lower real bidiagonal form B by a unitary transformation: Q*H A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_zgebd2 (m, n, a, lda, d, e, tauq, taup, work, info)
	ZGEBD2: reduces a complex general m by n matrix A to upper or lower real bidiagonal form B by a unitary transformation: Q*H A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_wgebd2 (m, n, a, lda, d, e, tauq, taup, work, info)
	WGEBD2: reduces a complex general m by n matrix A to upper or lower real bidiagonal form B by a unitary transformation: Q*H A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_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)
	CTGSJA: computes the generalized singular value decomposition (GSVD) of two complex 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 CGGSVP 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*H AQ = D1( 0 R ), V*H BQ = D2( 0 R ), where U, V and Q are unitary 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) ), C2 + S2 = 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) ), C2 + S2 = 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 unitary 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.

pure subroutine, public	la_lapack_svd_comp::la_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)
	ZTGSJA: computes the generalized singular value decomposition (GSVD) of two complex 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 ZGGSVP 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*H AQ = D1( 0 R ), V*H BQ = D2( 0 R ), where U, V and Q are unitary 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) ), C2 + S2 = 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) ), C2 + S2 = 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 unitary 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.

pure subroutine, public	la_lapack_svd_comp::la_wtgsja (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)
	WTGSJA: computes the generalized singular value decomposition (GSVD) of two complex 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 ZGGSVP 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*H AQ = D1( 0 R ), V*H BQ = D2( 0 R ), where U, V and Q are unitary 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) ), C2 + S2 = 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) ), C2 + S2 = 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 unitary 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.

pure subroutine, public	la_lapack_svd_comp::la_cgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, rwork, info)
	CGBBRD: reduces a complex general m-by-n band matrix A to real upper bidiagonal form B by a unitary transformation: Q*H A * P = B. The routine computes B, and optionally forms Q or PH, or computes QH*C for a given matrix C.

pure subroutine, public	la_lapack_svd_comp::la_zgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, rwork, info)
	ZGBBRD: reduces a complex general m-by-n band matrix A to real upper bidiagonal form B by a unitary transformation: Q*H A * P = B. The routine computes B, and optionally forms Q or PH, or computes QH*C for a given matrix C.

pure subroutine, public	la_lapack_svd_comp::la_wgbbrd (vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, rwork, info)
	WGBBRD: reduces a complex general m-by-n band matrix A to real upper bidiagonal form B by a unitary transformation: Q*H A * P = B. The routine computes B, and optionally forms Q or PH, or computes QH*C for a given matrix C.

pure subroutine, public	la_lapack_svd_comp::la_cgebrd (m, n, a, lda, d, e, tauq, taup, work, lwork, info)
	CGEBRD: reduces a general complex M-by-N matrix A to upper or lower bidiagonal form B by a unitary transformation: Q*H A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_zgebrd (m, n, a, lda, d, e, tauq, taup, work, lwork, info)
	ZGEBRD: reduces a general complex M-by-N matrix A to upper or lower bidiagonal form B by a unitary transformation: Q*H A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_wgebrd (m, n, a, lda, d, e, tauq, taup, work, lwork, info)
	WGEBRD: reduces a general complex M-by-N matrix A to upper or lower bidiagonal form B by a unitary transformation: Q*H A * P = B. If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.

pure subroutine, public	la_lapack_svd_comp::la_cungbr (vect, m, n, k, a, lda, tau, work, lwork, info)
	CUNGBR: generates one of the complex unitary matrices Q or P*H determined by CGEBRD when reducing a complex matrix A to bidiagonal form: A = Q B * PH. Q and PH are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and CUNGBR returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and CUNGBR returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and PH is of order N: if k < n, PH = G(k) . . . G(2) G(1) and CUNGBR returns the first m rows of PH, where n >= m >= k; if k >= n, PH = G(n-1) . . . G(2) G(1) and CUNGBR returns P**H as an N-by-N matrix.

pure subroutine, public	la_lapack_svd_comp::la_zungbr (vect, m, n, k, a, lda, tau, work, lwork, info)
	ZUNGBR: generates one of the complex unitary matrices Q or P*H determined by ZGEBRD when reducing a complex matrix A to bidiagonal form: A = Q B * PH. Q and PH are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and ZUNGBR returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and ZUNGBR returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and PH is of order N: if k < n, PH = G(k) . . . G(2) G(1) and ZUNGBR returns the first m rows of PH, where n >= m >= k; if k >= n, PH = G(n-1) . . . G(2) G(1) and ZUNGBR returns P**H as an N-by-N matrix.

pure subroutine, public	la_lapack_svd_comp::la_wungbr (vect, m, n, k, a, lda, tau, work, lwork, info)
	WUNGBR: generates one of the complex unitary matrices Q or P*H determined by WGEBRD when reducing a complex matrix A to bidiagonal form: A = Q B * PH. Q and PH are defined as products of elementary reflectors H(i) or G(i) respectively. If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q is of order M: if m >= k, Q = H(1) H(2) . . . H(k) and WUNGBR returns the first n columns of Q, where m >= n >= k; if m < k, Q = H(1) H(2) . . . H(m-1) and WUNGBR returns Q as an M-by-M matrix. If VECT = 'P', A is assumed to have been a K-by-N matrix, and PH is of order N: if k < n, PH = G(k) . . . G(2) G(1) and WUNGBR returns the first m rows of PH, where n >= m >= k; if k >= n, PH = G(n-1) . . . G(2) G(1) and WUNGBR returns P**H as an N-by-N matrix.

pure subroutine, public	la_lapack_svd_comp::la_cunmbr (vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
	If VECT = 'Q', CUNMBR: overwrites the general complex M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'C': Q*H C C * Q*H If VECT = 'P', CUNMBR overwrites the general complex M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P C C * P TRANS = 'C': P*H C C * PH Here Q and PH are the unitary matrices determined by CGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * PH. Q and PH are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the unitary matrix Q or P**H that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1).

pure subroutine, public	la_lapack_svd_comp::la_zunmbr (vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
	If VECT = 'Q', ZUNMBR: overwrites the general complex M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'C': Q*H C C * Q*H If VECT = 'P', ZUNMBR overwrites the general complex M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P C C * P TRANS = 'C': P*H C C * PH Here Q and PH are the unitary matrices determined by ZGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * PH. Q and PH are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the unitary matrix Q or P**H that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1).

pure subroutine, public	la_lapack_svd_comp::la_wunmbr (vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
	If VECT = 'Q', WUNMBR: overwrites the general complex M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': Q * C C * Q TRANS = 'C': Q*H C C * Q*H If VECT = 'P', WUNMBR overwrites the general complex M-by-N matrix C with SIDE = 'L' SIDE = 'R' TRANS = 'N': P C C * P TRANS = 'C': P*H C C * PH Here Q and PH are the unitary matrices determined by WGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * PH. Q and PH are defined as products of elementary reflectors H(i) and G(i) respectively. Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the order of the unitary matrix Q or P**H that is applied. If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: if nq >= k, Q = H(1) H(2) . . . H(k); if nq < k, Q = H(1) H(2) . . . H(nq-1). If VECT = 'P', A is assumed to have been a K-by-NQ matrix: if k < nq, P = G(1) G(2) . . . G(k); if k >= nq, P = G(1) G(2) . . . G(nq-1).

pure subroutine, public	la_lapack_svd_comp::la_cgsvj0 (jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	CGSVJ0: is called from CGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as CGESVJ does, but it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer.

pure subroutine, public	la_lapack_svd_comp::la_zgsvj0 (jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	ZGSVJ0: is called from ZGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer.

pure subroutine, public	la_lapack_svd_comp::la_wgsvj0 (jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	WGSVJ0: is called from WGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as WGESVJ does, but it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer.

pure subroutine, public	la_lapack_svd_comp::la_cgsvj1 (jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	CGSVJ1: is called from CGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as CGESVJ does, but it targets only particular pivots and it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer. Further Details.

pure subroutine, public	la_lapack_svd_comp::la_zgsvj1 (jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	ZGSVJ1: is called from ZGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but it targets only particular pivots and it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer. Further Details.

pure subroutine, public	la_lapack_svd_comp::la_wgsvj1 (jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
	WGSVJ1: is called from WGESVJ as a pre-processor and that is its main purpose. It applies Jacobi rotations in the same way as WGESVJ does, but it targets only particular pivots and it does not check convergence (stopping criterion). Few tuning parameters (marked by [TP]) are available for the implementer. Further Details.

Modules

Functions/Subroutines