version 1.2, 2010/12/21 13:53:58
|
version 1.9, 2014/01/27 09:24:37
|
Line 1
|
Line 1
|
|
*> \brief \b ZUNCSD |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZUNCSD + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zuncsd.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zuncsd.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zuncsd.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* RECURSIVE SUBROUTINE ZUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, |
|
* SIGNS, M, P, Q, X11, LDX11, X12, |
|
* LDX12, X21, LDX21, X22, LDX22, THETA, |
|
* U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, |
|
* LDV2T, WORK, LWORK, RWORK, LRWORK, |
|
* IWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS |
|
* INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, |
|
* $ LDX21, LDX22, LRWORK, LWORK, M, P, Q |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER IWORK( * ) |
|
* DOUBLE PRECISION THETA( * ) |
|
* DOUBLE PRECISION RWORK( * ) |
|
* COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), |
|
* $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), |
|
* $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, |
|
* $ * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZUNCSD computes the CS decomposition of an M-by-M partitioned |
|
*> unitary matrix X: |
|
*> |
|
*> [ I 0 0 | 0 0 0 ] |
|
*> [ 0 C 0 | 0 -S 0 ] |
|
*> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H |
|
*> X = [-----------] = [---------] [---------------------] [---------] . |
|
*> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ] |
|
*> [ 0 S 0 | 0 C 0 ] |
|
*> [ 0 0 I | 0 0 0 ] |
|
*> |
|
*> X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P, |
|
*> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are |
|
*> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in |
|
*> which R = MIN(P,M-P,Q,M-Q). |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] JOBU1 |
|
*> \verbatim |
|
*> JOBU1 is CHARACTER |
|
*> = 'Y': U1 is computed; |
|
*> otherwise: U1 is not computed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] JOBU2 |
|
*> \verbatim |
|
*> JOBU2 is CHARACTER |
|
*> = 'Y': U2 is computed; |
|
*> otherwise: U2 is not computed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] JOBV1T |
|
*> \verbatim |
|
*> JOBV1T is CHARACTER |
|
*> = 'Y': V1T is computed; |
|
*> otherwise: V1T is not computed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] JOBV2T |
|
*> \verbatim |
|
*> JOBV2T is CHARACTER |
|
*> = 'Y': V2T is computed; |
|
*> otherwise: V2T is not computed. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] TRANS |
|
*> \verbatim |
|
*> TRANS is CHARACTER |
|
*> = 'T': X, U1, U2, V1T, and V2T are stored in row-major |
|
*> order; |
|
*> otherwise: X, U1, U2, V1T, and V2T are stored in column- |
|
*> major order. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] SIGNS |
|
*> \verbatim |
|
*> SIGNS is CHARACTER |
|
*> = 'O': The lower-left block is made nonpositive (the |
|
*> "other" convention); |
|
*> otherwise: The upper-right block is made nonpositive (the |
|
*> "default" convention). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of rows and columns in X. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] P |
|
*> \verbatim |
|
*> P is INTEGER |
|
*> The number of rows in X11 and X12. 0 <= P <= M. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] Q |
|
*> \verbatim |
|
*> Q is INTEGER |
|
*> The number of columns in X11 and X21. 0 <= Q <= M. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] X11 |
|
*> \verbatim |
|
*> X11 is COMPLEX*16 array, dimension (LDX11,Q) |
|
*> On entry, part of the unitary matrix whose CSD is desired. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDX11 |
|
*> \verbatim |
|
*> LDX11 is INTEGER |
|
*> The leading dimension of X11. LDX11 >= MAX(1,P). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] X12 |
|
*> \verbatim |
|
*> X12 is COMPLEX*16 array, dimension (LDX12,M-Q) |
|
*> On entry, part of the unitary matrix whose CSD is desired. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDX12 |
|
*> \verbatim |
|
*> LDX12 is INTEGER |
|
*> The leading dimension of X12. LDX12 >= MAX(1,P). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] X21 |
|
*> \verbatim |
|
*> X21 is COMPLEX*16 array, dimension (LDX21,Q) |
|
*> On entry, part of the unitary matrix whose CSD is desired. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDX21 |
|
*> \verbatim |
|
*> LDX21 is INTEGER |
|
*> The leading dimension of X11. LDX21 >= MAX(1,M-P). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] X22 |
|
*> \verbatim |
|
*> X22 is COMPLEX*16 array, dimension (LDX22,M-Q) |
|
*> On entry, part of the unitary matrix whose CSD is desired. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDX22 |
|
*> \verbatim |
|
*> LDX22 is INTEGER |
|
*> The leading dimension of X11. LDX22 >= MAX(1,M-P). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] THETA |
|
*> \verbatim |
|
*> THETA is DOUBLE PRECISION array, dimension (R), in which R = |
|
*> MIN(P,M-P,Q,M-Q). |
|
*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and |
|
*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] U1 |
|
*> \verbatim |
|
*> U1 is COMPLEX*16 array, dimension (P) |
|
*> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDU1 |
|
*> \verbatim |
|
*> LDU1 is INTEGER |
|
*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= |
|
*> MAX(1,P). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] U2 |
|
*> \verbatim |
|
*> U2 is COMPLEX*16 array, dimension (M-P) |
|
*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary |
|
*> matrix U2. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDU2 |
|
*> \verbatim |
|
*> LDU2 is INTEGER |
|
*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= |
|
*> MAX(1,M-P). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] V1T |
|
*> \verbatim |
|
*> V1T is COMPLEX*16 array, dimension (Q) |
|
*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary |
|
*> matrix V1**H. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDV1T |
|
*> \verbatim |
|
*> LDV1T is INTEGER |
|
*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= |
|
*> MAX(1,Q). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] V2T |
|
*> \verbatim |
|
*> V2T is COMPLEX*16 array, dimension (M-Q) |
|
*> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary |
|
*> matrix V2**H. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDV2T |
|
*> \verbatim |
|
*> LDV2T is INTEGER |
|
*> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >= |
|
*> MAX(1,M-Q). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) |
|
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> The dimension of the array WORK. |
|
*> |
|
*> If LWORK = -1, then a workspace query is assumed; the routine |
|
*> only calculates the optimal size of the WORK array, returns |
|
*> this value as the first entry of the work array, and no error |
|
*> message related to LWORK is issued by XERBLA. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] RWORK |
|
*> \verbatim |
|
*> RWORK is DOUBLE PRECISION array, dimension MAX(1,LRWORK) |
|
*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. |
|
*> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1), |
|
*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), |
|
*> define the matrix in intermediate bidiagonal-block form |
|
*> remaining after nonconvergence. INFO specifies the number |
|
*> of nonzero PHI's. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LRWORK |
|
*> \verbatim |
|
*> LRWORK is INTEGER |
|
*> The dimension of the array RWORK. |
|
*> |
|
*> If LRWORK = -1, then a workspace query is assumed; the routine |
|
*> only calculates the optimal size of the RWORK array, returns |
|
*> this value as the first entry of the work array, and no error |
|
*> message related to LRWORK is issued by XERBLA. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit. |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
|
*> > 0: ZBBCSD did not converge. See the description of RWORK |
|
*> above for details. |
|
*> \endverbatim |
|
* |
|
*> \par References: |
|
* ================ |
|
*> |
|
*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. |
|
*> Algorithms, 50(1):33-65, 2009. |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date November 2013 |
|
* |
|
*> \ingroup complex16OTHERcomputational |
|
* |
|
* ===================================================================== |
RECURSIVE SUBROUTINE ZUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, |
RECURSIVE SUBROUTINE ZUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, |
$ SIGNS, M, P, Q, X11, LDX11, X12, |
$ SIGNS, M, P, Q, X11, LDX11, X12, |
$ LDX12, X21, LDX21, X22, LDX22, THETA, |
$ LDX12, X21, LDX21, X22, LDX22, THETA, |
$ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, |
$ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, |
$ LDV2T, WORK, LWORK, RWORK, LRWORK, |
$ LDV2T, WORK, LWORK, RWORK, LRWORK, |
$ IWORK, INFO ) |
$ IWORK, INFO ) |
IMPLICIT NONE |
|
* |
|
* -- LAPACK routine (version 3.3.0) -- |
|
* |
|
* -- Contributed by Brian Sutton of the Randolph-Macon College -- |
|
* -- November 2010 |
|
* |
* |
|
* -- LAPACK computational routine (version 3.5.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
|
* November 2013 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS |
CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS |
Line 29
|
Line 340
|
$ * ) |
$ * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZUNCSD computes the CS decomposition of an M-by-M partitioned |
|
* unitary matrix X: |
|
* |
|
* [ I 0 0 | 0 0 0 ] |
|
* [ 0 C 0 | 0 -S 0 ] |
|
* [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H |
|
* X = [-----------] = [---------] [---------------------] [---------] . |
|
* [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ] |
|
* [ 0 S 0 | 0 C 0 ] |
|
* [ 0 0 I | 0 0 0 ] |
|
* |
|
* X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P, |
|
* (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are |
|
* R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in |
|
* which R = MIN(P,M-P,Q,M-Q). |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* JOBU1 (input) CHARACTER |
|
* = 'Y': U1 is computed; |
|
* otherwise: U1 is not computed. |
|
* |
|
* JOBU2 (input) CHARACTER |
|
* = 'Y': U2 is computed; |
|
* otherwise: U2 is not computed. |
|
* |
|
* JOBV1T (input) CHARACTER |
|
* = 'Y': V1T is computed; |
|
* otherwise: V1T is not computed. |
|
* |
|
* JOBV2T (input) CHARACTER |
|
* = 'Y': V2T is computed; |
|
* otherwise: V2T is not computed. |
|
* |
|
* TRANS (input) CHARACTER |
|
* = 'T': X, U1, U2, V1T, and V2T are stored in row-major |
|
* order; |
|
* otherwise: X, U1, U2, V1T, and V2T are stored in column- |
|
* major order. |
|
* |
|
* SIGNS (input) CHARACTER |
|
* = 'O': The lower-left block is made nonpositive (the |
|
* "other" convention); |
|
* otherwise: The upper-right block is made nonpositive (the |
|
* "default" convention). |
|
* |
|
* M (input) INTEGER |
|
* The number of rows and columns in X. |
|
* |
|
* P (input) INTEGER |
|
* The number of rows in X11 and X12. 0 <= P <= M. |
|
* |
|
* Q (input) INTEGER |
|
* The number of columns in X11 and X21. 0 <= Q <= M. |
|
* |
|
* X (input/workspace) COMPLEX*16 array, dimension (LDX,M) |
|
* On entry, the unitary matrix whose CSD is desired. |
|
* |
|
* LDX (input) INTEGER |
|
* The leading dimension of X. LDX >= MAX(1,M). |
|
* |
|
* THETA (output) DOUBLE PRECISION array, dimension (R), in which R = |
|
* MIN(P,M-P,Q,M-Q). |
|
* C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and |
|
* S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). |
|
* |
|
* U1 (output) COMPLEX*16 array, dimension (P) |
|
* If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1. |
|
* |
|
* LDU1 (input) INTEGER |
|
* The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= |
|
* MAX(1,P). |
|
* |
|
* U2 (output) COMPLEX*16 array, dimension (M-P) |
|
* If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary |
|
* matrix U2. |
|
* |
|
* LDU2 (input) INTEGER |
|
* The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= |
|
* MAX(1,M-P). |
|
* |
|
* V1T (output) COMPLEX*16 array, dimension (Q) |
|
* If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary |
|
* matrix V1**H. |
|
* |
|
* LDV1T (input) INTEGER |
|
* The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= |
|
* MAX(1,Q). |
|
* |
|
* V2T (output) COMPLEX*16 array, dimension (M-Q) |
|
* If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary |
|
* matrix V2**H. |
|
* |
|
* LDV2T (input) INTEGER |
|
* The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >= |
|
* MAX(1,M-Q). |
|
* |
|
* WORK (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK)) |
|
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
* |
|
* LWORK (input) INTEGER |
|
* The dimension of the array WORK. |
|
* |
|
* If LWORK = -1, then a workspace query is assumed; the routine |
|
* only calculates the optimal size of the WORK array, returns |
|
* this value as the first entry of the work array, and no error |
|
* message related to LWORK is issued by XERBLA. |
|
* |
|
* RWORK (workspace) DOUBLE PRECISION array, dimension MAX(1,LRWORK) |
|
* On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. |
|
* If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1), |
|
* ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), |
|
* define the matrix in intermediate bidiagonal-block form |
|
* remaining after nonconvergence. INFO specifies the number |
|
* of nonzero PHI's. |
|
* |
|
* LRWORK (input) INTEGER |
|
* The dimension of the array RWORK. |
|
* |
|
* If LRWORK = -1, then a workspace query is assumed; the routine |
|
* only calculates the optimal size of the RWORK array, returns |
|
* this value as the first entry of the work array, and no error |
|
* message related to LRWORK is issued by XERBLA. |
|
* |
|
* IWORK (workspace) INTEGER array, dimension (M-Q) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit. |
|
* < 0: if INFO = -i, the i-th argument had an illegal value. |
|
* > 0: ZBBCSD did not converge. See the description of RWORK |
|
* above for details. |
|
* |
|
* Reference |
|
* ========= |
|
* |
|
* [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. |
|
* Algorithms, 50(1):33-65, 2009. |
|
* |
|
* =================================================================== |
* =================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
DOUBLE PRECISION REALONE |
COMPLEX*16 ONE, ZERO |
PARAMETER ( REALONE = 1.0D0 ) |
PARAMETER ( ONE = (1.0D0,0.0D0), |
COMPLEX*16 NEGONE, ONE, PIOVER2, ZERO |
|
PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0), |
|
$ PIOVER2 = 1.57079632679489662D0, |
|
$ ZERO = (0.0D0,0.0D0) ) |
$ ZERO = (0.0D0,0.0D0) ) |
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
Line 190
|
Line 356
|
$ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, |
$ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, |
$ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, |
$ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, |
$ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, |
$ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, |
$ LORGQRWORKOPT, LWORKMIN, LWORKOPT |
$ LORGQRWORKOPT, LWORKMIN, LWORKOPT, P1, Q1 |
LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, |
LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, |
$ WANTV1T, WANTV2T |
$ WANTV1T, WANTV2T |
INTEGER LRWORKMIN, LRWORKOPT |
INTEGER LRWORKMIN, LRWORKOPT |
Line 205
|
Line 371
|
EXTERNAL LSAME |
EXTERNAL LSAME |
* .. |
* .. |
* .. Intrinsic Functions |
* .. Intrinsic Functions |
INTRINSIC COS, INT, MAX, MIN, SIN |
INTRINSIC INT, MAX, MIN |
* .. |
* .. |
* .. Executable Statements .. |
* .. Executable Statements .. |
* |
* |
Line 226
|
Line 392
|
INFO = -8 |
INFO = -8 |
ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN |
ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN |
INFO = -9 |
INFO = -9 |
ELSE IF( ( COLMAJOR .AND. LDX11 .LT. MAX(1,P) ) .OR. |
ELSE IF ( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN |
$ ( .NOT.COLMAJOR .AND. LDX11 .LT. MAX(1,Q) ) ) THEN |
INFO = -11 |
INFO = -11 |
ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN |
|
INFO = -11 |
|
ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN |
|
INFO = -13 |
|
ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN |
|
INFO = -13 |
|
ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN |
|
INFO = -15 |
|
ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN |
|
INFO = -15 |
|
ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN |
|
INFO = -17 |
|
ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN |
|
INFO = -17 |
ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN |
ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN |
INFO = -14 |
INFO = -20 |
ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN |
ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN |
INFO = -16 |
INFO = -22 |
ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN |
ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN |
INFO = -18 |
INFO = -24 |
ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN |
ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN |
INFO = -20 |
INFO = -26 |
END IF |
END IF |
* |
* |
* Work with transpose if convenient |
* Work with transpose if convenient |
Line 292
|
Line 471
|
IB22D = IB21E + MAX( 1, Q - 1 ) |
IB22D = IB21E + MAX( 1, Q - 1 ) |
IB22E = IB22D + MAX( 1, Q ) |
IB22E = IB22D + MAX( 1, Q ) |
IBBCSD = IB22E + MAX( 1, Q - 1 ) |
IBBCSD = IB22E + MAX( 1, Q - 1 ) |
CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 0, |
CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, |
$ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, 0, |
$ THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, |
$ 0, 0, 0, 0, 0, 0, 0, RWORK, -1, CHILDINFO ) |
$ V2T, LDV2T, THETA, THETA, THETA, THETA, THETA, |
|
$ THETA, THETA, THETA, RWORK, -1, CHILDINFO ) |
LBBCSDWORKOPT = INT( RWORK(1) ) |
LBBCSDWORKOPT = INT( RWORK(1) ) |
LBBCSDWORKMIN = LBBCSDWORKOPT |
LBBCSDWORKMIN = LBBCSDWORKOPT |
LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1 |
LRWORKOPT = IBBCSD + LBBCSDWORKOPT - 1 |
Line 308
|
Line 488
|
ITAUQ1 = ITAUP2 + MAX( 1, M - P ) |
ITAUQ1 = ITAUP2 + MAX( 1, M - P ) |
ITAUQ2 = ITAUQ1 + MAX( 1, Q ) |
ITAUQ2 = ITAUQ1 + MAX( 1, Q ) |
IORGQR = ITAUQ2 + MAX( 1, M - Q ) |
IORGQR = ITAUQ2 + MAX( 1, M - Q ) |
CALL ZUNGQR( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, |
CALL ZUNGQR( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGQRWORKOPT = INT( WORK(1) ) |
LORGQRWORKOPT = INT( WORK(1) ) |
LORGQRWORKMIN = MAX( 1, M - Q ) |
LORGQRWORKMIN = MAX( 1, M - Q ) |
IORGLQ = ITAUQ2 + MAX( 1, M - Q ) |
IORGLQ = ITAUQ2 + MAX( 1, M - Q ) |
CALL ZUNGLQ( M-Q, M-Q, M-Q, 0, MAX(1,M-Q), 0, WORK, -1, |
CALL ZUNGLQ( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGLQWORKOPT = INT( WORK(1) ) |
LORGLQWORKOPT = INT( WORK(1) ) |
LORGLQWORKMIN = MAX( 1, M - Q ) |
LORGLQWORKMIN = MAX( 1, M - Q ) |
IORBDB = ITAUQ2 + MAX( 1, M - Q ) |
IORBDB = ITAUQ2 + MAX( 1, M - Q ) |
CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, |
CALL ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, |
$ X21, LDX21, X22, LDX22, 0, 0, 0, 0, 0, 0, WORK, |
$ X21, LDX21, X22, LDX22, THETA, THETA, U1, U2, |
$ -1, CHILDINFO ) |
$ V1T, V2T, WORK, -1, CHILDINFO ) |
LORBDBWORKOPT = INT( WORK(1) ) |
LORBDBWORKOPT = INT( WORK(1) ) |
LORBDBWORKMIN = LORBDBWORKOPT |
LORBDBWORKMIN = LORBDBWORKOPT |
LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, |
LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, |
$ IORBDB + LORBDBWORKOPT ) - 1 |
$ IORBDB + LORBDBWORKOPT ) - 1 |
LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, |
LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, |
$ IORBDB + LORBDBWORKMIN ) - 1 |
$ IORBDB + LORBDBWORKMIN ) - 1 |
WORK(1) = LWORKOPT |
WORK(1) = MAX(LWORKOPT,LWORKMIN) |
* |
* |
IF( LWORK .LT. LWORKMIN |
IF( LWORK .LT. LWORKMIN |
$ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN |
$ .AND. .NOT. ( LQUERY .OR. LRQUERY ) ) THEN |
Line 385
|
Line 565
|
END IF |
END IF |
IF( WANTV2T .AND. M-Q .GT. 0 ) THEN |
IF( WANTV2T .AND. M-Q .GT. 0 ) THEN |
CALL ZLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) |
CALL ZLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) |
CALL ZLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, |
IF( M-P .GT. Q) THEN |
$ V2T(P+1,P+1), LDV2T ) |
CALL ZLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, |
CALL ZUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), |
$ V2T(P+1,P+1), LDV2T ) |
$ WORK(IORGLQ), LORGLQWORK, INFO ) |
END IF |
|
IF( M .GT. Q ) THEN |
|
CALL ZUNGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), |
|
$ WORK(IORGLQ), LORGLQWORK, INFO ) |
|
END IF |
END IF |
END IF |
ELSE |
ELSE |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
Line 413
|
Line 597
|
$ WORK(IORGQR), LORGQRWORK, INFO ) |
$ WORK(IORGQR), LORGQRWORK, INFO ) |
END IF |
END IF |
IF( WANTV2T .AND. M-Q .GT. 0 ) THEN |
IF( WANTV2T .AND. M-Q .GT. 0 ) THEN |
|
P1 = MIN( P+1, M ) |
|
Q1 = MIN( Q+1, M ) |
CALL ZLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) |
CALL ZLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) |
CALL ZLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, |
IF( M .GT. P+Q ) THEN |
$ V2T(P+1,P+1), LDV2T ) |
CALL ZLACPY( 'L', M-P-Q, M-P-Q, X22(P1,Q1), LDX22, |
|
$ V2T(P+1,P+1), LDV2T ) |
|
END IF |
CALL ZUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), |
CALL ZUNGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), |
$ WORK(IORGQR), LORGQRWORK, INFO ) |
$ WORK(IORGQR), LORGQRWORK, INFO ) |
END IF |
END IF |