Handle corner cases of LWORK (Reference-LAPACK PR 942)

This commit is contained in:
Martin Kroeker
2023-12-23 19:37:03 +01:00
committed by GitHub
parent 5c11b2ff41
commit 0814491d96
47 changed files with 783 additions and 533 deletions

View File

@@ -215,7 +215,8 @@
*> \param[in] LWORK
*> \verbatim
*> LWORK is INTEGER
*> The dimension of the array WORK.
*> The dimension of the array WORK. LWORK >= MAX(1,2*N).
*> For good performance, LWORK must generally be larger.
*>
*> If LWORK = -1, then a workspace query is assumed; the routine
*> only calculates the optimal size of the WORK array, returns
@@ -260,7 +261,7 @@
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complexGEeigen
*> \ingroup gges3
*
* =====================================================================
SUBROUTINE CGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
@@ -300,7 +301,8 @@
LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
$ LQUERY, WANTST
INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
$ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT
$ ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKOPT,
$ LWKMIN
REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
$ PVSR, SMLNUM
* ..
@@ -310,13 +312,12 @@
* ..
* .. External Subroutines ..
EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHD3, CLAQZ0, CLACPY,
$ CLASCL, CLASET, CTGSEN, CUNGQR, CUNMQR, SLABAD,
$ XERBLA
$ CLASCL, CLASET, CTGSEN, CUNGQR, CUNMQR, XERBLA
* ..
* .. External Functions ..
LOGICAL LSAME
REAL CLANGE, SLAMCH
EXTERNAL LSAME, CLANGE, SLAMCH
REAL CLANGE, SLAMCH, SROUNDUP_LWORK
EXTERNAL LSAME, CLANGE, SLAMCH, SROUNDUP_LWORK
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, SQRT
@@ -353,6 +354,8 @@
*
INFO = 0
LQUERY = ( LWORK.EQ.-1 )
LWKMIN = MAX( 1, 2*N )
*
IF( IJOBVL.LE.0 ) THEN
INFO = -1
ELSE IF( IJOBVR.LE.0 ) THEN
@@ -369,7 +372,7 @@
INFO = -14
ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
INFO = -16
ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
INFO = -18
END IF
*
@@ -377,29 +380,33 @@
*
IF( INFO.EQ.0 ) THEN
CALL CGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
LWKOPT = MAX( 1, N + INT ( WORK( 1 ) ) )
LWKOPT = MAX( LWKMIN, N + INT( WORK( 1 ) ) )
CALL CUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK,
$ -1, IERR )
LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
LWKOPT = MAX( LWKOPT, N + INT( WORK( 1 ) ) )
IF( ILVSL ) THEN
CALL CUNGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1,
$ IERR )
LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
LWKOPT = MAX( LWKOPT, N + INT( WORK( 1 ) ) )
END IF
CALL CGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL,
$ LDVSL, VSR, LDVSR, WORK, -1, IERR )
LWKOPT = MAX( LWKOPT, N + INT ( WORK( 1 ) ) )
LWKOPT = MAX( LWKOPT, N + INT( WORK( 1 ) ) )
CALL CLAQZ0( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB,
$ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, -1,
$ RWORK, 0, IERR )
LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) )
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
IF( WANTST ) THEN
CALL CTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
$ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM,
$ PVSL, PVSR, DIF, WORK, -1, IDUM, 1, IERR )
LWKOPT = MAX( LWKOPT, INT ( WORK( 1 ) ) )
LWKOPT = MAX( LWKOPT, INT( WORK( 1 ) ) )
END IF
IF( N.EQ.0 ) THEN
WORK( 1 ) = 1
ELSE
WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
END IF
WORK( 1 ) = CMPLX( LWKOPT )
END IF
*
@@ -422,7 +429,6 @@
EPS = SLAMCH( 'P' )
SMLNUM = SLAMCH( 'S' )
BIGNUM = ONE / SMLNUM
CALL SLABAD( SMLNUM, BIGNUM )
SMLNUM = SQRT( SMLNUM ) / EPS
BIGNUM = ONE / SMLNUM
*
@@ -585,7 +591,7 @@
*
30 CONTINUE
*
WORK( 1 ) = CMPLX( LWKOPT )
WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
*
RETURN
*