From 9e9c5a541f27c23ce3764d6852aeeaf3d8ebdd0c Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Wed, 24 Jun 2026 20:35:34 +0200 Subject: [PATCH] Handle degenerate cases having R=0 (Reference-LAPACK PR 1291) --- lapack-netlib/SRC/cuncsd2by1.f | 78 +++++++++++++++++++++++++------- lapack-netlib/SRC/dorcsd2by1.f | 63 +++++++++++++++++++++----- lapack-netlib/SRC/sorcsd2by1.f | 65 +++++++++++++++++++++------ lapack-netlib/SRC/zuncsd2by1.f | 81 ++++++++++++++++++++++++++-------- 4 files changed, 227 insertions(+), 60 deletions(-) diff --git a/lapack-netlib/SRC/cuncsd2by1.f b/lapack-netlib/SRC/cuncsd2by1.f index 128e82cecf..f1f42379a0 100644 --- a/lapack-netlib/SRC/cuncsd2by1.f +++ b/lapack-netlib/SRC/cuncsd2by1.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CUNCSD2BY1 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -250,10 +248,12 @@ *> \ingroup uncsd2by1 * * ===================================================================== - SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, + SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, + $ LDX11, $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -293,7 +293,9 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, COMPLEX CDUM( 1, 1 ) * .. * .. External Subroutines .. - EXTERNAL CBBCSD, CCOPY, CLACPY, CLAPMR, CLAPMT, CUNBDB1, + EXTERNAL CBBCSD, CCOPY, CLACPY, CLAPMR, CLAPMT, + $ CLASET, + $ CUNBDB1, $ CUNBDB2, CUNBDB3, CUNBDB4, CUNGLQ, CUNGQR, $ XERBLA * .. @@ -412,17 +414,20 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q-1 ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL CBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, + CALL CBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, + $ THETA, $ DUM(1), 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 CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + CALL CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) LORBDB = INT( WORK(1) ) IF( WANTU1 .AND. P .GT. 0 ) THEN - CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, WORK(1), + CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, + $ WORK(1), $ -1, CHILDINFO ) LORGQRMIN = MAX( LORGQRMIN, P-1 ) LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) @@ -439,13 +444,15 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL CBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, + CALL CBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, + $ THETA, $ 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 CUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + CALL CUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) LORBDB = INT( WORK(1) ) IF( WANTU1 .AND. P .GT. 0 ) THEN @@ -472,7 +479,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, $ RWORK(1), -1, CHILDINFO ) LBBCSD = INT( RWORK(1) ) ELSE - CALL CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + CALL CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO $ ) LORBDB = M + INT( WORK(1) ) @@ -483,7 +491,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN - CALL CUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), -1, + CALL CUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), + $ -1, $ CHILDINFO ) LORGQRMIN = MAX( LORGQRMIN, M-P ) LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) @@ -502,7 +511,7 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF LRWORKMIN = IBBCSD+LBBCSD-1 LRWORKOPT = LRWORKMIN - RWORK(1) = LRWORKOPT + RWORK(1) = REAL( LRWORKOPT ) LWORKMIN = MAX( IORBDB+LORBDB-1, $ IORGQR+LORGQRMIN-1, $ IORGLQ+LORGLQMIN-1 ) @@ -525,6 +534,36 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF LORGQR = LWORK-IORGQR+1 LORGLQ = LWORK-IORGLQ+1 +* + IF( R .EQ. 0 ) THEN +* +* R = 0: C and S are empty. Handle the trivial CSD directly. +* + IF( Q .EQ. 0 ) THEN + IF( WANTU1 .AND. P .GT. 0 ) + $ CALL CLASET( 'A', P, P, ZERO, ONE, U1, LDU1 ) + IF( WANTU2 .AND. M-P .GT. 0 ) + $ CALL CLASET( 'A', M-P, M-P, ZERO, ONE, U2, LDU2 ) + RETURN + END IF +* + IF( P .EQ. 0 .AND. M .EQ. Q ) THEN + IF( WANTU2 ) + $ CALL CLACPY( 'A', M-P, Q, X21, LDX21, U2, LDU2 ) + IF( WANTV1T ) + $ CALL CLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + IF( P .EQ. M .AND. M .EQ. Q ) THEN + IF( WANTU1 ) + $ CALL CLACPY( 'A', P, Q, X11, LDX11, U1, LDU1 ) + IF( WANTV1T ) + $ CALL CLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + END IF * * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q, * in which R = MIN(P,M-P,Q,M-Q) @@ -543,7 +582,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -559,7 +599,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END DO CALL CLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2), $ LDV1T ) - CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), + CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, + $ WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF * @@ -602,7 +643,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, U1(1,J) = ZERO U1(J,1) = ZERO END DO - CALL CLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 ) + CALL CLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), + $ LDU1 ) CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1), $ WORK(IORGQR), LORGQR, CHILDINFO ) END IF @@ -652,7 +694,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -736,7 +779,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF IF( WANTV1T .AND. Q .GT. 0 ) THEN CALL CLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T ) - CALL CLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11, + CALL CLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), + $ LDX11, $ V1T(M-Q+1,M-Q+1), LDV1T ) CALL CLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21, $ V1T(P+1,P+1), LDV1T ) diff --git a/lapack-netlib/SRC/dorcsd2by1.f b/lapack-netlib/SRC/dorcsd2by1.f index 25fab0f33c..a9020958c2 100644 --- a/lapack-netlib/SRC/dorcsd2by1.f +++ b/lapack-netlib/SRC/dorcsd2by1.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DORCSD2BY1 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -224,12 +222,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHERcomputational +*> \ingroup uncsd2by1 * * ===================================================================== - 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, $ LDV1T, WORK, LWORK, IWORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine (3.5.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -266,7 +266,9 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, DOUBLE PRECISION DUM1(1), DUM2(1,1) * .. * .. External Subroutines .. - EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1, + EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, + $ DLASET, + $ DORBDB1, $ DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR, $ XERBLA * .. @@ -370,7 +372,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q-1 ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, + CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, + $ THETA, $ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T, $ DUM2, 1, DUM1, DUM1, DUM1, $ DUM1, DUM1, DUM1, DUM1, @@ -399,7 +402,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, + CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, + $ THETA, $ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1, $ U2, LDU2, DUM1, DUM1, DUM1, $ DUM1, DUM1, DUM1, DUM1, @@ -485,6 +489,36 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF LORGQR = LWORK-IORGQR+1 LORGLQ = LWORK-IORGLQ+1 +* + IF( R .EQ. 0 ) THEN +* +* R = 0: C and S are empty. Handle the trivial CSD directly. +* + IF( Q .EQ. 0 ) THEN + IF( WANTU1 .AND. P .GT. 0 ) + $ CALL DLASET( 'A', P, P, ZERO, ONE, U1, LDU1 ) + IF( WANTU2 .AND. M-P .GT. 0 ) + $ CALL DLASET( 'A', M-P, M-P, ZERO, ONE, U2, LDU2 ) + RETURN + END IF +* + IF( P .EQ. 0 .AND. M .EQ. Q ) THEN + IF( WANTU2 ) + $ CALL DLACPY( 'A', M-P, Q, X21, LDX21, U2, LDU2 ) + IF( WANTV1T ) + $ CALL DLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + IF( P .EQ. M .AND. M .EQ. Q ) THEN + IF( WANTU1 ) + $ CALL DLACPY( 'A', P, Q, X11, LDX11, U1, LDU1 ) + IF( WANTV1T ) + $ CALL DLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + END IF * * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q, * in which R = MIN(P,M-P,Q,M-Q) @@ -503,7 +537,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -519,7 +554,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END DO CALL DLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2), $ LDV1T ) - 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 ) END IF * @@ -562,7 +598,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, U1(1,J) = ZERO U1(J,1) = ZERO END DO - CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 ) + CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), + $ LDU1 ) CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1), $ WORK(IORGQR), LORGQR, CHILDINFO ) END IF @@ -612,7 +649,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -695,7 +733,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF IF( WANTV1T .AND. Q .GT. 0 ) THEN CALL DLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T ) - CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11, + CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), + $ LDX11, $ V1T(M-Q+1,M-Q+1), LDV1T ) CALL DLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21, $ V1T(P+1,P+1), LDV1T ) diff --git a/lapack-netlib/SRC/sorcsd2by1.f b/lapack-netlib/SRC/sorcsd2by1.f index 25c317f6f6..b5eb82f906 100644 --- a/lapack-netlib/SRC/sorcsd2by1.f +++ b/lapack-netlib/SRC/sorcsd2by1.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SORCSD2BY1 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -224,12 +222,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realOTHERcomputational +*> \ingroup uncsd2by1 * * ===================================================================== - SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, + SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, + $ LDX11, $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, $ LDV1T, WORK, LWORK, IWORK, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -266,7 +266,9 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, REAL DUM1(1), DUM2(1,1) * .. * .. External Subroutines .. - EXTERNAL SBBCSD, SCOPY, SLACPY, SLAPMR, SLAPMT, SORBDB1, + EXTERNAL SBBCSD, SCOPY, SLACPY, SLAPMR, SLAPMT, + $ SLASET, + $ SORBDB1, $ SORBDB2, SORBDB3, SORBDB4, SORGLQ, SORGQR, $ XERBLA * .. @@ -370,7 +372,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q-1 ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL SBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, + CALL SBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, + $ THETA, $ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T, DUM2, $ 1, DUM1, DUM1, DUM1, DUM1, DUM1, $ DUM1, DUM1, DUM1, WORK(1), -1, CHILDINFO @@ -399,7 +402,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL SBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, + CALL SBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, + $ THETA, $ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1, U2, $ LDU2, DUM1, DUM1, DUM1, DUM1, DUM1, $ DUM1, DUM1, DUM1, WORK(1), -1, CHILDINFO @@ -472,7 +476,7 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, $ IORGQR+LORGQROPT-1, $ IORGLQ+LORGLQOPT-1, $ IBBCSD+LBBCSD-1 ) - WORK(1) = LWORKOPT + WORK(1) = REAL( LWORKOPT ) IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN INFO = -19 END IF @@ -485,6 +489,36 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF LORGQR = LWORK-IORGQR+1 LORGLQ = LWORK-IORGLQ+1 +* + IF( R .EQ. 0 ) THEN +* +* R = 0: C and S are empty. Handle the trivial CSD directly. +* + IF( Q .EQ. 0 ) THEN + IF( WANTU1 .AND. P .GT. 0 ) + $ CALL SLASET( 'A', P, P, ZERO, ONE, U1, LDU1 ) + IF( WANTU2 .AND. M-P .GT. 0 ) + $ CALL SLASET( 'A', M-P, M-P, ZERO, ONE, U2, LDU2 ) + RETURN + END IF +* + IF( P .EQ. 0 .AND. M .EQ. Q ) THEN + IF( WANTU2 ) + $ CALL SLACPY( 'A', M-P, Q, X21, LDX21, U2, LDU2 ) + IF( WANTV1T ) + $ CALL SLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + IF( P .EQ. M .AND. M .EQ. Q ) THEN + IF( WANTU1 ) + $ CALL SLACPY( 'A', P, Q, X11, LDX11, U1, LDU1 ) + IF( WANTV1T ) + $ CALL SLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + END IF * * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q, * in which R = MIN(P,M-P,Q,M-Q) @@ -503,7 +537,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL SLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -519,7 +554,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END DO CALL SLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2), $ LDV1T ) - CALL SORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), + CALL SORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, + $ WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF * @@ -562,7 +598,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, U1(1,J) = ZERO U1(J,1) = ZERO END DO - CALL SLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 ) + CALL SLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), + $ LDU1 ) CALL SORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1), $ WORK(IORGQR), LORGQR, CHILDINFO ) END IF @@ -612,7 +649,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL SLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -695,7 +733,8 @@ SUBROUTINE SORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF IF( WANTV1T .AND. Q .GT. 0 ) THEN CALL SLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T ) - CALL SLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11, + CALL SLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), + $ LDX11, $ V1T(M-Q+1,M-Q+1), LDV1T ) CALL SLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21, $ V1T(P+1,P+1), LDV1T ) diff --git a/lapack-netlib/SRC/zuncsd2by1.f b/lapack-netlib/SRC/zuncsd2by1.f index 399b598be0..1b28e691fb 100644 --- a/lapack-netlib/SRC/zuncsd2by1.f +++ b/lapack-netlib/SRC/zuncsd2by1.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZUNCSD2BY1 + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -246,13 +244,15 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERcomputational +*> \ingroup uncsd2by1 * * ===================================================================== - SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, + SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, + $ LDX11, $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -292,7 +292,9 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, COMPLEX*16 CDUM( 1, 1 ) * .. * .. External Subroutines .. - EXTERNAL ZBBCSD, ZCOPY, ZLACPY, ZLAPMR, ZLAPMT, ZUNBDB1, + EXTERNAL ZBBCSD, ZCOPY, ZLACPY, ZLAPMR, ZLAPMT, + $ ZLASET, + $ ZUNBDB1, $ ZUNBDB2, ZUNBDB3, ZUNBDB4, ZUNGLQ, ZUNGQR, $ XERBLA * .. @@ -388,7 +390,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = 1 LORGLQOPT = 1 IF( R .EQ. Q ) THEN - CALL ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + CALL ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, WORK, -1, CHILDINFO ) LORBDB = INT( WORK(1) ) IF( WANTU1 .AND. P .GT. 0 ) THEN @@ -409,17 +412,20 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q-1 ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA, + CALL ZBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, + $ THETA, $ 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, DUM, + CALL ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) LORBDB = INT( WORK(1) ) IF( WANTU1 .AND. P .GT. 0 ) THEN - CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, WORK(1), + CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, + $ WORK(1), $ -1, CHILDINFO ) LORGQRMIN = MAX( LORGQRMIN, P-1 ) LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) @@ -436,13 +442,15 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, LORGLQMIN = MAX( LORGLQMIN, Q ) LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) ) END IF - CALL ZBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA, + CALL ZBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, + $ THETA, $ 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, DUM, + CALL ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO ) LORBDB = INT( WORK(1) ) IF( WANTU1 .AND. P .GT. 0 ) THEN @@ -469,7 +477,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, $ RWORK(1), -1, CHILDINFO ) LBBCSD = INT( RWORK(1) ) ELSE - CALL ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM, + CALL ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, + $ DUM, $ CDUM, CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO $ ) LORBDB = M + INT( WORK(1) ) @@ -480,7 +489,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, 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, + CALL ZUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), + $ -1, $ CHILDINFO ) LORGQRMIN = MAX( LORGQRMIN, M-P ) LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) ) @@ -522,6 +532,36 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF LORGQR = LWORK-IORGQR+1 LORGLQ = LWORK-IORGLQ+1 +* + IF( R .EQ. 0 ) THEN +* +* R = 0: C and S are empty. Handle the trivial CSD directly. +* + IF( Q .EQ. 0 ) THEN + IF( WANTU1 .AND. P .GT. 0 ) + $ CALL ZLASET( 'A', P, P, ZERO, ONE, U1, LDU1 ) + IF( WANTU2 .AND. M-P .GT. 0 ) + $ CALL ZLASET( 'A', M-P, M-P, ZERO, ONE, U2, LDU2 ) + RETURN + END IF +* + IF( P .EQ. 0 .AND. M .EQ. Q ) THEN + IF( WANTU2 ) + $ CALL ZLACPY( 'A', M-P, Q, X21, LDX21, U2, LDU2 ) + IF( WANTV1T ) + $ CALL ZLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + IF( P .EQ. M .AND. M .EQ. Q ) THEN + IF( WANTU1 ) + $ CALL ZLACPY( 'A', P, Q, X11, LDX11, U1, LDU1 ) + IF( WANTV1T ) + $ CALL ZLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T ) + RETURN + END IF +* + END IF * * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q, * in which R = MIN(P,M-P,Q,M-Q) @@ -540,7 +580,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL ZLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -556,7 +597,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END DO CALL ZLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2), $ LDV1T ) - CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), + CALL ZUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, + $ WORK(ITAUQ1), $ WORK(IORGLQ), LORGLQ, CHILDINFO ) END IF * @@ -599,7 +641,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, U1(1,J) = ZERO U1(J,1) = ZERO END DO - CALL ZLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 ) + CALL ZLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), + $ LDU1 ) CALL ZUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1), $ WORK(IORGQR), LORGQR, CHILDINFO ) END IF @@ -649,7 +692,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, * IF( WANTU1 .AND. P .GT. 0 ) THEN CALL ZLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) - CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), + CALL ZUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), + $ WORK(IORGQR), $ LORGQR, CHILDINFO ) END IF IF( WANTU2 .AND. M-P .GT. 0 ) THEN @@ -732,7 +776,8 @@ SUBROUTINE ZUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, END IF IF( WANTV1T .AND. Q .GT. 0 ) THEN CALL ZLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T ) - CALL ZLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11, + CALL ZLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), + $ LDX11, $ V1T(M-Q+1,M-Q+1), LDV1T ) CALL ZLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21, $ V1T(P+1,P+1), LDV1T )