SVD components: bidiagonal reduction and its orthogonal factors, Jacobi sweeps, generalized SVD. More...
Functions/Subroutines | |
| pure subroutine, public | 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 P**T, or computes Q**T*C for a given matrix C. | |
| pure subroutine, public | 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 P**T, or computes Q**T*C for a given matrix C. | |
| pure subroutine, public | 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 P**T, or computes Q**T*C for a given matrix C. | |
| pure subroutine, public | 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_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_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_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_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_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_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_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_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_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 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. | |
| pure subroutine, public | 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 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. | |
| pure subroutine, public | 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 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. | |
| pure subroutine, public | 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_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_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_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 * P**T. Q and P**T 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 P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and SORGBR returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and SORGBR returns P**T as an N-by-N matrix. | |
| pure subroutine, public | 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 * P**T. Q and P**T 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 P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as an N-by-N matrix. | |
| pure subroutine, public | 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 * P**T. Q and P**T 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 P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and QORGBR returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and QORGBR returns P**T as an N-by-N matrix. | |
| pure subroutine, public | 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 * P**T Here Q and P**T are the orthogonal matrices determined by SGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T 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_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 * P**T Here Q and P**T are the orthogonal matrices determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T 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_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 * P**T Here Q and P**T are the orthogonal matrices determined by QGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T 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_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_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_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_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 A*Q = D1( 0 R ), V**H B*Q = 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) ), 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 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_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 A*Q = D1( 0 R ), V**H B*Q = 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) ), 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 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_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 A*Q = D1( 0 R ), V**H B*Q = 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) ), 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 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_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 P**H, or computes Q**H*C for a given matrix C. | |
| pure subroutine, public | 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 P**H, or computes Q**H*C for a given matrix C. | |
| pure subroutine, public | 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 P**H, or computes Q**H*C for a given matrix C. | |
| pure subroutine, public | 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_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_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_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 * P**H. Q and P**H 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 P**H is of order N: if k < n, P**H = G(k) . . . G(2) G(1) and CUNGBR returns the first m rows of P**H, where n >= m >= k; if k >= n, P**H = G(n-1) . . . G(2) G(1) and CUNGBR returns P**H as an N-by-N matrix. | |
| pure subroutine, public | 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 * P**H. Q and P**H 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 P**H is of order N: if k < n, P**H = G(k) . . . G(2) G(1) and ZUNGBR returns the first m rows of P**H, where n >= m >= k; if k >= n, P**H = G(n-1) . . . G(2) G(1) and ZUNGBR returns P**H as an N-by-N matrix. | |
| pure subroutine, public | 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 * P**H. Q and P**H 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 P**H is of order N: if k < n, P**H = G(k) . . . G(2) G(1) and WUNGBR returns the first m rows of P**H, where n >= m >= k; if k >= n, P**H = G(n-1) . . . G(2) G(1) and WUNGBR returns P**H as an N-by-N matrix. | |
| pure subroutine, public | 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 * P**H Here Q and P**H are the unitary matrices determined by CGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q and P**H 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_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 * P**H Here Q and P**H are the unitary matrices determined by ZGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q and P**H 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_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 * P**H Here Q and P**H are the unitary matrices determined by WGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q and P**H 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_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_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_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_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_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_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. | |
Detailed Description
SVD components: bidiagonal reduction and its orthogonal factors, Jacobi sweeps, generalized SVD.
Function/Subroutine Documentation
◆ la_cgbbrd()
| pure subroutine, public la_lapack_svd_comp::la_cgbbrd | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | ncc, | ||
| integer(ilp), intent(in) | kl, | ||
| integer(ilp), intent(in) | ku, | ||
| complex(sp), dimension(ldab,*), intent(inout) | ab, | ||
| integer(ilp), intent(in) | ldab, | ||
| real(sp), dimension(*), intent(out) | d, | ||
| real(sp), dimension(*), intent(out) | e, | ||
| complex(sp), dimension(ldq,*), intent(out) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| complex(sp), dimension(ldpt,*), intent(out) | pt, | ||
| integer(ilp), intent(in) | ldpt, | ||
| complex(sp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| complex(sp), dimension(*), intent(out) | work, | ||
| real(sp), dimension(*), intent(out) | rwork, | ||
| integer(ilp), intent(out) | 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 P**H, or computes Q**H*C for a given matrix C.
◆ la_cgebd2()
| pure subroutine, public la_lapack_svd_comp::la_cgebd2 | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(*), intent(out) | d, | ||
| real(sp), dimension(*), intent(out) | e, | ||
| complex(sp), dimension(*), intent(out) | tauq, | ||
| complex(sp), dimension(*), intent(out) | taup, | ||
| complex(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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.
◆ la_cgebrd()
| pure subroutine, public la_lapack_svd_comp::la_cgebrd | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(*), intent(out) | d, | ||
| real(sp), dimension(*), intent(out) | e, | ||
| complex(sp), dimension(*), intent(out) | tauq, | ||
| complex(sp), dimension(*), intent(out) | taup, | ||
| complex(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_cgsvj0()
| pure subroutine, public la_lapack_svd_comp::la_cgsvj0 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(sp), dimension(n), intent(inout) | d, | ||
| real(sp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| complex(sp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(sp), intent(in) | eps, | ||
| real(sp), intent(in) | sfmin, | ||
| real(sp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| complex(sp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_cgsvj1()
| pure subroutine, public la_lapack_svd_comp::la_cgsvj1 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | n1, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(sp), dimension(n), intent(inout) | d, | ||
| real(sp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| complex(sp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(sp), intent(in) | eps, | ||
| real(sp), intent(in) | sfmin, | ||
| real(sp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| complex(sp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
CGSVJ1 applies few sweeps of Jacobi rotations in the column space of the input M-by-N matrix A. The pivot pairs are taken from the (1,2) off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The block-entries (tiles) of the (1,2) off-diagonal block are marked by the [x]'s in the following scheme: | * * * [x] [x] [x]| | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. |[x] [x] [x] * * * | |[x] [x] [x] * * * | |[x] [x] [x] * * * | In terms of the columns of A, the first N1 columns are rotated 'against' the remaining N-N1 columns, trying to increase the angle between the corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter. The number of sweeps is given in NSWEEP and the orthogonality threshold is given in TOL.
◆ la_ctgsja()
| pure subroutine, public la_lapack_svd_comp::la_ctgsja | ( | character, intent(in) | jobu, |
| character, intent(in) | jobv, | ||
| character, intent(in) | jobq, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | p, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| integer(ilp), intent(in) | l, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(sp), dimension(ldb,*), intent(inout) | b, | ||
| integer(ilp), intent(in) | ldb, | ||
| real(sp), intent(in) | tola, | ||
| real(sp), intent(in) | tolb, | ||
| real(sp), dimension(*), intent(out) | alpha, | ||
| real(sp), dimension(*), intent(out) | beta, | ||
| complex(sp), dimension(ldu,*), intent(inout) | u, | ||
| integer(ilp), intent(in) | ldu, | ||
| complex(sp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| complex(sp), dimension(ldq,*), intent(inout) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| complex(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | ncycle, | ||
| integer(ilp), intent(out) | 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 A*Q = D1( 0 R ), V**H B*Q = 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) ), 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 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.
◆ la_cungbr()
| pure subroutine, public la_lapack_svd_comp::la_cungbr | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(sp), dimension(*), intent(in) | tau, | ||
| complex(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**H. Q and P**H 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 P**H is of order N: if k < n, P**H = G(k) . . . G(2) G(1) and CUNGBR returns the first m rows of P**H, where n >= m >= k; if k >= n, P**H = G(n-1) . . . G(2) G(1) and CUNGBR returns P**H as an N-by-N matrix.
◆ la_cunmbr()
| pure subroutine, public la_lapack_svd_comp::la_cunmbr | ( | character, intent(in) | vect, |
| character, intent(in) | side, | ||
| character, intent(in) | trans, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| complex(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(sp), dimension(*), intent(in) | tau, | ||
| complex(sp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| complex(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**H Here Q and P**H are the unitary matrices determined by CGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q and P**H 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).
◆ la_dgbbrd()
| pure subroutine, public la_lapack_svd_comp::la_dgbbrd | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | ncc, | ||
| integer(ilp), intent(in) | kl, | ||
| integer(ilp), intent(in) | ku, | ||
| real(dp), dimension(ldab,*), intent(inout) | ab, | ||
| integer(ilp), intent(in) | ldab, | ||
| real(dp), dimension(*), intent(out) | d, | ||
| real(dp), dimension(*), intent(out) | e, | ||
| real(dp), dimension(ldq,*), intent(out) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| real(dp), dimension(ldpt,*), intent(out) | pt, | ||
| integer(ilp), intent(in) | ldpt, | ||
| real(dp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| real(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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 P**T, or computes Q**T*C for a given matrix C.
◆ la_dgebd2()
| pure subroutine, public la_lapack_svd_comp::la_dgebd2 | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(*), intent(out) | d, | ||
| real(dp), dimension(*), intent(out) | e, | ||
| real(dp), dimension(*), intent(out) | tauq, | ||
| real(dp), dimension(*), intent(out) | taup, | ||
| real(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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.
◆ la_dgebrd()
| pure subroutine, public la_lapack_svd_comp::la_dgebrd | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(*), intent(out) | d, | ||
| real(dp), dimension(*), intent(out) | e, | ||
| real(dp), dimension(*), intent(out) | tauq, | ||
| real(dp), dimension(*), intent(out) | taup, | ||
| real(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_dgsvj0()
| pure subroutine, public la_lapack_svd_comp::la_dgsvj0 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(n), intent(inout) | d, | ||
| real(dp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| real(dp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(dp), intent(in) | eps, | ||
| real(dp), intent(in) | sfmin, | ||
| real(dp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| real(dp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_dgsvj1()
| pure subroutine, public la_lapack_svd_comp::la_dgsvj1 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | n1, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(n), intent(inout) | d, | ||
| real(dp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| real(dp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(dp), intent(in) | eps, | ||
| real(dp), intent(in) | sfmin, | ||
| real(dp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| real(dp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
DGSVJ1 applies few sweeps of Jacobi rotations in the column space of the input M-by-N matrix A. The pivot pairs are taken from the (1,2) off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The block-entries (tiles) of the (1,2) off-diagonal block are marked by the [x]'s in the following scheme: | * * * [x] [x] [x]| | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. |[x] [x] [x] * * * | |[x] [x] [x] * * * | |[x] [x] [x] * * * | In terms of the columns of A, the first N1 columns are rotated 'against' the remaining N-N1 columns, trying to increase the angle between the corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter. The number of sweeps is given in NSWEEP and the orthogonality threshold is given in TOL.
◆ la_dorgbr()
| pure subroutine, public la_lapack_svd_comp::la_dorgbr | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(*), intent(in) | tau, | ||
| real(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**T. Q and P**T 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 P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as an N-by-N matrix.
◆ la_dormbr()
| pure subroutine, public la_lapack_svd_comp::la_dormbr | ( | character, intent(in) | vect, |
| character, intent(in) | side, | ||
| character, intent(in) | trans, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(*), intent(in) | tau, | ||
| real(dp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| real(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**T Here Q and P**T are the orthogonal matrices determined by DGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T 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).
◆ la_dtgsja()
| pure subroutine, public la_lapack_svd_comp::la_dtgsja | ( | character, intent(in) | jobu, |
| character, intent(in) | jobv, | ||
| character, intent(in) | jobq, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | p, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| integer(ilp), intent(in) | l, | ||
| real(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(ldb,*), intent(inout) | b, | ||
| integer(ilp), intent(in) | ldb, | ||
| real(dp), intent(in) | tola, | ||
| real(dp), intent(in) | tolb, | ||
| real(dp), dimension(*), intent(out) | alpha, | ||
| real(dp), dimension(*), intent(out) | beta, | ||
| real(dp), dimension(ldu,*), intent(inout) | u, | ||
| integer(ilp), intent(in) | ldu, | ||
| real(dp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(dp), dimension(ldq,*), intent(inout) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| real(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | ncycle, | ||
| integer(ilp), intent(out) | 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 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.
◆ la_qgbbrd()
| pure subroutine, public la_lapack_svd_comp::la_qgbbrd | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | ncc, | ||
| integer(ilp), intent(in) | kl, | ||
| integer(ilp), intent(in) | ku, | ||
| real(qp), dimension(ldab,*), intent(inout) | ab, | ||
| integer(ilp), intent(in) | ldab, | ||
| real(qp), dimension(*), intent(out) | d, | ||
| real(qp), dimension(*), intent(out) | e, | ||
| real(qp), dimension(ldq,*), intent(out) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| real(qp), dimension(ldpt,*), intent(out) | pt, | ||
| integer(ilp), intent(in) | ldpt, | ||
| real(qp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| real(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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 P**T, or computes Q**T*C for a given matrix C.
◆ la_qgebd2()
| pure subroutine, public la_lapack_svd_comp::la_qgebd2 | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(*), intent(out) | d, | ||
| real(qp), dimension(*), intent(out) | e, | ||
| real(qp), dimension(*), intent(out) | tauq, | ||
| real(qp), dimension(*), intent(out) | taup, | ||
| real(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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.
◆ la_qgebrd()
| pure subroutine, public la_lapack_svd_comp::la_qgebrd | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(*), intent(out) | d, | ||
| real(qp), dimension(*), intent(out) | e, | ||
| real(qp), dimension(*), intent(out) | tauq, | ||
| real(qp), dimension(*), intent(out) | taup, | ||
| real(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_qgsvj0()
| pure subroutine, public la_lapack_svd_comp::la_qgsvj0 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(n), intent(inout) | d, | ||
| real(qp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| real(qp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(qp), intent(in) | eps, | ||
| real(qp), intent(in) | sfmin, | ||
| real(qp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| real(qp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_qgsvj1()
| pure subroutine, public la_lapack_svd_comp::la_qgsvj1 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | n1, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(n), intent(inout) | d, | ||
| real(qp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| real(qp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(qp), intent(in) | eps, | ||
| real(qp), intent(in) | sfmin, | ||
| real(qp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| real(qp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
QGSVJ1 applies few sweeps of Jacobi rotations in the column space of the input M-by-N matrix A. The pivot pairs are taken from the (1,2) off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The block-entries (tiles) of the (1,2) off-diagonal block are marked by the [x]'s in the following scheme: | * * * [x] [x] [x]| | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. |[x] [x] [x] * * * | |[x] [x] [x] * * * | |[x] [x] [x] * * * | In terms of the columns of A, the first N1 columns are rotated 'against' the remaining N-N1 columns, trying to increase the angle between the corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter. The number of sweeps is given in NSWEEP and the orthogonality threshold is given in TOL.
◆ la_qorgbr()
| pure subroutine, public la_lapack_svd_comp::la_qorgbr | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(*), intent(in) | tau, | ||
| real(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**T. Q and P**T 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 P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and QORGBR returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and QORGBR returns P**T as an N-by-N matrix.
◆ la_qormbr()
| pure subroutine, public la_lapack_svd_comp::la_qormbr | ( | character, intent(in) | vect, |
| character, intent(in) | side, | ||
| character, intent(in) | trans, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(*), intent(in) | tau, | ||
| real(qp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| real(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**T Here Q and P**T are the orthogonal matrices determined by QGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T 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).
◆ la_qtgsja()
| pure subroutine, public la_lapack_svd_comp::la_qtgsja | ( | character, intent(in) | jobu, |
| character, intent(in) | jobv, | ||
| character, intent(in) | jobq, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | p, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| integer(ilp), intent(in) | l, | ||
| real(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(ldb,*), intent(inout) | b, | ||
| integer(ilp), intent(in) | ldb, | ||
| real(qp), intent(in) | tola, | ||
| real(qp), intent(in) | tolb, | ||
| real(qp), dimension(*), intent(out) | alpha, | ||
| real(qp), dimension(*), intent(out) | beta, | ||
| real(qp), dimension(ldu,*), intent(inout) | u, | ||
| integer(ilp), intent(in) | ldu, | ||
| real(qp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(qp), dimension(ldq,*), intent(inout) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| real(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | ncycle, | ||
| integer(ilp), intent(out) | 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 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.
◆ la_sgbbrd()
| pure subroutine, public la_lapack_svd_comp::la_sgbbrd | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | ncc, | ||
| integer(ilp), intent(in) | kl, | ||
| integer(ilp), intent(in) | ku, | ||
| real(sp), dimension(ldab,*), intent(inout) | ab, | ||
| integer(ilp), intent(in) | ldab, | ||
| real(sp), dimension(*), intent(out) | d, | ||
| real(sp), dimension(*), intent(out) | e, | ||
| real(sp), dimension(ldq,*), intent(out) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| real(sp), dimension(ldpt,*), intent(out) | pt, | ||
| integer(ilp), intent(in) | ldpt, | ||
| real(sp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| real(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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 P**T, or computes Q**T*C for a given matrix C.
◆ la_sgebd2()
| pure subroutine, public la_lapack_svd_comp::la_sgebd2 | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(*), intent(out) | d, | ||
| real(sp), dimension(*), intent(out) | e, | ||
| real(sp), dimension(*), intent(out) | tauq, | ||
| real(sp), dimension(*), intent(out) | taup, | ||
| real(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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.
◆ la_sgebrd()
| pure subroutine, public la_lapack_svd_comp::la_sgebrd | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(*), intent(out) | d, | ||
| real(sp), dimension(*), intent(out) | e, | ||
| real(sp), dimension(*), intent(out) | tauq, | ||
| real(sp), dimension(*), intent(out) | taup, | ||
| real(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_sgsvj0()
| pure subroutine, public la_lapack_svd_comp::la_sgsvj0 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(n), intent(inout) | d, | ||
| real(sp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| real(sp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(sp), intent(in) | eps, | ||
| real(sp), intent(in) | sfmin, | ||
| real(sp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| real(sp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_sgsvj1()
| pure subroutine, public la_lapack_svd_comp::la_sgsvj1 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | n1, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(n), intent(inout) | d, | ||
| real(sp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| real(sp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(sp), intent(in) | eps, | ||
| real(sp), intent(in) | sfmin, | ||
| real(sp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| real(sp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
SGSVJ1 applies few sweeps of Jacobi rotations in the column space of the input M-by-N matrix A. The pivot pairs are taken from the (1,2) off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The block-entries (tiles) of the (1,2) off-diagonal block are marked by the [x]'s in the following scheme: | * * * [x] [x] [x]| | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. |[x] [x] [x] * * * | |[x] [x] [x] * * * | |[x] [x] [x] * * * | In terms of the columns of A, the first N1 columns are rotated 'against' the remaining N-N1 columns, trying to increase the angle between the corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter. The number of sweeps is given in NSWEEP and the orthogonality threshold is given in TOL.
◆ la_sorgbr()
| pure subroutine, public la_lapack_svd_comp::la_sorgbr | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(*), intent(in) | tau, | ||
| real(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**T. Q and P**T 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 P**T is of order N: if k < n, P**T = G(k) . . . G(2) G(1) and SORGBR returns the first m rows of P**T, where n >= m >= k; if k >= n, P**T = G(n-1) . . . G(2) G(1) and SORGBR returns P**T as an N-by-N matrix.
◆ la_sormbr()
| pure subroutine, public la_lapack_svd_comp::la_sormbr | ( | character, intent(in) | vect, |
| character, intent(in) | side, | ||
| character, intent(in) | trans, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(*), intent(in) | tau, | ||
| real(sp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| real(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**T Here Q and P**T are the orthogonal matrices determined by SGEBRD when reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and P**T 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).
◆ la_stgsja()
| pure subroutine, public la_lapack_svd_comp::la_stgsja | ( | character, intent(in) | jobu, |
| character, intent(in) | jobv, | ||
| character, intent(in) | jobq, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | p, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| integer(ilp), intent(in) | l, | ||
| real(sp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(sp), dimension(ldb,*), intent(inout) | b, | ||
| integer(ilp), intent(in) | ldb, | ||
| real(sp), intent(in) | tola, | ||
| real(sp), intent(in) | tolb, | ||
| real(sp), dimension(*), intent(out) | alpha, | ||
| real(sp), dimension(*), intent(out) | beta, | ||
| real(sp), dimension(ldu,*), intent(inout) | u, | ||
| integer(ilp), intent(in) | ldu, | ||
| real(sp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(sp), dimension(ldq,*), intent(inout) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| real(sp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | ncycle, | ||
| integer(ilp), intent(out) | 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 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.
◆ la_wgbbrd()
| pure subroutine, public la_lapack_svd_comp::la_wgbbrd | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | ncc, | ||
| integer(ilp), intent(in) | kl, | ||
| integer(ilp), intent(in) | ku, | ||
| complex(qp), dimension(ldab,*), intent(inout) | ab, | ||
| integer(ilp), intent(in) | ldab, | ||
| real(qp), dimension(*), intent(out) | d, | ||
| real(qp), dimension(*), intent(out) | e, | ||
| complex(qp), dimension(ldq,*), intent(out) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| complex(qp), dimension(ldpt,*), intent(out) | pt, | ||
| integer(ilp), intent(in) | ldpt, | ||
| complex(qp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| complex(qp), dimension(*), intent(out) | work, | ||
| real(qp), dimension(*), intent(out) | rwork, | ||
| integer(ilp), intent(out) | 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 P**H, or computes Q**H*C for a given matrix C.
◆ la_wgebd2()
| pure subroutine, public la_lapack_svd_comp::la_wgebd2 | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(*), intent(out) | d, | ||
| real(qp), dimension(*), intent(out) | e, | ||
| complex(qp), dimension(*), intent(out) | tauq, | ||
| complex(qp), dimension(*), intent(out) | taup, | ||
| complex(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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.
◆ la_wgebrd()
| pure subroutine, public la_lapack_svd_comp::la_wgebrd | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(qp), dimension(*), intent(out) | d, | ||
| real(qp), dimension(*), intent(out) | e, | ||
| complex(qp), dimension(*), intent(out) | tauq, | ||
| complex(qp), dimension(*), intent(out) | taup, | ||
| complex(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_wgsvj0()
| pure subroutine, public la_lapack_svd_comp::la_wgsvj0 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(qp), dimension(n), intent(inout) | d, | ||
| real(qp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| complex(qp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(qp), intent(in) | eps, | ||
| real(qp), intent(in) | sfmin, | ||
| real(qp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| complex(qp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_wgsvj1()
| pure subroutine, public la_lapack_svd_comp::la_wgsvj1 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | n1, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(qp), dimension(n), intent(inout) | d, | ||
| real(qp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| complex(qp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(qp), intent(in) | eps, | ||
| real(qp), intent(in) | sfmin, | ||
| real(qp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| complex(qp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
WGSVJ1 applies few sweeps of Jacobi rotations in the column space of the input M-by-N matrix A. The pivot pairs are taken from the (1,2) off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The block-entries (tiles) of the (1,2) off-diagonal block are marked by the [x]'s in the following scheme: | * * * [x] [x] [x]| | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. |[x] [x] [x] * * * | |[x] [x] [x] * * * | |[x] [x] [x] * * * | In terms of the columns of A, the first N1 columns are rotated 'against' the remaining N-N1 columns, trying to increase the angle between the corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter. The number of sweeps is given in NSWEEP and the orthogonality threshold is given in TOL.
◆ la_wtgsja()
| pure subroutine, public la_lapack_svd_comp::la_wtgsja | ( | character, intent(in) | jobu, |
| character, intent(in) | jobv, | ||
| character, intent(in) | jobq, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | p, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| integer(ilp), intent(in) | l, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(qp), dimension(ldb,*), intent(inout) | b, | ||
| integer(ilp), intent(in) | ldb, | ||
| real(qp), intent(in) | tola, | ||
| real(qp), intent(in) | tolb, | ||
| real(qp), dimension(*), intent(out) | alpha, | ||
| real(qp), dimension(*), intent(out) | beta, | ||
| complex(qp), dimension(ldu,*), intent(inout) | u, | ||
| integer(ilp), intent(in) | ldu, | ||
| complex(qp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| complex(qp), dimension(ldq,*), intent(inout) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| complex(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | ncycle, | ||
| integer(ilp), intent(out) | 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 A*Q = D1( 0 R ), V**H B*Q = 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) ), 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 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.
◆ la_wungbr()
| pure subroutine, public la_lapack_svd_comp::la_wungbr | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(qp), dimension(*), intent(in) | tau, | ||
| complex(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**H. Q and P**H 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 P**H is of order N: if k < n, P**H = G(k) . . . G(2) G(1) and WUNGBR returns the first m rows of P**H, where n >= m >= k; if k >= n, P**H = G(n-1) . . . G(2) G(1) and WUNGBR returns P**H as an N-by-N matrix.
◆ la_wunmbr()
| pure subroutine, public la_lapack_svd_comp::la_wunmbr | ( | character, intent(in) | vect, |
| character, intent(in) | side, | ||
| character, intent(in) | trans, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| complex(qp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(qp), dimension(*), intent(in) | tau, | ||
| complex(qp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| complex(qp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**H Here Q and P**H are the unitary matrices determined by WGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q and P**H 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).
◆ la_zgbbrd()
| pure subroutine, public la_lapack_svd_comp::la_zgbbrd | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | ncc, | ||
| integer(ilp), intent(in) | kl, | ||
| integer(ilp), intent(in) | ku, | ||
| complex(dp), dimension(ldab,*), intent(inout) | ab, | ||
| integer(ilp), intent(in) | ldab, | ||
| real(dp), dimension(*), intent(out) | d, | ||
| real(dp), dimension(*), intent(out) | e, | ||
| complex(dp), dimension(ldq,*), intent(out) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| complex(dp), dimension(ldpt,*), intent(out) | pt, | ||
| integer(ilp), intent(in) | ldpt, | ||
| complex(dp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| complex(dp), dimension(*), intent(out) | work, | ||
| real(dp), dimension(*), intent(out) | rwork, | ||
| integer(ilp), intent(out) | 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 P**H, or computes Q**H*C for a given matrix C.
◆ la_zgebd2()
| pure subroutine, public la_lapack_svd_comp::la_zgebd2 | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(*), intent(out) | d, | ||
| real(dp), dimension(*), intent(out) | e, | ||
| complex(dp), dimension(*), intent(out) | tauq, | ||
| complex(dp), dimension(*), intent(out) | taup, | ||
| complex(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | 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.
◆ la_zgebrd()
| pure subroutine, public la_lapack_svd_comp::la_zgebrd | ( | integer(ilp), intent(in) | m, |
| integer(ilp), intent(in) | n, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| real(dp), dimension(*), intent(out) | d, | ||
| real(dp), dimension(*), intent(out) | e, | ||
| complex(dp), dimension(*), intent(out) | tauq, | ||
| complex(dp), dimension(*), intent(out) | taup, | ||
| complex(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_zgsvj0()
| pure subroutine, public la_lapack_svd_comp::la_zgsvj0 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(dp), dimension(n), intent(inout) | d, | ||
| real(dp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| complex(dp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(dp), intent(in) | eps, | ||
| real(dp), intent(in) | sfmin, | ||
| real(dp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| complex(dp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
◆ la_zgsvj1()
| pure subroutine, public la_lapack_svd_comp::la_zgsvj1 | ( | character, intent(in) | jobv, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | n1, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(dp), dimension(n), intent(inout) | d, | ||
| real(dp), dimension(n), intent(inout) | sva, | ||
| integer(ilp), intent(in) | mv, | ||
| complex(dp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| real(dp), intent(in) | eps, | ||
| real(dp), intent(in) | sfmin, | ||
| real(dp), intent(in) | tol, | ||
| integer(ilp), intent(in) | nsweep, | ||
| complex(dp), dimension(lwork), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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.
ZGSVJ1 applies few sweeps of Jacobi rotations in the column space of the input M-by-N matrix A. The pivot pairs are taken from the (1,2) off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The block-entries (tiles) of the (1,2) off-diagonal block are marked by the [x]'s in the following scheme: | * * * [x] [x] [x]| | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks. | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block. |[x] [x] [x] * * * | |[x] [x] [x] * * * | |[x] [x] [x] * * * | In terms of the columns of A, the first N1 columns are rotated 'against' the remaining N-N1 columns, trying to increase the angle between the corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter. The number of sweeps is given in NSWEEP and the orthogonality threshold is given in TOL.
◆ la_ztgsja()
| pure subroutine, public la_lapack_svd_comp::la_ztgsja | ( | character, intent(in) | jobu, |
| character, intent(in) | jobv, | ||
| character, intent(in) | jobq, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | p, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| integer(ilp), intent(in) | l, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(dp), dimension(ldb,*), intent(inout) | b, | ||
| integer(ilp), intent(in) | ldb, | ||
| real(dp), intent(in) | tola, | ||
| real(dp), intent(in) | tolb, | ||
| real(dp), dimension(*), intent(out) | alpha, | ||
| real(dp), dimension(*), intent(out) | beta, | ||
| complex(dp), dimension(ldu,*), intent(inout) | u, | ||
| integer(ilp), intent(in) | ldu, | ||
| complex(dp), dimension(ldv,*), intent(inout) | v, | ||
| integer(ilp), intent(in) | ldv, | ||
| complex(dp), dimension(ldq,*), intent(inout) | q, | ||
| integer(ilp), intent(in) | ldq, | ||
| complex(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(out) | ncycle, | ||
| integer(ilp), intent(out) | 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 A*Q = D1( 0 R ), V**H B*Q = 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) ), 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 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.
◆ la_zungbr()
| pure subroutine, public la_lapack_svd_comp::la_zungbr | ( | character, intent(in) | vect, |
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(dp), dimension(*), intent(in) | tau, | ||
| complex(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**H. Q and P**H 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 P**H is of order N: if k < n, P**H = G(k) . . . G(2) G(1) and ZUNGBR returns the first m rows of P**H, where n >= m >= k; if k >= n, P**H = G(n-1) . . . G(2) G(1) and ZUNGBR returns P**H as an N-by-N matrix.
◆ la_zunmbr()
| pure subroutine, public la_lapack_svd_comp::la_zunmbr | ( | character, intent(in) | vect, |
| character, intent(in) | side, | ||
| character, intent(in) | trans, | ||
| integer(ilp), intent(in) | m, | ||
| integer(ilp), intent(in) | n, | ||
| integer(ilp), intent(in) | k, | ||
| complex(dp), dimension(lda,*), intent(inout) | a, | ||
| integer(ilp), intent(in) | lda, | ||
| complex(dp), dimension(*), intent(in) | tau, | ||
| complex(dp), dimension(ldc,*), intent(inout) | c, | ||
| integer(ilp), intent(in) | ldc, | ||
| complex(dp), dimension(*), intent(out) | work, | ||
| integer(ilp), intent(in) | lwork, | ||
| integer(ilp), intent(out) | 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 * P**H Here Q and P**H are the unitary matrices determined by ZGEBRD when reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q and P**H 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).
