handle and document corner cases of lwork in lapack, align all precisions

Author: dklyuchinskiy

SRC/cgebrd.f

@@ -252,7 +252,7 @@ SUBROUTINE CGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
       ELSE
          LWKMIN = MAX( M, N )
          NB = MAX( 1, ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 ) )
-         LWKOPT = MAX( 1, ( M+N )*NB )
+         LWKOPT = ( M+N )*NB
       END IF
       WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
       LQUERY = ( LWORK.EQ.-1 )
@@ -292,7 +292,7 @@ SUBROUTINE CGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
 *        Determine when to switch from blocked to unblocked code.
 *
          IF( NX.LT.MINMN ) THEN
-            WS = ( M+N )*NB
+            WS = LWKOPT
             IF( LWORK.LT.WS ) THEN
 *
 *              Not enough work space for the optimal NB, consider using

SRC/cgehrd.f

@@ -222,12 +222,12 @@ SUBROUTINE CGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
          INFO = -8
       END IF
 *
+      NH = IHI - ILO + 1
       IF( INFO.EQ.0 ) THEN
 *
 *        Compute the workspace requirements
 *
-
-         IF( N.EQ.0 ) THEN
+         IF( NH.LE.1 ) THEN
             LWKOPT = 1
          ELSE
             NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI,
@@ -255,7 +255,6 @@ SUBROUTINE CGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
 *
 *     Quick return if possible
 *
-      NH = IHI - ILO + 1
       IF( NH.LE.1 ) THEN
          WORK( 1 ) = 1
          RETURN

SRC/cgelqf.f

@@ -185,7 +185,7 @@ SUBROUTINE CGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
          INFO = -2
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
-      ELSE IF ( .NOT.LQUERY ) THEN
+      ELSE IF( .NOT.LQUERY ) THEN
          IF( LWORK.LE.0 .OR. ( N.GT.0 .AND. LWORK.LT.MAX( 1, M ) ) )
      $      INFO = -7
       END IF

SRC/cgemlq.f

@@ -110,8 +110,8 @@
 *>
 *> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK))
-*>         On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
+*>          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
 *> \endverbatim
 *>
 *> \param[in] LWORK
@@ -227,7 +227,7 @@ SUBROUTINE CGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
       ELSE
          LWMIN = MAX( 1, LW )
       END IF
-
+*
       IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
         IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
           NBLCKS = ( MN - K ) / ( NB - K )

SRC/cgemqr.f

@@ -111,8 +111,8 @@
 *>
 *> \param[out] WORK
 *> \verbatim
-*>         (workspace) COMPLEX array, dimension (MAX(1,LWORK))
-*>         On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
+*>          (workspace) COMPLEX array, dimension (MAX(1,LWORK))
+*>          On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
 *> \endverbatim
 *>
 *> \param[in] LWORK

SRC/cgeqlf.f

@@ -192,7 +192,7 @@ SUBROUTINE CGEQLF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
 *
          IF( .NOT.LQUERY ) THEN
             IF( LWORK.LE.0 .OR. ( M.GT.0 .AND. LWORK.LT.MAX( 1, N ) ) )
-     $          INFO = -7
+     $         INFO = -7
          END IF
       END IF
 *

SRC/cgeqp3rk.f

@@ -428,7 +428,7 @@
 *> \verbatim
 *>          LWORK is INTEGER
 *>          The dimension of the array WORK.
-*>          LWORK >= 1, if MIN(M,N) = 0,
+*>          LWORK >= 1, if MIN(M,N) = 0, and
 *>          LWORK >= N+NRHS-1, otherwise.
 *>          For optimal performance LWORK >= NB*( N+NRHS+1 ),
 *>          where NB is the optimal block size for CGEQP3RK returned
@@ -628,8 +628,9 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
 *     .. External Functions ..
       LOGICAL            SISNAN
       INTEGER            ISAMAX, ILAENV
-      REAL               SLAMCH, SCNRM2
-      EXTERNAL           SISNAN, SLAMCH, SCNRM2, ISAMAX, ILAENV
+      REAL               SLAMCH, SCNRM2, SROUNDUP_LWORK
+      EXTERNAL           SISNAN, SLAMCH, SCNRM2, ISAMAX, ILAENV,
+     $                   SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          CMPLX, MAX, MIN
@@ -704,7 +705,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
 *
             LWKOPT = 2*N + NB*( N+NRHS+1 )
          END IF
-         WORK( 1 ) = CMPLX( LWKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
 *
          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
             INFO = -15
@@ -727,7 +728,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
          K = 0
          MAXC2NRMK = ZERO
          RELMAXC2NRMK = ZERO
-         WORK( 1 ) = CMPLX( LWKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
          RETURN
       END IF
 *
@@ -779,7 +780,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
 *
 *        Array TAU is not set and contains undefined elements.
 *
-         WORK( 1 ) = CMPLX( LWKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
          RETURN
       END IF
 *
@@ -798,7 +799,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
             TAU( J ) = CZERO
          END DO
 *
-         WORK( 1 ) = CMPLX( LWKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
          RETURN
 *
       END IF
@@ -829,7 +830,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
          DO J = 1, MINMN
             TAU( J ) = CZERO
          END DO
-         WORK( 1 ) = CMPLX( LWKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
          RETURN
       END IF
 *
@@ -874,7 +875,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
             TAU( J ) = CZERO
          END DO
 *
-         WORK( 1 ) = CMPLX( LWKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
          RETURN
       END IF
 *
@@ -992,7 +993,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
 *
 *              Return from the routine.
 *
-               WORK( 1 ) = CMPLX( LWKOPT )
+               WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
 *
                RETURN
 *
@@ -1083,7 +1084,7 @@ SUBROUTINE CGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA,
 *
       END IF
 *
-      WORK( 1 ) = CMPLX( LWKOPT )
+      WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
 *
       RETURN
 *

SRC/cgeqrfp.f

@@ -184,12 +184,12 @@ SUBROUTINE CGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO )
       INFO = 0
       NB = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
       K = MIN( M, N )
-      IF ( K.EQ.0 ) THEN
-          LWKMIN = 1
-          LWKOPT = 1
+      IF( K.EQ.0 ) THEN
+         LWKMIN = 1
+         LWKOPT = 1
       ELSE
-          LWKMIN = N
-          LWKOPT = N*NB
+         LWKMIN = N
+         LWKOPT = N*NB
       END IF
       WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
 *

SRC/cgesvj.f

@@ -208,16 +208,17 @@
 *> \verbatim
 *>          CWORK is COMPLEX array, dimension (max(1,LWORK))
 *>          Used as workspace.
-*>          If on entry LWORK = -1, then a workspace query is assumed and
-*>          no computation is done; CWORK(1) is set to the minial (and optimal)
-*>          length of CWORK.
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
 *>          LWORK is INTEGER.
 *>          Length of CWORK.
-*>          LWORK >= 1, if MIN(M,N) = 0, and LWORK >= MAX(1,M+N), otherwise.
+*>          LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M+N, otherwise.
+*>
+*>          If on entry LWORK = -1, then a workspace query is assumed and
+*>          no computation is done; CWORK(1) is set to the minial (and optimal)
+*>          length of CWORK.
 *> \endverbatim
 *>
 *> \param[in,out] RWORK
@@ -248,15 +249,17 @@
 *>          RWORK(6) = the largest absolute value over all sines of the
 *>                    Jacobi rotation angles in the last sweep. It can be
 *>                    useful for a post festum analysis.
-*>         If on entry LRWORK = -1, then a workspace query is assumed and
-*>         no computation is done; RWORK(1) is set to the minial (and optimal)
-*>         length of RWORK.
 *> \endverbatim
 *>
 *> \param[in] LRWORK
 *> \verbatim
 *>         LRWORK is INTEGER
-*>         Length of RWORK, LRWORK >= MAX(6,N).
+*>         Length of RWORK.
+*>         LRWORK >= 1, if MIN(M,N) = 0, and LRWORK >= MAX(6,N), otherwise
+*>
+*>         If on entry LRWORK = -1, then a workspace query is assumed and
+*>         no computation is done; RWORK(1) is set to the minial (and optimal)
+*>         length of RWORK.
 *> \endverbatim
 *>
 *> \param[out] INFO
@@ -400,8 +403,8 @@ SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
       INTEGER            ISAMAX
       EXTERNAL           ISAMAX
 *     from LAPACK
-      REAL               SLAMCH
-      EXTERNAL           SLAMCH
+      REAL               SLAMCH, SROUNDUP_LWORK
+      EXTERNAL           SLAMCH, SROUNDUP_LWORK
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     ..
@@ -423,19 +426,17 @@ SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
       APPLV = LSAME( JOBV, 'A' )
       UPPER = LSAME( JOBA, 'U' )
       LOWER = LSAME( JOBA, 'L' )
-
+*
       MINMN = MIN( M, N )
       IF( MINMN.EQ.0 ) THEN
-         LWMIN = 1
-         LRWMIN = 6
+         LWMIN  = 1
+         LRWMIN = 1
       ELSE
-         LWMIN = M + N
+         LWMIN  = M + N
          LRWMIN = MAX( 6, N )
       END IF
-      CWORK(1) = LWMIN
-      RWORK(1) = LRWMIN
 *
-      LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
+      LQUERY = ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 )
       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
          INFO = -1
       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
@@ -467,7 +468,9 @@ SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
       IF( INFO.NE.0 ) THEN
          CALL XERBLA( 'CGESVJ', -INFO )
          RETURN
-      ELSE IF ( LQUERY ) THEN
+      ELSE IF( LQUERY ) THEN
+         CWORK( 1 ) = SROUNDUP_LWORK( LWMIN )
+         RWORK( 1 ) = SROUNDUP_LWORK( LRWMIN )
          RETURN
       END IF
 *

SRC/cgetsqrhrt.f

@@ -131,13 +131,15 @@
 *> \param[in] LWORK
 *> \verbatim
 *>          The dimension of the array WORK.
+*>          If MIN(M,N) = 0, LWORK >= 1, else
 *>          LWORK >= MAX( 1, LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) ),
 *>          where
 *>             NUM_ALL_ROW_BLOCKS = CEIL((M-N)/(MB1-N)),
 *>             NB1LOCAL = MIN(NB1,N).
 *>             LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL,
 *>             LW1 = NB1LOCAL * N,
 *>             LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) ).
+*>
 *>          If LWORK = -1, then a workspace query is assumed.
 *>          The routine only calculates the optimal size of the WORK
 *>          array, returns this value as the first entry of the WORK
@@ -200,6 +202,10 @@ SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
       INTEGER            I, IINFO, J, LW1, LW2, LWT, LDWT, LWORKOPT,
      $                   NB1LOCAL, NB2LOCAL, NUM_ALL_ROW_BLOCKS
 *     ..
+*     .. External Functions ..
+      REAL               SROUNDUP_LWORK
+      EXTERNAL           SROUNDUP_LWORK
+*     ..
 *     .. External Subroutines ..
       EXTERNAL           CCOPY, CLATSQR, CUNGTSQR_ROW, CUNHR_COL,
      $                   XERBLA
@@ -212,7 +218,7 @@ SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
 *     Test the input arguments
 *
       INFO = 0
-      LQUERY  = ( LWORK.EQ.-1 )
+      LQUERY = ( LWORK.EQ.-1 )
       IF( M.LT.0 ) THEN
          INFO = -1
       ELSE IF( N.LT.0 .OR. M.LT.N ) THEN
@@ -225,7 +231,7 @@ SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
          INFO = -5
       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -7
-      ELSE IF( LDT.LT.MAX( 1,  MIN( NB2, N ) ) ) THEN
+      ELSE IF( LDT.LT.MAX( 1, MIN( NB2, N ) ) ) THEN
          INFO = -9
       ELSE
 *
@@ -278,14 +284,14 @@ SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
          CALL XERBLA( 'CGETSQRHRT', -INFO )
          RETURN
       ELSE IF ( LQUERY ) THEN
-         WORK( 1 ) = CMPLX( LWORKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWORKOPT )
          RETURN
       END IF
 *
 *     Quick return if possible
 *
       IF( MIN( M, N ).EQ.0 ) THEN
-         WORK( 1 ) = CMPLX( LWORKOPT )
+         WORK( 1 ) = SROUNDUP_LWORK( LWORKOPT )
          RETURN
       END IF
 *
@@ -342,7 +348,7 @@ SUBROUTINE CGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
          END IF
       END DO
 *
-      WORK( 1 ) = CMPLX( LWORKOPT )
+      WORK( 1 ) = SROUNDUP_LWORK( LWORKOPT )
       RETURN
 *
 *     End of CGETSQRHRT