version 1.1, 2014/01/27 09:24:35
|
version 1.10, 2023/08/07 08:39:01
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DORCSD2BY1 + dependencies |
*> Download DORCSD2BY1 + dependencies |
Line 21
|
Line 21
|
* SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, |
* SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, |
* X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, |
* X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, |
* LDV1T, WORK, LWORK, IWORK, INFO ) |
* LDV1T, WORK, LWORK, IWORK, INFO ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER JOBU1, JOBU2, JOBV1T |
* CHARACTER JOBU1, JOBU2, JOBV1T |
* INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21, |
* INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21, |
Line 33
|
Line 33
|
* $ X11(LDX11,*), X21(LDX21,*) |
* $ X11(LDX11,*), X21(LDX21,*) |
* INTEGER IWORK(*) |
* INTEGER IWORK(*) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
*> ============= |
* ============= |
*> |
*> |
*>\verbatim |
*>\verbatim |
*> Purpose: |
|
*> ======== |
|
*> |
*> |
*> DORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with |
*> DORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with |
*> orthonormal columns that has been partitioned into a 2-by-1 block |
*> orthonormal columns that has been partitioned into a 2-by-1 block |
*> structure: |
*> structure: |
*> |
*> |
*> [ I 0 0 ] |
*> [ I1 0 0 ] |
*> [ 0 C 0 ] |
*> [ 0 C 0 ] |
*> [ X11 ] [ U1 | ] [ 0 0 0 ] |
*> [ X11 ] [ U1 | ] [ 0 0 0 ] |
*> X = [-----] = [---------] [----------] V1**T . |
*> X = [-----] = [---------] [----------] V1**T . |
*> [ X21 ] [ | U2 ] [ 0 0 0 ] |
*> [ X21 ] [ | U2 ] [ 0 0 0 ] |
*> [ 0 S 0 ] |
*> [ 0 S 0 ] |
*> [ 0 0 I ] |
*> [ 0 0 I2] |
*> |
|
*> X11 is P-by-Q. The orthogonal 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 |
*> X11 is P-by-Q. The orthogonal matrices U1, U2, and V1 are P-by-P, |
|
*> (M-P)-by-(M-P), and Q-by-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). I1 is a K1-by-K1 identity matrix and I2 is a |
|
*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0). |
|
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
* ========== |
* ========== |
Line 67
|
Line 65
|
*> \param[in] JOBU1 |
*> \param[in] JOBU1 |
*> \verbatim |
*> \verbatim |
*> JOBU1 is CHARACTER |
*> JOBU1 is CHARACTER |
*> = 'Y': U1 is computed; |
*> = 'Y': U1 is computed; |
*> otherwise: U1 is not computed. |
*> otherwise: U1 is not computed. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] JOBU2 |
*> \param[in] JOBU2 |
*> \verbatim |
*> \verbatim |
*> JOBU2 is CHARACTER |
*> JOBU2 is CHARACTER |
*> = 'Y': U2 is computed; |
*> = 'Y': U2 is computed; |
*> otherwise: U2 is not computed. |
*> otherwise: U2 is not computed. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] JOBV1T |
*> \param[in] JOBV1T |
*> \verbatim |
*> \verbatim |
*> JOBV1T is CHARACTER |
*> JOBV1T is CHARACTER |
*> = 'Y': V1T is computed; |
*> = 'Y': V1T is computed; |
*> otherwise: V1T is not computed. |
*> otherwise: V1T is not computed. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] M |
*> \param[in] M |
*> \verbatim |
*> \verbatim |
*> M is INTEGER |
*> M is INTEGER |
*> The number of rows and columns in X. |
*> The number of rows in X. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] P |
*> \param[in] P |
*> \verbatim |
*> \verbatim |
*> P is INTEGER |
*> P is INTEGER |
*> The number of rows in X11 and X12. 0 <= P <= M. |
*> The number of rows in X11. 0 <= P <= M. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] Q |
*> \param[in] Q |
*> \verbatim |
*> \verbatim |
*> Q is INTEGER |
*> Q is INTEGER |
*> The number of columns in X11 and X21. 0 <= Q <= M. |
*> The number of columns in X11 and X21. 0 <= Q <= M. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] X11 |
*> \param[in,out] X11 |
*> \verbatim |
*> \verbatim |
*> X11 is DOUBLE PRECISION array, dimension (LDX11,Q) |
*> X11 is DOUBLE PRECISION array, dimension (LDX11,Q) |
*> On entry, part of the orthogonal matrix whose CSD is |
*> On entry, part of the orthogonal matrix whose CSD is desired. |
*> desired. |
|
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDX11 |
*> \param[in] LDX11 |
*> \verbatim |
*> \verbatim |
*> LDX11 is INTEGER |
*> LDX11 is INTEGER |
*> The leading dimension of X11. LDX11 >= MAX(1,P). |
*> The leading dimension of X11. LDX11 >= MAX(1,P). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] X21 |
*> \param[in,out] X21 |
*> \verbatim |
*> \verbatim |
*> X21 is DOUBLE PRECISION array, dimension (LDX21,Q) |
*> X21 is DOUBLE PRECISION array, dimension (LDX21,Q) |
*> On entry, part of the orthogonal matrix whose CSD is |
*> On entry, part of the orthogonal matrix whose CSD is desired. |
*> desired. |
|
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDX21 |
*> \param[in] LDX21 |
*> \verbatim |
*> \verbatim |
*> LDX21 is INTEGER |
*> LDX21 is INTEGER |
*> The leading dimension of X21. LDX21 >= MAX(1,M-P). |
*> The leading dimension of X21. LDX21 >= MAX(1,M-P). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] THETA |
*> \param[out] THETA |
*> \verbatim |
*> \verbatim |
*> THETA is DOUBLE PRECISION array, dimension (R), in which R = |
*> THETA is DOUBLE PRECISION array, dimension (R), in which R = |
*> MIN(P,M-P,Q,M-Q). |
*> MIN(P,M-P,Q,M-Q). |
*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and |
*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and |
*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). |
*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] U1 |
*> \param[out] U1 |
*> \verbatim |
*> \verbatim |
*> U1 is DOUBLE PRECISION array, dimension (P) |
*> U1 is DOUBLE PRECISION array, dimension (P) |
*> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1. |
*> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDU1 |
*> \param[in] LDU1 |
*> \verbatim |
*> \verbatim |
*> LDU1 is INTEGER |
*> LDU1 is INTEGER |
*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= |
*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >= |
*> MAX(1,P). |
*> MAX(1,P). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] U2 |
*> \param[out] U2 |
*> \verbatim |
*> \verbatim |
*> U2 is DOUBLE PRECISION array, dimension (M-P) |
*> U2 is DOUBLE PRECISION array, dimension (M-P) |
*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal |
*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal |
*> matrix U2. |
*> matrix U2. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDU2 |
*> \param[in] LDU2 |
*> \verbatim |
*> \verbatim |
*> LDU2 is INTEGER |
*> LDU2 is INTEGER |
*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= |
*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >= |
*> MAX(1,M-P). |
*> MAX(1,M-P). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] V1T |
*> \param[out] V1T |
*> \verbatim |
*> \verbatim |
*> V1T is DOUBLE PRECISION array, dimension (Q) |
*> V1T is DOUBLE PRECISION array, dimension (Q) |
*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal |
*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal |
*> matrix V1**T. |
*> matrix V1**T. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDV1T |
*> \param[in] LDV1T |
*> \verbatim |
*> \verbatim |
*> LDV1T is INTEGER |
*> LDV1T is INTEGER |
*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= |
*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= |
*> MAX(1,Q). |
*> MAX(1,Q). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[out] WORK |
*> \param[out] WORK |
*> \verbatim |
*> \verbatim |
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
*> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1), |
*> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1), |
*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), |
*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R), |
*> define the matrix in intermediate bidiagonal-block form |
*> define the matrix in intermediate bidiagonal-block form |
*> remaining after nonconvergence. INFO specifies the number |
*> remaining after nonconvergence. INFO specifies the number |
*> of nonzero PHI's. |
*> of nonzero PHI's. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LWORK |
*> \param[in] LWORK |
*> \verbatim |
*> \verbatim |
*> LWORK is INTEGER |
*> LWORK is INTEGER |
*> The dimension of the array WORK. |
*> 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 |
*> \endverbatim |
*> |
*> |
*> 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. |
|
*> \param[out] IWORK |
*> \param[out] IWORK |
*> \verbatim |
*> \verbatim |
*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) |
*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) |
Line 207
|
Line 204
|
*> \param[out] INFO |
*> \param[out] INFO |
*> \verbatim |
*> \verbatim |
*> INFO is INTEGER |
*> INFO is INTEGER |
*> = 0: successful exit. |
*> = 0: successful exit. |
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
*> > 0: DBBCSD did not converge. See the description of WORK |
*> > 0: DBBCSD did not converge. See the description of WORK |
*> above for details. |
*> above for details. |
*> \endverbatim |
*> \endverbatim |
* |
* |
|
*> \par References: |
|
* ================ |
|
*> |
|
*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. |
|
*> Algorithms, 50(1):33-65, 2009. |
|
* |
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
|
*> \date July 2012 |
|
* |
* |
*> \ingroup doubleOTHERcomputational |
*> \ingroup doubleOTHERcomputational |
* |
* |
*> \par References: |
|
* ================ |
|
*> |
|
*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer. |
|
*> Algorithms, 50(1):33-65, 2009. |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, |
SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, |
$ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, |
$ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, |
Line 240
|
Line 234
|
* -- LAPACK computational routine (3.5.0) -- |
* -- LAPACK computational routine (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..-- |
* July 2012 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER JOBU1, JOBU2, JOBV1T |
CHARACTER JOBU1, JOBU2, JOBV1T |
Line 253
|
Line 246
|
$ X11(LDX11,*), X21(LDX21,*) |
$ X11(LDX11,*), X21(LDX21,*) |
INTEGER IWORK(*) |
INTEGER IWORK(*) |
* .. |
* .. |
* |
* |
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 269
|
Line 262
|
$ LWORKMIN, LWORKOPT, R |
$ LWORKMIN, LWORKOPT, R |
LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T |
LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T |
* .. |
* .. |
|
* .. Local Arrays .. |
|
DOUBLE PRECISION DUM1(1), DUM2(1,1) |
|
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1, |
EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1, |
$ DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR, |
$ DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR, |
Line 301
|
Line 297
|
INFO = -8 |
INFO = -8 |
ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN |
ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN |
INFO = -10 |
INFO = -10 |
ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN |
ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN |
INFO = -13 |
INFO = -13 |
ELSE IF( WANTU2 .AND. LDU2 .LT. M - P ) THEN |
ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN |
INFO = -15 |
INFO = -15 |
ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN |
ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN |
INFO = -17 |
INFO = -17 |
END IF |
END IF |
* |
* |
Line 347
|
Line 343
|
IORBDB = ITAUQ1 + MAX( 1, Q ) |
IORBDB = ITAUQ1 + MAX( 1, Q ) |
IORGQR = ITAUQ1 + MAX( 1, Q ) |
IORGQR = ITAUQ1 + MAX( 1, Q ) |
IORGLQ = ITAUQ1 + MAX( 1, Q ) |
IORGLQ = ITAUQ1 + MAX( 1, Q ) |
|
LORGQRMIN = 1 |
|
LORGQROPT = 1 |
|
LORGLQMIN = 1 |
|
LORGLQOPT = 1 |
IF( R .EQ. Q ) THEN |
IF( R .EQ. Q ) THEN |
CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, |
CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, |
$ 0, 0, WORK, -1, CHILDINFO ) |
$ DUM1, DUM1, DUM1, DUM1, WORK, |
|
$ -1, CHILDINFO ) |
LORBDB = INT( WORK(1) ) |
LORBDB = INT( WORK(1) ) |
IF( P .GE. M-P ) THEN |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
CALL DORGQR( P, P, Q, U1, LDU1, 0, WORK(1), -1, |
CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGQRMIN = MAX( 1, P ) |
LORGQRMIN = MAX( LORGQRMIN, P ) |
LORGQROPT = INT( WORK(1) ) |
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
ELSE |
ENDIF |
CALL DORGQR( M-P, M-P, Q, U2, LDU2, 0, WORK(1), -1, |
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
$ CHILDINFO ) |
CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1), |
LORGQRMIN = MAX( 1, M-P ) |
$ -1, CHILDINFO ) |
LORGQROPT = INT( WORK(1) ) |
LORGQRMIN = MAX( LORGQRMIN, M-P ) |
|
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
|
END IF |
|
IF( WANTV1T .AND. Q .GT. 0 ) THEN |
|
CALL DORGLQ( Q-1, Q-1, Q-1, V1T, LDV1T, |
|
$ DUM1, WORK(1), -1, CHILDINFO ) |
|
LORGLQMIN = MAX( LORGLQMIN, Q-1 ) |
|
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) |
END IF |
END IF |
CALL DORGLQ( MAX(0,Q-1), MAX(0,Q-1), MAX(0,Q-1), V1T, LDV1T, |
|
$ 0, WORK(1), -1, CHILDINFO ) |
|
LORGLQMIN = MAX( 1, Q-1 ) |
|
LORGLQOPT = INT( WORK(1) ) |
|
CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, |
CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, |
$ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, 0, 1, 0, 0, |
$ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T, |
$ 0, 0, 0, 0, 0, 0, WORK(1), -1, CHILDINFO ) |
$ DUM2, 1, DUM1, DUM1, DUM1, |
|
$ DUM1, DUM1, DUM1, DUM1, |
|
$ DUM1, WORK(1), -1, CHILDINFO ) |
LBBCSD = INT( WORK(1) ) |
LBBCSD = INT( WORK(1) ) |
ELSE IF( R .EQ. P ) THEN |
ELSE IF( R .EQ. P ) THEN |
CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, |
CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, |
$ 0, 0, WORK(1), -1, CHILDINFO ) |
$ DUM1, DUM1, DUM1, DUM1, |
|
$ WORK(1), -1, CHILDINFO ) |
LORBDB = INT( WORK(1) ) |
LORBDB = INT( WORK(1) ) |
IF( P-1 .GE. M-P ) THEN |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, 0, WORK(1), |
CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, DUM1, |
|
$ WORK(1), -1, CHILDINFO ) |
|
LORGQRMIN = MAX( LORGQRMIN, P-1 ) |
|
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
|
END IF |
|
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
|
CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1), |
$ -1, CHILDINFO ) |
$ -1, CHILDINFO ) |
LORGQRMIN = MAX( 1, P-1 ) |
LORGQRMIN = MAX( LORGQRMIN, M-P ) |
LORGQROPT = INT( WORK(1) ) |
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
ELSE |
END IF |
CALL DORGQR( M-P, M-P, Q, U2, LDU2, 0, WORK(1), -1, |
IF( WANTV1T .AND. Q .GT. 0 ) THEN |
|
CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGQRMIN = MAX( 1, M-P ) |
LORGLQMIN = MAX( LORGLQMIN, Q ) |
LORGQROPT = INT( WORK(1) ) |
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) |
END IF |
END IF |
CALL DORGLQ( Q, Q, R, V1T, LDV1T, 0, WORK(1), -1, |
|
$ CHILDINFO ) |
|
LORGLQMIN = MAX( 1, Q ) |
|
LORGLQOPT = INT( WORK(1) ) |
|
CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, |
CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, |
$ 0, V1T, LDV1T, 0, 1, U1, LDU1, U2, LDU2, 0, 0, |
$ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1, |
$ 0, 0, 0, 0, 0, 0, WORK(1), -1, CHILDINFO ) |
$ U2, LDU2, DUM1, DUM1, DUM1, |
|
$ DUM1, DUM1, DUM1, DUM1, |
|
$ DUM1, WORK(1), -1, CHILDINFO ) |
LBBCSD = INT( WORK(1) ) |
LBBCSD = INT( WORK(1) ) |
ELSE IF( R .EQ. M-P ) THEN |
ELSE IF( R .EQ. M-P ) THEN |
CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, |
CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, |
$ 0, 0, WORK(1), -1, CHILDINFO ) |
$ DUM1, DUM1, DUM1, DUM1, |
|
$ WORK(1), -1, CHILDINFO ) |
LORBDB = INT( WORK(1) ) |
LORBDB = INT( WORK(1) ) |
IF( P .GE. M-P-1 ) THEN |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
CALL DORGQR( P, P, Q, U1, LDU1, 0, WORK(1), -1, |
CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGQRMIN = MAX( 1, P ) |
LORGQRMIN = MAX( LORGQRMIN, P ) |
LORGQROPT = INT( WORK(1) ) |
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
ELSE |
END IF |
CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, 0, |
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
$ WORK(1), -1, CHILDINFO ) |
CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, |
LORGQRMIN = MAX( 1, M-P-1 ) |
$ DUM1, WORK(1), -1, CHILDINFO ) |
LORGQROPT = INT( WORK(1) ) |
LORGQRMIN = MAX( LORGQRMIN, M-P-1 ) |
|
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
|
END IF |
|
IF( WANTV1T .AND. Q .GT. 0 ) THEN |
|
CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1, |
|
$ CHILDINFO ) |
|
LORGLQMIN = MAX( LORGLQMIN, Q ) |
|
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) |
END IF |
END IF |
CALL DORGLQ( Q, Q, R, V1T, LDV1T, 0, WORK(1), -1, |
|
$ CHILDINFO ) |
|
LORGLQMIN = MAX( 1, Q ) |
|
LORGLQOPT = INT( WORK(1) ) |
|
CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, |
CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, |
$ THETA, 0, 0, 1, V1T, LDV1T, U2, LDU2, U1, LDU1, |
$ THETA, DUM1, DUM2, 1, V1T, LDV1T, U2, |
$ 0, 0, 0, 0, 0, 0, 0, 0, WORK(1), -1, |
$ LDU2, U1, LDU1, DUM1, DUM1, DUM1, |
$ CHILDINFO ) |
$ DUM1, DUM1, DUM1, DUM1, |
|
$ DUM1, WORK(1), -1, CHILDINFO ) |
LBBCSD = INT( WORK(1) ) |
LBBCSD = INT( WORK(1) ) |
ELSE |
ELSE |
CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, |
CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, |
$ 0, 0, 0, WORK(1), -1, CHILDINFO ) |
$ DUM1, DUM1, DUM1, DUM1, |
|
$ DUM1, WORK(1), -1, CHILDINFO ) |
LORBDB = M + INT( WORK(1) ) |
LORBDB = M + INT( WORK(1) ) |
IF( P .GE. M-P ) THEN |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
CALL DORGQR( P, P, M-Q, U1, LDU1, 0, WORK(1), -1, |
CALL DORGQR( P, P, M-Q, U1, LDU1, DUM1, WORK(1), -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGQRMIN = MAX( 1, P ) |
LORGQRMIN = MAX( LORGQRMIN, P ) |
LORGQROPT = INT( WORK(1) ) |
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
ELSE |
END IF |
CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, 0, WORK(1), -1, |
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
|
CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, DUM1, WORK(1), |
|
$ -1, CHILDINFO ) |
|
LORGQRMIN = MAX( LORGQRMIN, M-P ) |
|
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) |
|
END IF |
|
IF( WANTV1T .AND. Q .GT. 0 ) THEN |
|
CALL DORGLQ( Q, Q, Q, V1T, LDV1T, DUM1, WORK(1), -1, |
$ CHILDINFO ) |
$ CHILDINFO ) |
LORGQRMIN = MAX( 1, M-P ) |
LORGLQMIN = MAX( LORGLQMIN, Q ) |
LORGQROPT = INT( WORK(1) ) |
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) |
END IF |
END IF |
CALL DORGLQ( Q, Q, Q, V1T, LDV1T, 0, WORK(1), -1, |
|
$ CHILDINFO ) |
|
LORGLQMIN = MAX( 1, Q ) |
|
LORGLQOPT = INT( WORK(1) ) |
|
CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, |
CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, |
$ THETA, 0, U2, LDU2, U1, LDU1, 0, 1, V1T, LDV1T, |
$ THETA, DUM1, U2, LDU2, U1, LDU1, DUM2, |
$ 0, 0, 0, 0, 0, 0, 0, 0, WORK(1), -1, |
$ 1, V1T, LDV1T, DUM1, DUM1, DUM1, |
$ CHILDINFO ) |
$ DUM1, DUM1, DUM1, DUM1, |
|
$ DUM1, WORK(1), -1, CHILDINFO ) |
LBBCSD = INT( WORK(1) ) |
LBBCSD = INT( WORK(1) ) |
END IF |
END IF |
LWORKMIN = MAX( IORBDB+LORBDB-1, |
LWORKMIN = MAX( IORBDB+LORBDB-1, |
Line 500
|
Line 522
|
CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), |
CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
END IF |
END IF |
* |
* |
* Simultaneously diagonalize X11 and X21. |
* Simultaneously diagonalize X11 and X21. |
* |
* |
CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, |
CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, |
$ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, 0, 1, |
$ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, |
$ WORK(IB11D), WORK(IB11E), WORK(IB12D), |
$ DUM2, 1, WORK(IB11D), WORK(IB11E), |
$ WORK(IB12E), WORK(IB21D), WORK(IB21E), |
$ WORK(IB12D), WORK(IB12E), WORK(IB21D), |
$ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD, |
$ WORK(IB21E), WORK(IB22D), WORK(IB22E), |
$ CHILDINFO ) |
$ WORK(IBBCSD), LBBCSD, CHILDINFO ) |
* |
* |
* Permute rows and columns to place zero submatrices in |
* Permute rows and columns to place zero submatrices in |
* preferred positions |
* preferred positions |
* |
* |
Line 554
|
Line 576
|
CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), |
CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
END IF |
END IF |
* |
* |
* Simultaneously diagonalize X11 and X21. |
* Simultaneously diagonalize X11 and X21. |
* |
* |
CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, |
CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, |
$ WORK(IPHI), V1T, LDV1T, 0, 1, U1, LDU1, U2, LDU2, |
$ WORK(IPHI), V1T, LDV1T, DUM1, 1, U1, LDU1, U2, |
$ WORK(IB11D), WORK(IB11E), WORK(IB12D), |
$ LDU2, WORK(IB11D), WORK(IB11E), WORK(IB12D), |
$ WORK(IB12E), WORK(IB21D), WORK(IB21E), |
$ WORK(IB12E), WORK(IB21D), WORK(IB21E), |
$ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD, |
$ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD, |
$ CHILDINFO ) |
$ CHILDINFO ) |
* |
* |
* Permute rows and columns to place identity submatrices in |
* Permute rows and columns to place identity submatrices in |
* preferred positions |
* preferred positions |
* |
* |
Line 609
|
Line 631
|
CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), |
CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
END IF |
END IF |
* |
* |
* Simultaneously diagonalize X11 and X21. |
* Simultaneously diagonalize X11 and X21. |
* |
* |
CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, |
CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, |
$ THETA, WORK(IPHI), 0, 1, V1T, LDV1T, U2, LDU2, U1, |
$ THETA, WORK(IPHI), DUM1, 1, V1T, LDV1T, U2, |
$ LDU1, WORK(IB11D), WORK(IB11E), WORK(IB12D), |
$ LDU2, U1, LDU1, WORK(IB11D), WORK(IB11E), |
$ WORK(IB12E), WORK(IB21D), WORK(IB21E), |
$ WORK(IB12D), WORK(IB12E), WORK(IB21D), |
$ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD, |
$ WORK(IB21E), WORK(IB22D), WORK(IB22E), |
$ CHILDINFO ) |
$ WORK(IBBCSD), LBBCSD, CHILDINFO ) |
* |
* |
* Permute rows and columns to place identity submatrices in |
* Permute rows and columns to place identity submatrices in |
* preferred positions |
* preferred positions |
* |
* |
Line 649
|
Line 671
|
* |
* |
* Accumulate Householder reflectors |
* Accumulate Householder reflectors |
* |
* |
|
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
|
CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 ) |
|
END IF |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
IF( WANTU1 .AND. P .GT. 0 ) THEN |
CALL DCOPY( P, WORK(IORBDB), 1, U1, 1 ) |
CALL DCOPY( P, WORK(IORBDB), 1, U1, 1 ) |
DO J = 2, P |
DO J = 2, P |
Line 660
|
Line 685
|
$ WORK(IORGQR), LORGQR, CHILDINFO ) |
$ WORK(IORGQR), LORGQR, CHILDINFO ) |
END IF |
END IF |
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
IF( WANTU2 .AND. M-P .GT. 0 ) THEN |
CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 ) |
|
DO J = 2, M-P |
DO J = 2, M-P |
U2(1,J) = ZERO |
U2(1,J) = ZERO |
END DO |
END DO |
Line 678
|
Line 702
|
CALL DORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1), |
CALL DORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1), |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
$ WORK(IORGLQ), LORGLQ, CHILDINFO ) |
END IF |
END IF |
* |
* |
* Simultaneously diagonalize X11 and X21. |
* Simultaneously diagonalize X11 and X21. |
* |
* |
CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, |
CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, |
$ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, 0, 1, V1T, |
$ THETA, WORK(IPHI), U2, LDU2, U1, LDU1, DUM1, |
$ LDV1T, WORK(IB11D), WORK(IB11E), WORK(IB12D), |
$ 1, V1T, LDV1T, WORK(IB11D), WORK(IB11E), |
$ WORK(IB12E), WORK(IB21D), WORK(IB21E), |
$ WORK(IB12D), WORK(IB12E), WORK(IB21D), |
$ WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD, |
$ WORK(IB21E), WORK(IB22D), WORK(IB22E), |
$ CHILDINFO ) |
$ WORK(IBBCSD), LBBCSD, CHILDINFO ) |
* |
* |
* Permute rows and columns to place identity submatrices in |
* Permute rows and columns to place identity submatrices in |
* preferred positions |
* preferred positions |
* |
* |