Skip to content

Handle degenerate cases in OR/UNCSD2BY1 (Reference-LAPACK PR 1291) #5860

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
martin-frbg merged 1 commit into from
Jun 25, 2026
+227 −60
Merged

Handle degenerate cases in OR/UNCSD2BY1 (Reference-LAPACK PR 1291) #5860

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 61 additions & 17 deletions lapack-netlib/SRC/cuncsd2by1.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,13 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CUNCSD2BY1 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd2by1.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd2by1.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd2by1.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
Expand Down Expand Up @@ -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, --
Expand Down Expand Up @@ -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
* ..
Expand Down Expand Up @@ -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) ) )
Expand All @@ -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
Expand All @@ -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) )
Expand All @@ -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) ) )
Expand All @@ -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 )
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
*
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 )
Expand Down
63 changes: 51 additions & 12 deletions lapack-netlib/SRC/dorcsd2by1.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,13 @@
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DORCSD2BY1 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
Expand Down Expand Up @@ -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, --
Expand Down Expand Up @@ -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
* ..
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
*
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 )
Expand Down
Loading
Loading