Skip to content

Commit 10c2ebd

Browse files
wernsaarwernsaar
committed
BUGFIX: removed fixes for bugs #148 and #149, because info for xerbla is wrong
1 parent 26b3b3a commit 10c2ebd

4 files changed

Lines changed: 124 additions & 241 deletions

File tree

‎lapack-netlib/SRC/cgesvdx.f‎

Lines changed: 33 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@
170170
*> vectors, stored columnwise) as specified by RANGE; if
171171
*> JOBU = 'N', U is not referenced.
172172
*> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
173-
*> the exact value of NS is not known in advance and an upper
173+
*> the exact value of NS is not known ILQFin advance and an upper
174174
*> bound must be used.
175175
*> \endverbatim
176176
*>
@@ -294,8 +294,8 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
294294
CHARACTER JOBZ, RNGTGK
295295
LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
296296
INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
297-
$ ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ,
298-
$ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR
297+
$ ITAU, ITAUP, ITAUQ, ITEMP, ITGKZ, IUTGK,
298+
$ J, K, MAXWRK, MINMN, MINWRK, MNTHR
299299
REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
300300
* ..
301301
* .. Local Arrays ..
@@ -367,14 +367,8 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
367367
IF( INFO.EQ.0 ) THEN
368368
IF( WANTU .AND. LDU.LT.M ) THEN
369369
INFO = -15
370-
ELSE IF( WANTVT ) THEN
371-
IF( INDS ) THEN
372-
IF( LDVT.LT.IU-IL+1 ) THEN
373-
INFO = -17
374-
END IF
375-
ELSE IF( LDVT.LT.MINMN ) THEN
376-
INFO = -17
377-
END IF
370+
ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
371+
INFO = -16
378372
END IF
379373
END IF
380374
END IF
@@ -396,50 +390,37 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
396390
*
397391
* Path 1 (M much larger than N)
398392
*
399-
MINWRK = N*(N+5)
400-
MAXWRK = N + N*ILAENV(1,'CGEQRF',' ',M,N,-1,-1)
401-
MAXWRK = MAX(MAXWRK,
402-
$ N*N+2*N+2*N*ILAENV(1,'CGEBRD',' ',N,N,-1,-1))
403-
IF (WANTU .OR. WANTVT) THEN
404-
MAXWRK = MAX(MAXWRK,
405-
$ N*N+2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1))
406-
END IF
393+
MAXWRK = N + N*
394+
$ ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
395+
MAXWRK = MAX( MAXWRK, N*N + N + 2*N*
396+
$ ILAENV( 1, 'SGEBRD', ' ', N, N, -1, -1 ) )
397+
MINWRK = N*(N+4)
407398
ELSE
408399
*
409400
* Path 2 (M at least N, but not much larger)
410401
*
411-
MINWRK = 3*N + M
412-
MAXWRK = 2*N + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1)
413-
IF (WANTU .OR. WANTVT) THEN
414-
MAXWRK = MAX(MAXWRK,
415-
$ 2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1))
416-
END IF
402+
MAXWRK = 2*N + ( M+N )*
403+
$ ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 )
404+
MINWRK = 2*N + M
417405
END IF
418406
ELSE
419407
MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
420408
IF( N.GE.MNTHR ) THEN
421409
*
422410
* Path 1t (N much larger than M)
423411
*
424-
MINWRK = M*(M+5)
425-
MAXWRK = M + M*ILAENV(1,'CGELQF',' ',M,N,-1,-1)
426-
MAXWRK = MAX(MAXWRK,
427-
$ M*M+2*M+2*M*ILAENV(1,'CGEBRD',' ',M,M,-1,-1))
428-
IF (WANTU .OR. WANTVT) THEN
429-
MAXWRK = MAX(MAXWRK,
430-
$ M*M+2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1))
431-
END IF
412+
MAXWRK = M + M*
413+
$ ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
414+
MAXWRK = MAX( MAXWRK, M*M + M + 2*M*
415+
$ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
416+
MINWRK = M*(M+4)
432417
ELSE
433418
*
434419
* Path 2t (N greater than M, but not much larger)
435420
*
436-
*
437-
MINWRK = 3*M + N
438-
MAXWRK = 2*M + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1)
439-
IF (WANTU .OR. WANTVT) THEN
440-
MAXWRK = MAX(MAXWRK,
441-
$ 2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1))
442-
END IF
421+
MAXWRK = M*(M*2+19) + ( M+N )*
422+
$ ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 )
423+
MINWRK = 2*M + N
443424
END IF
444425
END IF
445426
END IF
@@ -537,14 +518,14 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
537518
CALL CGEBRD( N, N, WORK( IQRF ), N, RWORK( ID ),
538519
$ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
539520
$ WORK( ITEMP ), LWORK-ITEMP+1, INFO )
540-
ITEMPR = ITGKZ + N*(N*2+1)
521+
ITEMP = ITGKZ + N*(N*2+1)
541522
*
542523
* Solve eigenvalue problem TGK*Z=Z*S.
543524
* (Workspace: need 2*N*N+14*N)
544525
*
545526
CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
546527
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
547-
$ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
528+
$ RWORK( ITGKZ ), N*2, RWORK( ITEMP ),
548529
$ IWORK, INFO)
549530
*
550531
* If needed, compute left singular vectors.
@@ -558,7 +539,7 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
558539
END DO
559540
K = K + N
560541
END DO
561-
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
542+
CALL CLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU )
562543
*
563544
* Call CUNMBR to compute QB*UB.
564545
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@@ -613,14 +594,14 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
613594
CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),
614595
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
615596
$ LWORK-ITEMP+1, INFO )
616-
ITEMPR = ITGKZ + N*(N*2+1)
597+
ITEMP = ITGKZ + N*(N*2+1)
617598
*
618599
* Solve eigenvalue problem TGK*Z=Z*S.
619600
* (Workspace: need 2*N*N+14*N)
620601
*
621602
CALL SBDSVDX( 'U', JOBZ, RNGTGK, N, RWORK( ID ),
622603
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
623-
$ RWORK( ITGKZ ), N*2, RWORK( ITEMPR ),
604+
$ RWORK( ITGKZ ), N*2, RWORK( ITEMP ),
624605
$ IWORK, INFO)
625606
*
626607
* If needed, compute left singular vectors.
@@ -634,7 +615,7 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
634615
END DO
635616
K = K + N
636617
END DO
637-
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
618+
CALL CLASET( 'A', M-N, N, CZERO, CZERO, U( N+1,1 ), LDU )
638619
*
639620
* Call CUNMBR to compute QB*UB.
640621
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@@ -700,14 +681,14 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
700681
CALL CGEBRD( M, M, WORK( ILQF ), M, RWORK( ID ),
701682
$ RWORK( IE ), WORK( ITAUQ ), WORK( ITAUP ),
702683
$ WORK( ITEMP ), LWORK-ITEMP+1, INFO )
703-
ITEMPR = ITGKZ + M*(M*2+1)
684+
ITEMP = ITGKZ + M*(M*2+1)
704685
*
705686
* Solve eigenvalue problem TGK*Z=Z*S.
706687
* (Workspace: need 2*M*M+14*M)
707688
*
708689
CALL SBDSVDX( 'U', JOBZ, RNGTGK, M, RWORK( ID ),
709690
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
710-
$ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
691+
$ RWORK( ITGKZ ), M*2, RWORK( ITEMP ),
711692
$ IWORK, INFO)
712693
*
713694
* If needed, compute left singular vectors.
@@ -741,7 +722,7 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
741722
END DO
742723
K = K + M
743724
END DO
744-
CALL CLASET( 'A', NS, N-M, CZERO, CZERO,
725+
CALL CLASET( 'A', M, N-M, CZERO, CZERO,
745726
$ VT( 1,M+1 ), LDVT )
746727
*
747728
* Call CUNMBR to compute (VB**T)*(PB**T)
@@ -777,14 +758,14 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
777758
CALL CGEBRD( M, N, A, LDA, RWORK( ID ), RWORK( IE ),
778759
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( ITEMP ),
779760
$ LWORK-ITEMP+1, INFO )
780-
ITEMPR = ITGKZ + M*(M*2+1)
761+
ITEMP = ITGKZ + M*(M*2+1)
781762
*
782763
* Solve eigenvalue problem TGK*Z=Z*S.
783764
* (Workspace: need 2*M*M+14*M)
784765
*
785766
CALL SBDSVDX( 'L', JOBZ, RNGTGK, M, RWORK( ID ),
786767
$ RWORK( IE ), VL, VU, ILTGK, IUTGK, NS, S,
787-
$ RWORK( ITGKZ ), M*2, RWORK( ITEMPR ),
768+
$ RWORK( ITGKZ ), M*2, RWORK( ITEMP ),
788769
$ IWORK, INFO)
789770
*
790771
* If needed, compute left singular vectors.
@@ -818,7 +799,7 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
818799
END DO
819800
K = K + M
820801
END DO
821-
CALL CLASET( 'A', NS, N-M, CZERO, CZERO,
802+
CALL CLASET( 'A', M, N-M, CZERO, CZERO,
822803
$ VT( 1,M+1 ), LDVT )
823804
*
824805
* Call CUNMBR to compute VB**T * PB**T

‎lapack-netlib/SRC/dgesvdx.f‎

Lines changed: 18 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -169,7 +169,7 @@
169169
*> vectors, stored columnwise) as specified by RANGE; if
170170
*> JOBU = 'N', U is not referenced.
171171
*> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
172-
*> the exact value of NS is not known in advance and an upper
172+
*> the exact value of NS is not known ILQFin advance and an upper
173173
*> bound must be used.
174174
*> \endverbatim
175175
*>
@@ -357,14 +357,8 @@ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
357357
IF( INFO.EQ.0 ) THEN
358358
IF( WANTU .AND. LDU.LT.M ) THEN
359359
INFO = -15
360-
ELSE IF( WANTVT ) THEN
361-
IF( INDS ) THEN
362-
IF( LDVT.LT.IU-IL+1 ) THEN
363-
INFO = -17
364-
END IF
365-
ELSE IF( LDVT.LT.MINMN ) THEN
366-
INFO = -17
367-
END IF
360+
ELSE IF( WANTVT .AND. LDVT.LT.MINMN ) THEN
361+
INFO = -16
368362
END IF
369363
END IF
370364
END IF
@@ -386,69 +380,37 @@ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
386380
*
387381
* Path 1 (M much larger than N)
388382
*
389-
MAXWRK = N +
383+
MAXWRK = N*(N*2+16) +
390384
$ N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
391-
MAXWRK = MAX( MAXWRK, N*(N+5) + 2*N*
385+
MAXWRK = MAX( MAXWRK, N*(N*2+20) + 2*N*
392386
$ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
393-
IF (WANTU) THEN
394-
MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
395-
$ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
396-
END IF
397-
IF (WANTVT) THEN
398-
MAXWRK = MAX(MAXWRK,N*(N*3+6)+N*
399-
$ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
400-
END IF
401-
MINWRK = N*(N*3+20)
387+
MINWRK = N*(N*2+21)
402388
ELSE
403389
*
404390
* Path 2 (M at least N, but not much larger)
405391
*
406-
MAXWRK = 4*N + ( M+N )*
392+
MAXWRK = N*(N*2+19) + ( M+N )*
407393
$ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
408-
IF (WANTU) THEN
409-
MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
410-
$ ILAENV( 1, 'DORMQR', ' ', N, N, -1, -1 ) )
411-
END IF
412-
IF (WANTVT) THEN
413-
MAXWRK = MAX(MAXWRK,N*(N*2+5)+N*
414-
$ ILAENV( 1, 'DORMLQ', ' ', N, N, -1, -1 ) )
415-
END IF
416-
MINWRK = MAX(N*(N*2+19),4*N+M)
394+
MINWRK = N*(N*2+20) + M
417395
END IF
418396
ELSE
419397
MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
420398
IF( N.GE.MNTHR ) THEN
421399
*
422400
* Path 1t (N much larger than M)
423401
*
424-
MAXWRK = M +
402+
MAXWRK = M*(M*2+16) +
425403
$ M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
426-
MAXWRK = MAX( MAXWRK, M*(M+5) + 2*M*
404+
MAXWRK = MAX( MAXWRK, M*(M*2+20) + 2*M*
427405
$ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
428-
IF (WANTU) THEN
429-
MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
430-
$ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
431-
END IF
432-
IF (WANTVT) THEN
433-
MAXWRK = MAX(MAXWRK,M*(M*3+6)+M*
434-
$ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
435-
END IF
436-
MINWRK = M*(M*3+20)
406+
MINWRK = M*(M*2+21)
437407
ELSE
438408
*
439-
* Path 2t (N at least M, but not much larger)
409+
* Path 2t (N greater than M, but not much larger)
440410
*
441-
MAXWRK = 4*M + ( M+N )*
411+
MAXWRK = M*(M*2+19) + ( M+N )*
442412
$ ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 )
443-
IF (WANTU) THEN
444-
MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
445-
$ ILAENV( 1, 'DORMQR', ' ', M, M, -1, -1 ) )
446-
END IF
447-
IF (WANTVT) THEN
448-
MAXWRK = MAX(MAXWRK,M*(M*2+5)+M*
449-
$ ILAENV( 1, 'DORMLQ', ' ', M, M, -1, -1 ) )
450-
END IF
451-
MINWRK = MAX(M*(M*2+19),4*M+N)
413+
MINWRK = M*(M*2+20) + N
452414
END IF
453415
END IF
454416
END IF
@@ -560,7 +522,7 @@ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
560522
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
561523
J = J + N*2
562524
END DO
563-
CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
525+
CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
564526
*
565527
* Call DORMBR to compute QB*UB.
566528
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@@ -629,7 +591,7 @@ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
629591
CALL DCOPY( N, WORK( J ), 1, U( 1,I ), 1 )
630592
J = J + N*2
631593
END DO
632-
CALL DLASET( 'A', M-N, NS, ZERO, ZERO, U( N+1,1 ), LDU )
594+
CALL DLASET( 'A', M-N, N, ZERO, ZERO, U( N+1,1 ), LDU )
633595
*
634596
* Call DORMBR to compute QB*UB.
635597
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
@@ -725,7 +687,7 @@ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
725687
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
726688
J = J + M*2
727689
END DO
728-
CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
690+
CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
729691
*
730692
* Call DORMBR to compute (VB**T)*(PB**T)
731693
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
@@ -794,7 +756,7 @@ SUBROUTINE DGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
794756
CALL DCOPY( M, WORK( J ), 1, VT( I,1 ), LDVT )
795757
J = J + M*2
796758
END DO
797-
CALL DLASET( 'A', NS, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT)
759+
CALL DLASET( 'A', M, N-M, ZERO, ZERO, VT( 1,M+1 ), LDVT )
798760
*
799761
* Call DORMBR to compute VB**T * PB**T
800762
* (Workspace in WORK( ITEMP ): need M, prefer M*NB)

0 commit comments

Comments
 (0)