Skip to content

Commit 3cdf5dc

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5860 from martin-frbg/lapack1291
Handle degenerate cases in OR/UNCSD2BY1 (Reference-LAPACK PR 1291)
2 parents 56b0893 + 9e9c5a5 commit 3cdf5dc

4 files changed

Lines changed: 227 additions & 60 deletions

File tree

‎lapack-netlib/SRC/cuncsd2by1.f‎

Lines changed: 61 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download CUNCSD2BY1 + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd2by1.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd2by1.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd2by1.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -250,10 +248,12 @@
250248
*> \ingroup uncsd2by1
251249
*
252250
* =====================================================================
253-
SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
251+
SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11,
252+
$ LDX11,
254253
$ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
255254
$ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
256255
$ INFO )
256+
IMPLICIT NONE
257257
*
258258
* -- LAPACK computational routine --
259259
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -293,7 +293,9 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
293293
COMPLEX CDUM( 1, 1 )
294294
* ..
295295
* .. External Subroutines ..
296-
EXTERNAL CBBCSD, CCOPY, CLACPY, CLAPMR, CLAPMT, CUNBDB1,
296+
EXTERNAL CBBCSD, CCOPY, CLACPY, CLAPMR, CLAPMT,
297+
$ CLASET,
298+
$ CUNBDB1,
297299
$ CUNBDB2, CUNBDB3, CUNBDB4, CUNGLQ, CUNGQR,
298300
$ XERBLA
299301
* ..
@@ -412,17 +414,20 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
412414
LORGLQMIN = MAX( LORGLQMIN, Q-1 )
413415
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
414416
END IF
415-
CALL CBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
417+
CALL CBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q,
418+
$ THETA,
416419
$ DUM(1), U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM,
417420
$ 1, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM,
418421
$ RWORK(1), -1, CHILDINFO )
419422
LBBCSD = INT( RWORK(1) )
420423
ELSE IF( R .EQ. P ) THEN
421-
CALL CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM,
424+
CALL CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
425+
$ DUM,
422426
$ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO )
423427
LORBDB = INT( WORK(1) )
424428
IF( WANTU1 .AND. P .GT. 0 ) THEN
425-
CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, WORK(1),
429+
CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM,
430+
$ WORK(1),
426431
$ -1, CHILDINFO )
427432
LORGQRMIN = MAX( LORGQRMIN, P-1 )
428433
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
@@ -439,13 +444,15 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
439444
LORGLQMIN = MAX( LORGLQMIN, Q )
440445
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
441446
END IF
442-
CALL CBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
447+
CALL CBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P,
448+
$ THETA,
443449
$ DUM, V1T, LDV1T, CDUM, 1, U1, LDU1, U2, LDU2,
444450
$ DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM,
445451
$ RWORK(1), -1, CHILDINFO )
446452
LBBCSD = INT( RWORK(1) )
447453
ELSE IF( R .EQ. M-P ) THEN
448-
CALL CUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM,
454+
CALL CUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
455+
$ DUM,
449456
$ CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO )
450457
LORBDB = INT( WORK(1) )
451458
IF( WANTU1 .AND. P .GT. 0 ) THEN
@@ -472,7 +479,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
472479
$ RWORK(1), -1, CHILDINFO )
473480
LBBCSD = INT( RWORK(1) )
474481
ELSE
475-
CALL CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM,
482+
CALL CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
483+
$ DUM,
476484
$ CDUM, CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO
477485
$ )
478486
LORBDB = M + INT( WORK(1) )
@@ -483,7 +491,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
483491
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
484492
END IF
485493
IF( WANTU2 .AND. M-P .GT. 0 ) THEN
486-
CALL CUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), -1,
494+
CALL CUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1),
495+
$ -1,
487496
$ CHILDINFO )
488497
LORGQRMIN = MAX( LORGQRMIN, M-P )
489498
LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
@@ -502,7 +511,7 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
502511
END IF
503512
LRWORKMIN = IBBCSD+LBBCSD-1
504513
LRWORKOPT = LRWORKMIN
505-
RWORK(1) = LRWORKOPT
514+
RWORK(1) = REAL( LRWORKOPT )
506515
LWORKMIN = MAX( IORBDB+LORBDB-1,
507516
$ IORGQR+LORGQRMIN-1,
508517
$ IORGLQ+LORGLQMIN-1 )
@@ -525,6 +534,36 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
525534
END IF
526535
LORGQR = LWORK-IORGQR+1
527536
LORGLQ = LWORK-IORGLQ+1
537+
*
538+
IF( R .EQ. 0 ) THEN
539+
*
540+
* R = 0: C and S are empty. Handle the trivial CSD directly.
541+
*
542+
IF( Q .EQ. 0 ) THEN
543+
IF( WANTU1 .AND. P .GT. 0 )
544+
$ CALL CLASET( 'A', P, P, ZERO, ONE, U1, LDU1 )
545+
IF( WANTU2 .AND. M-P .GT. 0 )
546+
$ CALL CLASET( 'A', M-P, M-P, ZERO, ONE, U2, LDU2 )
547+
RETURN
548+
END IF
549+
*
550+
IF( P .EQ. 0 .AND. M .EQ. Q ) THEN
551+
IF( WANTU2 )
552+
$ CALL CLACPY( 'A', M-P, Q, X21, LDX21, U2, LDU2 )
553+
IF( WANTV1T )
554+
$ CALL CLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T )
555+
RETURN
556+
END IF
557+
*
558+
IF( P .EQ. M .AND. M .EQ. Q ) THEN
559+
IF( WANTU1 )
560+
$ CALL CLACPY( 'A', P, Q, X11, LDX11, U1, LDU1 )
561+
IF( WANTV1T )
562+
$ CALL CLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T )
563+
RETURN
564+
END IF
565+
*
566+
END IF
528567
*
529568
* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
530569
* in which R = MIN(P,M-P,Q,M-Q)
@@ -543,7 +582,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
543582
*
544583
IF( WANTU1 .AND. P .GT. 0 ) THEN
545584
CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
546-
CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
585+
CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1),
586+
$ WORK(IORGQR),
547587
$ LORGQR, CHILDINFO )
548588
END IF
549589
IF( WANTU2 .AND. M-P .GT. 0 ) THEN
@@ -559,7 +599,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
559599
END DO
560600
CALL CLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
561601
$ LDV1T )
562-
CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
602+
CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T,
603+
$ WORK(ITAUQ1),
563604
$ WORK(IORGLQ), LORGLQ, CHILDINFO )
564605
END IF
565606
*
@@ -602,7 +643,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
602643
U1(1,J) = ZERO
603644
U1(J,1) = ZERO
604645
END DO
605-
CALL CLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
646+
CALL CLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2),
647+
$ LDU1 )
606648
CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
607649
$ WORK(IORGQR), LORGQR, CHILDINFO )
608650
END IF
@@ -652,7 +694,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
652694
*
653695
IF( WANTU1 .AND. P .GT. 0 ) THEN
654696
CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
655-
CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
697+
CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1),
698+
$ WORK(IORGQR),
656699
$ LORGQR, CHILDINFO )
657700
END IF
658701
IF( WANTU2 .AND. M-P .GT. 0 ) THEN
@@ -736,7 +779,8 @@ SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
736779
END IF
737780
IF( WANTV1T .AND. Q .GT. 0 ) THEN
738781
CALL CLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
739-
CALL CLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
782+
CALL CLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1),
783+
$ LDX11,
740784
$ V1T(M-Q+1,M-Q+1), LDV1T )
741785
CALL CLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
742786
$ V1T(P+1,P+1), LDV1T )

‎lapack-netlib/SRC/dorcsd2by1.f‎

Lines changed: 51 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download DORCSD2BY1 + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -224,12 +222,14 @@
224222
*> \author Univ. of Colorado Denver
225223
*> \author NAG Ltd.
226224
*
227-
*> \ingroup doubleOTHERcomputational
225+
*> \ingroup uncsd2by1
228226
*
229227
* =====================================================================
230-
SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
228+
SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11,
229+
$ LDX11,
231230
$ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
232231
$ LDV1T, WORK, LWORK, IWORK, INFO )
232+
IMPLICIT NONE
233233
*
234234
* -- LAPACK computational routine (3.5.0) --
235235
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -266,7 +266,9 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
266266
DOUBLE PRECISION DUM1(1), DUM2(1,1)
267267
* ..
268268
* .. External Subroutines ..
269-
EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1,
269+
EXTERNAL DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT,
270+
$ DLASET,
271+
$ DORBDB1,
270272
$ DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR,
271273
$ XERBLA
272274
* ..
@@ -370,7 +372,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
370372
LORGLQMIN = MAX( LORGLQMIN, Q-1 )
371373
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
372374
END IF
373-
CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
375+
CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q,
376+
$ THETA,
374377
$ DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T,
375378
$ DUM2, 1, DUM1, DUM1, DUM1,
376379
$ DUM1, DUM1, DUM1, DUM1,
@@ -399,7 +402,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
399402
LORGLQMIN = MAX( LORGLQMIN, Q )
400403
LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
401404
END IF
402-
CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
405+
CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P,
406+
$ THETA,
403407
$ DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1,
404408
$ U2, LDU2, DUM1, DUM1, DUM1,
405409
$ DUM1, DUM1, DUM1, DUM1,
@@ -485,6 +489,36 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
485489
END IF
486490
LORGQR = LWORK-IORGQR+1
487491
LORGLQ = LWORK-IORGLQ+1
492+
*
493+
IF( R .EQ. 0 ) THEN
494+
*
495+
* R = 0: C and S are empty. Handle the trivial CSD directly.
496+
*
497+
IF( Q .EQ. 0 ) THEN
498+
IF( WANTU1 .AND. P .GT. 0 )
499+
$ CALL DLASET( 'A', P, P, ZERO, ONE, U1, LDU1 )
500+
IF( WANTU2 .AND. M-P .GT. 0 )
501+
$ CALL DLASET( 'A', M-P, M-P, ZERO, ONE, U2, LDU2 )
502+
RETURN
503+
END IF
504+
*
505+
IF( P .EQ. 0 .AND. M .EQ. Q ) THEN
506+
IF( WANTU2 )
507+
$ CALL DLACPY( 'A', M-P, Q, X21, LDX21, U2, LDU2 )
508+
IF( WANTV1T )
509+
$ CALL DLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T )
510+
RETURN
511+
END IF
512+
*
513+
IF( P .EQ. M .AND. M .EQ. Q ) THEN
514+
IF( WANTU1 )
515+
$ CALL DLACPY( 'A', P, Q, X11, LDX11, U1, LDU1 )
516+
IF( WANTV1T )
517+
$ CALL DLASET( 'A', Q, Q, ZERO, ONE, V1T, LDV1T )
518+
RETURN
519+
END IF
520+
*
521+
END IF
488522
*
489523
* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
490524
* in which R = MIN(P,M-P,Q,M-Q)
@@ -503,7 +537,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
503537
*
504538
IF( WANTU1 .AND. P .GT. 0 ) THEN
505539
CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
506-
CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
540+
CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1),
541+
$ WORK(IORGQR),
507542
$ LORGQR, CHILDINFO )
508543
END IF
509544
IF( WANTU2 .AND. M-P .GT. 0 ) THEN
@@ -519,7 +554,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
519554
END DO
520555
CALL DLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
521556
$ LDV1T )
522-
CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
557+
CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T,
558+
$ WORK(ITAUQ1),
523559
$ WORK(IORGLQ), LORGLQ, CHILDINFO )
524560
END IF
525561
*
@@ -562,7 +598,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
562598
U1(1,J) = ZERO
563599
U1(J,1) = ZERO
564600
END DO
565-
CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
601+
CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2),
602+
$ LDU1 )
566603
CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
567604
$ WORK(IORGQR), LORGQR, CHILDINFO )
568605
END IF
@@ -612,7 +649,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
612649
*
613650
IF( WANTU1 .AND. P .GT. 0 ) THEN
614651
CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
615-
CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
652+
CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1),
653+
$ WORK(IORGQR),
616654
$ LORGQR, CHILDINFO )
617655
END IF
618656
IF( WANTU2 .AND. M-P .GT. 0 ) THEN
@@ -695,7 +733,8 @@ SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
695733
END IF
696734
IF( WANTV1T .AND. Q .GT. 0 ) THEN
697735
CALL DLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
698-
CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
736+
CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1),
737+
$ LDX11,
699738
$ V1T(M-Q+1,M-Q+1), LDV1T )
700739
CALL DLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
701740
$ V1T(P+1,P+1), LDV1T )

0 commit comments

Comments
 (0)