Skip to content

Avoid intermediate overflows in GEJSV,GESVJ,GESVDX,GEDMD (Reference-LAPACK PR 1283) #5855

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

Open
martin-frbg wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

Avoid intermediate overflows in GEJSV,GESVJ,GESVDX,GEDMD (Reference-LAPACK PR 1283) #5855

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
2,305 changes: 1,154 additions & 1,151 deletions lapack-netlib/SRC/cgedmd.f90
View file
Open in desktop

Large diffs are not rendered by default.

147 changes: 94 additions & 53 deletions lapack-netlib/SRC/cgejsv.f
View file
Open in desktop

Large diffs are not rendered by default.

50 changes: 31 additions & 19 deletions lapack-netlib/SRC/cgesvdx.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 CGESVDX + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvdx.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvdx.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvdx.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
Expand Down Expand Up @@ -267,6 +265,7 @@
SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
$ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
$ LWORK, RWORK, IWORK, INFO )
IMPLICIT NONE
*
* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
Expand Down Expand Up @@ -305,15 +304,17 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
REAL DUM( 1 )
* ..
* .. External Subroutines ..
EXTERNAL CGEBRD, CGELQF, CGEQRF, CLASCL, CLASET,
EXTERNAL CGEBRD, CGELQF, CGEQRF, CLASCL,
$ CLASET,
$ CUNMBR, CUNMQR, CUNMLQ, CLACPY,
$ SBDSVDX, SLASCL, XERBLA
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
REAL SLAMCH, CLANGE, SROUNDUP_LWORK
EXTERNAL LSAME, ILAENV, SLAMCH, CLANGE, SROUNDUP_LWORK
EXTERNAL LSAME, ILAENV, SLAMCH,
$ CLANGE, SROUNDUP_LWORK
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN, SQRT
Expand Down Expand Up @@ -395,51 +396,59 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
MAXWRK = 1
IF( MINMN.GT.0 ) THEN
IF( M.GE.N ) THEN
MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0,
$ 0 )
IF( M.GE.MNTHR ) THEN
*
* Path 1 (M much larger than N)
*
MINWRK = N*(N+5)
MAXWRK = N + N*ILAENV(1,'CGEQRF',' ',M,N,-1,-1)
MAXWRK = MAX(MAXWRK,
$ N*N+2*N+2*N*ILAENV(1,'CGEBRD',' ',N,N,-1,-1))
$ N*N+2*N+2*N*ILAENV(1,'CGEBRD',' ',N,N,-1,
$ -1))
IF (WANTU .OR. WANTVT) THEN
MAXWRK = MAX(MAXWRK,
$ N*N+2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1))
$ N*N+2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,
$ -1))
END IF
ELSE
*
* Path 2 (M at least N, but not much larger)
*
MINWRK = 3*N + M
MAXWRK = 2*N + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1)
MAXWRK = 2*N + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,
$ -1)
IF (WANTU .OR. WANTVT) THEN
MAXWRK = MAX(MAXWRK,
$ 2*N+N*ILAENV(1,'CUNMQR','LN',N,N,N,-1))
END IF
END IF
ELSE
MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0, 0 )
MNTHR = ILAENV( 6, 'CGESVD', JOBU // JOBVT, M, N, 0,
$ 0 )
IF( N.GE.MNTHR ) THEN
*
* Path 1t (N much larger than M)
*
MINWRK = M*(M+5)
MAXWRK = M + M*ILAENV(1,'CGELQF',' ',M,N,-1,-1)
MAXWRK = MAX(MAXWRK,
$ M*M+2*M+2*M*ILAENV(1,'CGEBRD',' ',M,M,-1,-1))
$ M*M+2*M+2*M*ILAENV(1,'CGEBRD',' ',M,M,-1,
$ -1))
IF (WANTU .OR. WANTVT) THEN
MAXWRK = MAX(MAXWRK,
$ M*M+2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1))
$ M*M+2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,
$ -1))
END IF
ELSE
*
* Path 2t (N greater than M, but not much larger)
*
*
MINWRK = 3*M + N
MAXWRK = 2*M + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,-1)
MAXWRK = 2*M + (M+N)*ILAENV(1,'CGEBRD',' ',M,N,-1,
$ -1)
IF (WANTU .OR. WANTVT) THEN
MAXWRK = MAX(MAXWRK,
$ 2*M+M*ILAENV(1,'CUNMQR','LN',M,M,M,-1))
Expand Down Expand Up @@ -560,7 +569,8 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
END DO
K = K + N
END DO
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ),
$ LDU)
*
* Call CUNMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
Expand Down Expand Up @@ -636,7 +646,8 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
END DO
K = K + N
END DO
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ), LDU)
CALL CLASET( 'A', M-N, NS, CZERO, CZERO, U( N+1,1 ),
$ LDU)
*
* Call CUNMBR to compute QB*UB.
* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
Expand Down Expand Up @@ -835,13 +846,14 @@ SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
*
* Undo scaling if necessary
*
IF( ISCL.EQ.1 ) THEN
IF( ANRM.GT.BIGNUM )
$ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1,
IF( ISCL.EQ.1 .AND. NS.GT.0 ) THEN
IF( ANRM.GT.BIGNUM ) THEN
CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, NS, 1,
$ S, MINMN, INFO )
IF( ANRM.LT.SMLNUM )
$ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1,
ELSE IF( ANRM.LT.SMLNUM ) THEN
CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, NS, 1,
$ S, MINMN, INFO )
ENDIF
END IF
*
* Return optimal workspace in WORK(1)
Expand Down
Loading
Loading