--- rpl/lapack/lapack/zuncsd2by1.f 2014/01/27 09:24:37 1.1 +++ rpl/lapack/lapack/zuncsd2by1.f 2023/08/07 08:39:43 1.10 @@ -2,8 +2,8 @@ * * =========== DOCUMENTATION =========== * -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download ZUNCSD2BY1 + dependencies @@ -22,7 +22,7 @@ * X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, * LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, * INFO ) -* +* * .. Scalar Arguments .. * CHARACTER JOBU1, JOBU2, JOBV1T * INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21, @@ -36,10 +36,10 @@ * $ X11(LDX11,*), X21(LDX21,*) * INTEGER IWORK(*) * .. -* -* +* +* *> \par Purpose: -*> ============= +* ============= *> *>\verbatim *> @@ -47,20 +47,20 @@ *> orthonormal columns that has been partitioned into a 2-by-1 block *> structure: *> -*> [ I 0 0 ] +*> [ I1 0 0 ] *> [ 0 C 0 ] *> [ X11 ] [ U1 | ] [ 0 0 0 ] *> X = [-----] = [---------] [----------] V1**T . *> [ X21 ] [ | U2 ] [ 0 0 0 ] *> [ 0 S 0 ] -*> [ 0 0 I ] -*> -*> 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). +*> [ 0 0 I2] *> -*>\endverbatim +*> X11 is P-by-Q. The unitary 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: * ========== @@ -68,172 +68,167 @@ *> \param[in] JOBU1 *> \verbatim *> JOBU1 is CHARACTER -*> = 'Y': U1 is computed; -*> otherwise: U1 is not computed. +*> = '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. +*> = '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. +*> = 'Y': V1T is computed; +*> otherwise: V1T is not computed. *> \endverbatim *> *> \param[in] M *> \verbatim *> M is INTEGER -*> The number of rows and columns in X. +*> The number of rows in X. *> \endverbatim *> *> \param[in] P *> \verbatim *> P is INTEGER -*> The number of rows in X11 and X12. 0 <= P <= M. +*> The number of rows in X11. 0 <= P <= M. *> \endverbatim *> *> \param[in] Q *> \verbatim *> 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 *> *> \param[in,out] X11 *> \verbatim *> X11 is COMPLEX*16 array, dimension (LDX11,Q) -*> On entry, part of the unitary matrix whose CSD is -*> desired. +*> 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). +*> The leading dimension of X11. LDX11 >= 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. +*> On entry, part of the unitary matrix whose CSD is desired. *> \endverbatim *> *> \param[in] LDX21 *> \verbatim *> LDX21 is INTEGER -*> The leading dimension of X21. LDX21 >= MAX(1,M-P). +*> The leading dimension of X21. LDX21 >= MAX(1,M-P). *> \endverbatim *> *> \param[out] THETA *> \verbatim -*> THETA is COMPLEX*16 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)) ). +*> 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. +*> 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). +*> 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. +*> 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). +*> 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**T. +*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary +*> matrix V1**T. *> \endverbatim *> *> \param[in] LDV1T *> \verbatim *> LDV1T is INTEGER -*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= -*> MAX(1,Q). +*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >= +*> MAX(1,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. -*> If INFO > 0 on exit, WORK(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. +*> 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. -*> \endverbatim -*> \verbatim -*> 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. +*> 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 and RWORK +*> arrays, returns this value as the first entry of the WORK +*> and RWORK array, respectively, and no error message related +*> to LWORK or LRWORK 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. +*> 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. +*> The dimension of the array RWORK. +*> +*> If LRWORK=-1, then a workspace query is assumed; the routine +*> only calculates the optimal size of the WORK and RWORK +*> arrays, returns this value as the first entry of the WORK +*> and RWORK array, respectively, and no error message related +*> to LWORK or LRWORK is issued by XERBLA. +*> \endverbatim +* *> \param[out] IWORK *> \verbatim *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q)) *> \endverbatim -*> \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 WORK +*> = 0: successful exit. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> > 0: ZBBCSD did not converge. See the description of WORK *> above for details. *> \endverbatim * @@ -246,12 +241,10 @@ * Authors: * ======== * -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \date July 2012 +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. * *> \ingroup complex16OTHERcomputational * @@ -261,10 +254,9 @@ $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, $ INFO ) * -* -- LAPACK computational routine (version 3.5.0) -- +* -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* July 2012 * * .. Scalar Arguments .. CHARACTER JOBU1, JOBU2, JOBV1T @@ -279,7 +271,7 @@ $ X11(LDX11,*), X21(LDX21,*) INTEGER IWORK(*) * .. -* +* * ===================================================================== * * .. Parameters .. @@ -295,6 +287,10 @@ $ LWORKMIN, LWORKOPT, R LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T * .. +* .. Local Arrays .. + DOUBLE PRECISION DUM( 1 ) + COMPLEX*16 CDUM( 1, 1 ) +* .. * .. External Subroutines .. EXTERNAL ZBBCSD, ZCOPY, ZLACPY, ZLAPMR, ZLAPMT, ZUNBDB1, $ ZUNBDB2, ZUNBDB3, ZUNBDB4, ZUNGLQ, ZUNGQR, @@ -315,7 +311,7 @@ WANTU1 = LSAME( JOBU1, 'Y' ) WANTU2 = LSAME( JOBU2, 'Y' ) WANTV1T = LSAME( JOBV1T, 'Y' ) - LQUERY = LWORK .EQ. -1 + LQUERY = ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 ) * IF( M .LT. 0 ) THEN INFO = -4 @@ -327,11 +323,11 @@ INFO = -8 ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN INFO = -10 - ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN + ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN INFO = -13 - ELSE IF( WANTU2 .AND. LDU2 .LT. M - P ) THEN + ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN INFO = -15 - ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN + ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN INFO = -17 END IF * @@ -387,99 +383,118 @@ IORBDB = ITAUQ1 + MAX( 1, Q ) IORGQR = ITAUQ1 + MAX( 1, Q ) IORGLQ = ITAUQ1 + MAX( 1, Q ) + LORGQRMIN = 1 + LORGQROPT = 1 + LORGLQMIN = 1 + LORGLQOPT = 1 IF( R .EQ. Q ) THEN - CALL ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, - $ 0, 0, WORK, -1, CHILDINFO ) + CALL ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + $ CDUM, CDUM, CDUM, WORK, -1, CHILDINFO ) LORBDB = INT( WORK(1) ) - IF( P .GE. M-P ) THEN - CALL ZUNGQR( P, P, Q, U1, LDU1, 0, WORK(1), -1, + IF( WANTU1 .AND. P .GT. 0 ) THEN + CALL ZUNGQR( P, P, Q, U1, LDU1, CDUM, WORK(1), -1, $ CHILDINFO ) - LORGQRMIN = MAX( 1, P ) - LORGQROPT = INT( WORK(1) ) - ELSE - CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, 0, WORK(1), -1, + LORGQRMIN = MAX( LORGQRMIN, P ) + LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) + ENDIF + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, CDUM, WORK(1), -1, $ CHILDINFO ) - LORGQRMIN = MAX( 1, M-P ) - 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 ZUNGLQ( Q-1, Q-1, Q-1, V1T, LDV1T, + $ CDUM, WORK(1), -1, CHILDINFO ) + LORGLQMIN = MAX( LORGLQMIN, Q-1 ) + LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL ZUNGLQ( 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 ZBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, - $ 0, U1, LDU1, U2, LDU2, V1T, LDV1T, 0, 1, 0, 0, - $ 0, 0, 0, 0, 0, 0, RWORK(1), -1, CHILDINFO ) + $ DUM, U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM, 1, + $ DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, + $ RWORK(1), -1, CHILDINFO ) LBBCSD = INT( RWORK(1) ) ELSE IF( R .EQ. P ) THEN - CALL ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, - $ 0, 0, WORK(1), -1, CHILDINFO ) + CALL ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) LORBDB = INT( WORK(1) ) - IF( P-1 .GE. M-P ) THEN - CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, 0, WORK(1), + IF( WANTU1 .AND. P .GT. 0 ) THEN + CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, WORK(1), $ -1, CHILDINFO ) - LORGQRMIN = MAX( 1, P-1 ) - LORGQROPT = INT( WORK(1) ) - ELSE - CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, 0, WORK(1), -1, + LORGQRMIN = MAX( LORGQRMIN, P-1 ) + LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) + END IF + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL ZUNGQR( M-P, M-P, Q, U2, LDU2, CDUM, WORK(1), -1, $ CHILDINFO ) - LORGQRMIN = MAX( 1, M-P ) - 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 ZUNGLQ( Q, Q, R, V1T, LDV1T, CDUM, WORK(1), -1, + $ CHILDINFO ) + LORGLQMIN = MAX( LORGLQMIN, Q ) + LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, 0, WORK(1), -1, - $ CHILDINFO ) - LORGLQMIN = MAX( 1, Q ) - LORGLQOPT = INT( WORK(1) ) CALL ZBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, - $ 0, V1T, LDV1T, 0, 1, U1, LDU1, U2, LDU2, 0, 0, - $ 0, 0, 0, 0, 0, 0, RWORK(1), -1, CHILDINFO ) + $ DUM, V1T, LDV1T, CDUM, 1, U1, LDU1, U2, LDU2, + $ DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, + $ RWORK(1), -1, CHILDINFO ) LBBCSD = INT( RWORK(1) ) ELSE IF( R .EQ. M-P ) THEN - CALL ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, - $ 0, 0, WORK(1), -1, CHILDINFO ) + CALL ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) LORBDB = INT( WORK(1) ) - IF( P .GE. M-P-1 ) THEN - CALL ZUNGQR( P, P, Q, U1, LDU1, 0, WORK(1), -1, + IF( WANTU1 .AND. P .GT. 0 ) THEN + CALL ZUNGQR( P, P, Q, U1, LDU1, CDUM, WORK(1), -1, $ CHILDINFO ) - LORGQRMIN = MAX( 1, P ) - LORGQROPT = INT( WORK(1) ) - ELSE - CALL ZUNGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, 0, + LORGQRMIN = MAX( LORGQRMIN, P ) + LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) + END IF + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL ZUNGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, CDUM, $ WORK(1), -1, CHILDINFO ) - LORGQRMIN = MAX( 1, M-P-1 ) - 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 ZUNGLQ( Q, Q, R, V1T, LDV1T, CDUM, WORK(1), -1, + $ CHILDINFO ) + LORGLQMIN = MAX( LORGLQMIN, Q ) + LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, 0, WORK(1), -1, - $ CHILDINFO ) - LORGLQMIN = MAX( 1, Q ) - LORGLQOPT = INT( WORK(1) ) CALL ZBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, - $ THETA, 0, 0, 1, V1T, LDV1T, U2, LDU2, U1, LDU1, - $ 0, 0, 0, 0, 0, 0, 0, 0, RWORK(1), -1, - $ CHILDINFO ) + $ THETA, DUM, CDUM, 1, V1T, LDV1T, U2, LDU2, U1, + $ LDU1, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, + $ RWORK(1), -1, CHILDINFO ) LBBCSD = INT( RWORK(1) ) ELSE - CALL ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, 0, 0, - $ 0, 0, 0, WORK(1), -1, CHILDINFO ) + CALL ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + $ CDUM, CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO + $ ) LORBDB = M + INT( WORK(1) ) - IF( P .GE. M-P ) THEN - CALL ZUNGQR( P, P, M-Q, U1, LDU1, 0, WORK(1), -1, + IF( WANTU1 .AND. P .GT. 0 ) THEN + CALL ZUNGQR( P, P, M-Q, U1, LDU1, CDUM, WORK(1), -1, + $ CHILDINFO ) + LORGQRMIN = MAX( LORGQRMIN, P ) + LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) + END IF + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL ZUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), -1, $ CHILDINFO ) - LORGQRMIN = MAX( 1, P ) - LORGQROPT = INT( WORK(1) ) - ELSE - CALL ZUNGQR( M-P, M-P, M-Q, U2, LDU2, 0, WORK(1), -1, + LORGQRMIN = MAX( LORGQRMIN, M-P ) + LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) + END IF + IF( WANTV1T .AND. Q .GT. 0 ) THEN + CALL ZUNGLQ( Q, Q, Q, V1T, LDV1T, CDUM, WORK(1), -1, $ CHILDINFO ) - LORGQRMIN = MAX( 1, M-P ) - LORGQROPT = INT( WORK(1) ) + LORGLQMIN = MAX( LORGLQMIN, Q ) + LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL ZUNGLQ( Q, Q, Q, V1T, LDV1T, 0, WORK(1), -1, - $ CHILDINFO ) - LORGLQMIN = MAX( 1, Q ) - LORGLQOPT = INT( WORK(1) ) CALL ZBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, - $ THETA, 0, U2, LDU2, U1, LDU1, 0, 1, V1T, LDV1T, - $ 0, 0, 0, 0, 0, 0, 0, 0, RWORK(1), -1, - $ CHILDINFO ) + $ THETA, DUM, U2, LDU2, U1, LDU1, CDUM, 1, V1T, + $ LDV1T, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM, + $ RWORK(1), -1, CHILDINFO ) LBBCSD = INT( RWORK(1) ) END IF LRWORKMIN = IBBCSD+LBBCSD-1 @@ -495,6 +510,9 @@ IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN INFO = -19 END IF + IF( LRWORK .LT. LRWORKMIN .AND. .NOT.LQUERY ) THEN + INFO = -21 + END IF END IF IF( INFO .NE. 0 ) THEN CALL XERBLA( 'ZUNCSD2BY1', -INFO ) @@ -541,16 +559,16 @@ CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF -* +* * Simultaneously diagonalize X11 and X21. -* +* CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, - $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, 0, 1, - $ RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), + $ RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM, + $ 1, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), - $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), LBBCSD, - $ CHILDINFO ) -* + $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), + $ LRWORK-IBBCSD+1, CHILDINFO ) +* * Permute rows and columns to place zero submatrices in * preferred positions * @@ -595,16 +613,16 @@ CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF -* +* * Simultaneously diagonalize X11 and X21. -* +* CALL ZBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, - $ RWORK(IPHI), V1T, LDV1T, 0, 1, U1, LDU1, U2, LDU2, - $ RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), + $ RWORK(IPHI), V1T, LDV1T, CDUM, 1, U1, LDU1, U2, + $ LDU2, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), LBBCSD, $ CHILDINFO ) -* +* * Permute rows and columns to place identity submatrices in * preferred positions * @@ -650,16 +668,16 @@ CALL ZUNGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF -* +* * Simultaneously diagonalize X11 and X21. -* +* CALL ZBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P, - $ THETA, RWORK(IPHI), 0, 1, V1T, LDV1T, U2, LDU2, + $ THETA, RWORK(IPHI), CDUM, 1, V1T, LDV1T, U2, LDU2, $ U1, LDU1, RWORK(IB11D), RWORK(IB11E), $ RWORK(IB12D), RWORK(IB12E), RWORK(IB21D), $ RWORK(IB21E), RWORK(IB22D), RWORK(IB22E), $ RWORK(IBBCSD), LBBCSD, CHILDINFO ) -* +* * Permute rows and columns to place identity submatrices in * preferred positions * @@ -690,6 +708,9 @@ * * Accumulate Householder reflectors * + IF( WANTU2 .AND. M-P .GT. 0 ) THEN + CALL ZCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 ) + END IF IF( WANTU1 .AND. P .GT. 0 ) THEN CALL ZCOPY( P, WORK(IORBDB), 1, U1, 1 ) DO J = 2, P @@ -701,7 +722,6 @@ $ WORK(IORGQR), LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN - CALL ZCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 ) DO J = 2, M-P U2(1,J) = ZERO END DO @@ -719,16 +739,16 @@ CALL ZUNGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF -* +* * Simultaneously diagonalize X11 and X21. -* +* CALL ZBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q, - $ THETA, RWORK(IPHI), U2, LDU2, U1, LDU1, 0, 1, V1T, - $ LDV1T, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D), - $ RWORK(IB12E), RWORK(IB21D), RWORK(IB21E), - $ RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), LBBCSD, - $ CHILDINFO ) -* + $ THETA, RWORK(IPHI), U2, LDU2, U1, LDU1, CDUM, 1, + $ V1T, LDV1T, RWORK(IB11D), RWORK(IB11E), + $ RWORK(IB12D), RWORK(IB12E), RWORK(IB21D), + $ RWORK(IB21E), RWORK(IB22D), RWORK(IB22E), + $ RWORK(IBBCSD), LBBCSD, CHILDINFO ) +* * Permute rows and columns to place identity submatrices in * preferred positions *