Skip to content

Commit 6032a6b

Browse files
dklyuchinskiydklyuchinskiy
committed
handle and document corner cases of lwork in lapack, double precision
1 parent a4298d5 commit 6032a6b

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

50 files changed

+443
-253
lines changed

‎SRC/dgebrd.f‎

Lines changed: 16 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,8 @@
122122
*> \param[in] LWORK
123123
*> \verbatim
124124
*> LWORK is INTEGER
125-
*> The length of the array WORK. LWORK >= max(1,M,N).
125+
*> The length of the array WORK.
126+
*> LWORK >= 1, if MIN(M,N) = 0, and LWORK >= MAX(M,N), otherwise.
126127
*> For optimum performance LWORK >= (M+N)*NB, where NB
127128
*> is the optimal blocksize.
128129
*>
@@ -223,8 +224,8 @@ SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
223224
* ..
224225
* .. Local Scalars ..
225226
LOGICAL LQUERY
226-
INTEGER I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
227-
$ NBMIN, NX, WS
227+
INTEGER I, IINFO, J, LDWRKX, LDWRKY, LWKMIN, LWKOPT,
228+
$ MINMN, NB, NBMIN, NX, WS
228229
* ..
229230
* .. External Subroutines ..
230231
EXTERNAL DGEBD2, DGEMM, DLABRD, XERBLA
@@ -241,17 +242,25 @@ SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
241242
* Test the input parameters
242243
*
243244
INFO = 0
244-
NB = MAX( 1, ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) )
245-
LWKOPT = ( M+N )*NB
245+
MINMN = MIN( M, N )
246+
IF( MINMN.EQ.0 ) THEN
247+
LWKMIN = 1
248+
LWKOPT = 1
249+
ELSE
250+
LWKMIN = MAX( M, N )
251+
NB = MAX( 1, ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) )
252+
LWKOPT = ( M+N )*NB
253+
ENDIF
246254
WORK( 1 ) = DBLE( LWKOPT )
255+
*
247256
LQUERY = ( LWORK.EQ.-1 )
248257
IF( M.LT.0 ) THEN
249258
INFO = -1
250259
ELSE IF( N.LT.0 ) THEN
251260
INFO = -2
252261
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
253262
INFO = -4
254-
ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
263+
ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
255264
INFO = -10
256265
END IF
257266
IF( INFO.LT.0 ) THEN
@@ -263,7 +272,6 @@ SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
263272
*
264273
* Quick return if possible
265274
*
266-
MINMN = MIN( M, N )
267275
IF( MINMN.EQ.0 ) THEN
268276
WORK( 1 ) = 1
269277
RETURN
@@ -282,7 +290,7 @@ SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
282290
* Determine when to switch from blocked to unblocked code.
283291
*
284292
IF( NX.LT.MINMN ) THEN
285-
WS = ( M+N )*NB
293+
WS = LWKOPT
286294
IF( LWORK.LT.WS ) THEN
287295
*
288296
* Not enough work space for the optimal NB, consider using

‎SRC/dgehrd.f‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -89,7 +89,7 @@
8989
*>
9090
*> \param[out] WORK
9191
*> \verbatim
92-
*> WORK is DOUBLE PRECISION array, dimension (LWORK)
92+
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
9393
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
9494
*> \endverbatim
9595
*>
@@ -225,8 +225,13 @@ SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
225225
*
226226
* Compute the workspace requirements
227227
*
228-
NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI, -1 ) )
229-
LWKOPT = N*NB + TSIZE
228+
IF( N.EQ.0 ) THEN
229+
LWKOPT = 1
230+
ELSE
231+
NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI,
232+
$ -1 ) )
233+
LWKOPT = N*NB + TSIZE
234+
ENDIF
230235
WORK( 1 ) = LWKOPT
231236
END IF
232237
*
@@ -344,6 +349,7 @@ SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
344349
* Use unblocked code to reduce the rest of the matrix
345350
*
346351
CALL DGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
352+
*
347353
WORK( 1 ) = LWKOPT
348354
*
349355
RETURN

‎SRC/dgelq.f‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -98,7 +98,7 @@
9898
*> \param[in] LWORK
9999
*> \verbatim
100100
*> LWORK is INTEGER
101-
*> The dimension of the array WORK.
101+
*> The dimension of the array WORK. LWORK >= 1.
102102
*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
103103
*> only calculates the sizes of the T and WORK arrays, returns these
104104
*> values as the first entries of the T and WORK arrays, and no error

‎SRC/dgelqf.f‎

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -93,7 +93,8 @@
9393
*> \param[in] LWORK
9494
*> \verbatim
9595
*> LWORK is INTEGER
96-
*> The dimension of the array WORK. LWORK >= max(1,M).
96+
*> The dimension of the array WORK.
97+
*> LWORK >= 1, if MIN(M,N) = 0, and LWORK >= M, otherwise.
9798
*> For optimum performance LWORK >= M*NB, where NB is the
9899
*> optimal blocksize.
99100
*>
@@ -174,29 +175,34 @@ SUBROUTINE DGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
174175
* Test the input arguments
175176
*
176177
INFO = 0
178+
K = MIN( M, N )
177179
NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
178-
LWKOPT = M*NB
179-
WORK( 1 ) = LWKOPT
180180
LQUERY = ( LWORK.EQ.-1 )
181181
IF( M.LT.0 ) THEN
182182
INFO = -1
183183
ELSE IF( N.LT.0 ) THEN
184184
INFO = -2
185185
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
186186
INFO = -4
187-
ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
188-
INFO = -7
187+
ELSE IF( .NOT.LQUERY ) THEN
188+
IF( LWORK.LE.0 .OR. ( N.GT.0 .AND. LWORK.LT.MAX( 1, M ) ) )
189+
$ INFO = -7
189190
END IF
190191
IF( INFO.NE.0 ) THEN
191192
CALL XERBLA( 'DGELQF', -INFO )
192193
RETURN
193194
ELSE IF( LQUERY ) THEN
195+
IF( K.EQ.0 ) THEN
196+
LWKOPT = 1
197+
ELSE
198+
LWKOPT = M*NB
199+
END IF
200+
WORK( 1 ) = LWKOPT
194201
RETURN
195202
END IF
196203
*
197204
* Quick return if possible
198205
*
199-
K = MIN( M, N )
200206
IF( K.EQ.0 ) THEN
201207
WORK( 1 ) = 1
202208
RETURN

‎SRC/dgelsd.f‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -276,7 +276,7 @@ SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
276276
$ LOG( TWO ) ) + 1, 0 )
277277
*
278278
IF( INFO.EQ.0 ) THEN
279-
MAXWRK = 0
279+
MAXWRK = 1
280280
LIWORK = 3*MINMN*NLVL + 11*MINMN
281281
MM = M
282282
IF( M.GE.N .AND. M.GE.MNTHR ) THEN

‎SRC/dgemlq.f‎

Lines changed: 16 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -111,13 +111,14 @@
111111
*>
112112
*> \param[out] WORK
113113
*> \verbatim
114-
*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
114+
*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
115+
*> On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
115116
*> \endverbatim
116117
*>
117118
*> \param[in] LWORK
118119
*> \verbatim
119120
*> LWORK is INTEGER
120-
*> The dimension of the array WORK.
121+
*> The dimension of the array WORK. LWORK >= 1.
121122
*> If LWORK = -1, then a workspace query is assumed. The routine
122123
*> only calculates the size of the WORK array, returns this
123124
*> value as WORK(1), and no error message related to WORK
@@ -188,7 +189,7 @@ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
188189
* ..
189190
* .. Local Scalars ..
190191
LOGICAL LEFT, RIGHT, TRAN, NOTRAN, LQUERY
191-
INTEGER MB, NB, LW, NBLCKS, MN
192+
INTEGER MB, NB, LW, NBLCKS, MN, MINMNK, LWMIN
192193
* ..
193194
* .. External Functions ..
194195
LOGICAL LSAME
@@ -204,7 +205,7 @@ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
204205
*
205206
* Test the input arguments
206207
*
207-
LQUERY = LWORK.EQ.-1
208+
LQUERY = ( LWORK.EQ.-1 )
208209
NOTRAN = LSAME( TRANS, 'N' )
209210
TRAN = LSAME( TRANS, 'T' )
210211
LEFT = LSAME( SIDE, 'L' )
@@ -219,6 +220,13 @@ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
219220
LW = M * MB
220221
MN = N
221222
END IF
223+
*
224+
MINMNK = MIN( M, N, K )
225+
IF( MINMNK.EQ.0 ) THEN
226+
LWMIN = 1
227+
ELSE
228+
LWMIN = MAX( 1, LW )
229+
END IF
222230
*
223231
IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
224232
IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
@@ -247,12 +255,12 @@ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
247255
INFO = -9
248256
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
249257
INFO = -11
250-
ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
258+
ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
251259
INFO = -13
252260
END IF
253261
*
254262
IF( INFO.EQ.0 ) THEN
255-
WORK( 1 ) = LW
263+
WORK( 1 ) = LWMIN
256264
END IF
257265
*
258266
IF( INFO.NE.0 ) THEN
@@ -264,7 +272,7 @@ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
264272
*
265273
* Quick return if possible
266274
*
267-
IF( MIN( M, N, K ).EQ.0 ) THEN
275+
IF( MINMNK.EQ.0 ) THEN
268276
RETURN
269277
END IF
270278
*
@@ -277,7 +285,7 @@ SUBROUTINE DGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
277285
$ MB, C, LDC, WORK, LWORK, INFO )
278286
END IF
279287
*
280-
WORK( 1 ) = LW
288+
WORK( 1 ) = LWMIN
281289
*
282290
RETURN
283291
*

‎SRC/dgemqr.f‎

Lines changed: 16 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -111,13 +111,14 @@
111111
*>
112112
*> \param[out] WORK
113113
*> \verbatim
114-
*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
114+
*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
115+
*> On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
115116
*> \endverbatim
116117
*>
117118
*> \param[in] LWORK
118119
*> \verbatim
119120
*> LWORK is INTEGER
120-
*> The dimension of the array WORK.
121+
*> The dimension of the array WORK. LWORK >= 1.
121122
*> If LWORK = -1, then a workspace query is assumed. The routine
122123
*> only calculates the size of the WORK array, returns this
123124
*> value as WORK(1), and no error message related to WORK
@@ -189,7 +190,7 @@ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
189190
* ..
190191
* .. Local Scalars ..
191192
LOGICAL LEFT, RIGHT, TRAN, NOTRAN, LQUERY
192-
INTEGER MB, NB, LW, NBLCKS, MN
193+
INTEGER MB, NB, LW, NBLCKS, MN, MINMNK, LWMIN
193194
* ..
194195
* .. External Functions ..
195196
LOGICAL LSAME
@@ -205,7 +206,7 @@ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
205206
*
206207
* Test the input arguments
207208
*
208-
LQUERY = LWORK.EQ.-1
209+
LQUERY = ( LWORK.EQ.-1 )
209210
NOTRAN = LSAME( TRANS, 'N' )
210211
TRAN = LSAME( TRANS, 'T' )
211212
LEFT = LSAME( SIDE, 'L' )
@@ -220,6 +221,13 @@ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
220221
LW = MB * NB
221222
MN = N
222223
END IF
224+
*
225+
MINMNK = MIN( M, N, K )
226+
IF( MINMNK.EQ.0 ) THEN
227+
LWMIN = 1
228+
ELSE
229+
LWMIN = MAX( 1, LW )
230+
END IF
223231
*
224232
IF( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
225233
IF( MOD( MN - K, MB - K ).EQ.0 ) THEN
@@ -248,12 +256,12 @@ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
248256
INFO = -9
249257
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
250258
INFO = -11
251-
ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
259+
ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
252260
INFO = -13
253261
END IF
254262
*
255263
IF( INFO.EQ.0 ) THEN
256-
WORK( 1 ) = LW
264+
WORK( 1 ) = LWMIN
257265
END IF
258266
*
259267
IF( INFO.NE.0 ) THEN
@@ -265,7 +273,7 @@ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
265273
*
266274
* Quick return if possible
267275
*
268-
IF( MIN( M, N, K ).EQ.0 ) THEN
276+
IF( MINMNK.EQ.0 ) THEN
269277
RETURN
270278
END IF
271279
*
@@ -278,7 +286,7 @@ SUBROUTINE DGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
278286
$ NB, C, LDC, WORK, LWORK, INFO )
279287
END IF
280288
*
281-
WORK( 1 ) = LW
289+
WORK( 1 ) = LWMIN
282290
*
283291
RETURN
284292
*

‎SRC/dgeqlf.f‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,8 @@
8888
*> \param[in] LWORK
8989
*> \verbatim
9090
*> LWORK is INTEGER
91-
*> The dimension of the array WORK. LWORK >= max(1,N).
91+
*> The dimension of the array WORK.
92+
*> LWORK >= 1, if MIN(M,N) = 0, and LWORK >= N, otherwise.
9293
*> For optimum performance LWORK >= N*NB, where NB is the
9394
*> optimal blocksize.
9495
*>
@@ -188,8 +189,9 @@ SUBROUTINE DGEQLF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
188189
END IF
189190
WORK( 1 ) = LWKOPT
190191
*
191-
IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
192-
INFO = -7
192+
IF( .NOT.LQUERY ) THEN
193+
IF( LWORK.LE.0 .OR. ( M.GT.0 .AND. LWORK.LT.MAX( 1, N ) ) )
194+
$ INFO = -7
193195
END IF
194196
END IF
195197
*

‎SRC/dgeqp3rk.f‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -427,7 +427,8 @@
427427
*> \verbatim
428428
*> LWORK is INTEGER
429429
*> The dimension of the array WORK.
430-
*. LWORK >= (3*N + NRHS - 1)
430+
*> LWORK >= 1, if MIN(M,N) = 0,
431+
*> LWORK >= (3*N + NRHS - 1), otherwise.
431432
*> For optimal performance LWORK >= (2*N + NB*( N+NRHS+1 )),
432433
*> where NB is the optimal block size for DGEQP3RK returned
433434
*> by ILAENV. Minimal block size MINNB=2.

‎SRC/dgeqr.f‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -99,7 +99,7 @@
9999
*> \param[in] LWORK
100100
*> \verbatim
101101
*> LWORK is INTEGER
102-
*> The dimension of the array WORK.
102+
*> The dimension of the array WORK. LWORK >= 1.
103103
*> If LWORK = -1 or -2, then a workspace query is assumed. The routine
104104
*> only calculates the sizes of the T and WORK arrays, returns these
105105
*> values as the first entries of the T and WORK arrays, and no error

0 commit comments

Comments
 (0)