Skip to content

Commit 0c1e0c3

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

39 files changed

+471
-275
lines changed

‎SRC/sgebrd.f‎

Lines changed: 17 additions & 9 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 SGEBRD( 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 SGEBD2, SGEMM, SLABRD, XERBLA
@@ -242,17 +243,24 @@ SUBROUTINE SGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
242243
* Test the input parameters
243244
*
244245
INFO = 0
245-
NB = MAX( 1, ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) )
246-
LWKOPT = ( M+N )*NB
247-
WORK( 1 ) = SROUNDUP_LWORK(LWKOPT)
246+
MINMN = MIN( M, N )
247+
IF( MINMN.EQ.0 ) THEN
248+
LWKMIN = 1
249+
LWKOPT = 1
250+
ELSE
251+
LWKMIN = MAX( M, N )
252+
NB = MAX( 1, ILAENV( 1, 'SGEBRD', ' ', M, N, -1, -1 ) )
253+
LWKOPT = ( M+N )*NB
254+
ENDIF
255+
WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
248256
LQUERY = ( LWORK.EQ.-1 )
249257
IF( M.LT.0 ) THEN
250258
INFO = -1
251259
ELSE IF( N.LT.0 ) THEN
252260
INFO = -2
253261
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
254262
INFO = -4
255-
ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
263+
ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
256264
INFO = -10
257265
END IF
258266
IF( INFO.LT.0 ) THEN
@@ -264,7 +272,6 @@ SUBROUTINE SGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
264272
*
265273
* Quick return if possible
266274
*
267-
MINMN = MIN( M, N )
268275
IF( MINMN.EQ.0 ) THEN
269276
WORK( 1 ) = 1
270277
RETURN
@@ -342,7 +349,8 @@ SUBROUTINE SGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
342349
*
343350
CALL SGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
344351
$ TAUQ( I ), TAUP( I ), WORK, IINFO )
345-
WORK( 1 ) = SROUNDUP_LWORK(WS)
352+
*
353+
WORK( 1 ) = SROUNDUP_LWORK( WS )
346354
RETURN
347355
*
348356
* End of SGEBRD

‎SRC/sgehrd.f‎

Lines changed: 14 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -89,7 +89,7 @@
8989
*>
9090
*> \param[out] WORK
9191
*> \verbatim
92-
*> WORK is REAL array, dimension (LWORK)
92+
*> WORK is REAL array, dimension (MAX(1,LWORK))
9393
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
9494
*> \endverbatim
9595
*>
@@ -173,7 +173,7 @@ SUBROUTINE SGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
173173
INTEGER IHI, ILO, INFO, LDA, LWORK, N
174174
* ..
175175
* .. Array Arguments ..
176-
REAL A( LDA, * ), TAU( * ), WORK( * )
176+
REAL A( LDA, * ), TAU( * ), WORK( * )
177177
* ..
178178
*
179179
* =====================================================================
@@ -182,15 +182,15 @@ SUBROUTINE SGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
182182
INTEGER NBMAX, LDT, TSIZE
183183
PARAMETER ( NBMAX = 64, LDT = NBMAX+1,
184184
$ TSIZE = LDT*NBMAX )
185-
REAL ZERO, ONE
185+
REAL ZERO, ONE
186186
PARAMETER ( ZERO = 0.0E+0,
187187
$ ONE = 1.0E+0 )
188188
* ..
189189
* .. Local Scalars ..
190190
LOGICAL LQUERY
191191
INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
192192
$ NBMIN, NH, NX
193-
REAL EI
193+
REAL EI
194194
* ..
195195
* .. External Subroutines ..
196196
EXTERNAL SAXPY, SGEHD2, SGEMM, SLAHR2, SLARFB, STRMM,
@@ -226,9 +226,14 @@ SUBROUTINE SGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
226226
*
227227
* Compute the workspace requirements
228228
*
229-
NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI, -1 ) )
230-
LWKOPT = N*NB + TSIZE
231-
WORK( 1 ) = SROUNDUP_LWORK(LWKOPT)
229+
IF( N.EQ.0 ) THEN
230+
LWKOPT = 1
231+
ELSE
232+
NB = MIN( NBMAX, ILAENV( 1, 'SGEHRD', ' ', N, ILO, IHI,
233+
$ -1 ) )
234+
LWKOPT = N*NB + TSIZE
235+
ENDIF
236+
WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
232237
END IF
233238
*
234239
IF( INFO.NE.0 ) THEN
@@ -345,7 +350,8 @@ SUBROUTINE SGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
345350
* Use unblocked code to reduce the rest of the matrix
346351
*
347352
CALL SGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
348-
WORK( 1 ) = SROUNDUP_LWORK(LWKOPT)
353+
*
354+
WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
349355
*
350356
RETURN
351357
*

‎SRC/sgelq.f‎

Lines changed: 4 additions & 4 deletions
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
@@ -295,9 +295,9 @@ SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
295295
T( 2 ) = MB
296296
T( 3 ) = NB
297297
IF( MINW ) THEN
298-
WORK( 1 ) = SROUNDUP_LWORK(LWMIN)
298+
WORK( 1 ) = SROUNDUP_LWORK( LWMIN )
299299
ELSE
300-
WORK( 1 ) = SROUNDUP_LWORK(LWREQ)
300+
WORK( 1 ) = SROUNDUP_LWORK( LWREQ )
301301
END IF
302302
END IF
303303
IF( INFO.NE.0 ) THEN
@@ -322,7 +322,7 @@ SUBROUTINE SGELQ( M, N, A, LDA, T, TSIZE, WORK, LWORK,
322322
$ LWORK, INFO )
323323
END IF
324324
*
325-
WORK( 1 ) = SROUNDUP_LWORK(LWREQ)
325+
WORK( 1 ) = SROUNDUP_LWORK( LWREQ )
326326
RETURN
327327
*
328328
* End of SGELQ

‎SRC/sgelqf.f‎

Lines changed: 13 additions & 7 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
*>
@@ -175,29 +176,34 @@ SUBROUTINE SGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
175176
* Test the input arguments
176177
*
177178
INFO = 0
179+
K = MIN( M, N )
178180
NB = ILAENV( 1, 'SGELQF', ' ', M, N, -1, -1 )
179-
LWKOPT = M*NB
180-
WORK( 1 ) = SROUNDUP_LWORK(LWKOPT)
181181
LQUERY = ( LWORK.EQ.-1 )
182182
IF( M.LT.0 ) THEN
183183
INFO = -1
184184
ELSE IF( N.LT.0 ) THEN
185185
INFO = -2
186186
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
187187
INFO = -4
188-
ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
189-
INFO = -7
188+
ELSE IF( .NOT.LQUERY ) THEN
189+
IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY )
190+
$ INFO = -7
190191
END IF
191192
IF( INFO.NE.0 ) THEN
192193
CALL XERBLA( 'SGELQF', -INFO )
193194
RETURN
194195
ELSE IF( LQUERY ) THEN
196+
IF( K.EQ.0 ) THEN
197+
LWKOPT = 1
198+
ELSE
199+
LWKOPT = M*NB
200+
END IF
201+
WORK( 1 ) = SROUNDUP_LWORK( LWKOPT )
195202
RETURN
196203
END IF
197204
*
198205
* Quick return if possible
199206
*
200-
K = MIN( M, N )
201207
IF( K.EQ.0 ) THEN
202208
WORK( 1 ) = 1
203209
RETURN
@@ -267,7 +273,7 @@ SUBROUTINE SGELQF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
267273
$ CALL SGELQ2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK,
268274
$ IINFO )
269275
*
270-
WORK( 1 ) = SROUNDUP_LWORK(IWS)
276+
WORK( 1 ) = SROUNDUP_LWORK( IWS )
271277
RETURN
272278
*
273279
* End of SGELQF

‎SRC/sgemlq.f‎

Lines changed: 16 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -110,13 +110,14 @@
110110
*>
111111
*> \param[out] WORK
112112
*> \verbatim
113-
*> (workspace) REAL array, dimension (MAX(1,LWORK))
113+
*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
114+
*> On exit, if INFO = 0, WORK(1) returns the minimal LWORK.
114115
*> \endverbatim
115116
*>
116117
*> \param[in] LWORK
117118
*> \verbatim
118119
*> LWORK is INTEGER
119-
*> The dimension of the array WORK.
120+
*> The dimension of the array WORK. LWORK >= 1.
120121
*> If LWORK = -1, then a workspace query is assumed. The routine
121122
*> only calculates the size of the WORK array, returns this
122123
*> value as WORK(1), and no error message related to WORK
@@ -187,7 +188,7 @@ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
187188
* ..
188189
* .. Local Scalars ..
189190
LOGICAL LEFT, RIGHT, TRAN, NOTRAN, LQUERY
190-
INTEGER MB, NB, LW, NBLCKS, MN
191+
INTEGER MB, NB, LW, NBLCKS, MN, MINMNK, LWMIN
191192
* ..
192193
* .. External Functions ..
193194
LOGICAL LSAME
@@ -207,7 +208,7 @@ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
207208
*
208209
* Test the input arguments
209210
*
210-
LQUERY = LWORK.EQ.-1
211+
LQUERY = ( LWORK.EQ.-1 )
211212
NOTRAN = LSAME( TRANS, 'N' )
212213
TRAN = LSAME( TRANS, 'T' )
213214
LEFT = LSAME( SIDE, 'L' )
@@ -222,6 +223,13 @@ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
222223
LW = M * MB
223224
MN = N
224225
END IF
226+
*
227+
MINMNK = MIN( M, N, K )
228+
IF( MINMNK.EQ.0 ) THEN
229+
LWMIN = 1
230+
ELSE
231+
LWMIN = MAX( 1, LW )
232+
END IF
225233
*
226234
IF( ( NB.GT.K ) .AND. ( MN.GT.K ) ) THEN
227235
IF( MOD( MN - K, NB - K ) .EQ. 0 ) THEN
@@ -250,12 +258,12 @@ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
250258
INFO = -9
251259
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
252260
INFO = -11
253-
ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
261+
ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
254262
INFO = -13
255263
END IF
256264
*
257265
IF( INFO.EQ.0 ) THEN
258-
WORK( 1 ) = SROUNDUP_LWORK( LW )
266+
WORK( 1 ) = SROUNDUP_LWORK( LWMIN )
259267
END IF
260268
*
261269
IF( INFO.NE.0 ) THEN
@@ -267,7 +275,7 @@ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
267275
*
268276
* Quick return if possible
269277
*
270-
IF( MIN( M, N, K ).EQ.0 ) THEN
278+
IF( MINMNK.EQ.0 ) THEN
271279
RETURN
272280
END IF
273281
*
@@ -280,7 +288,7 @@ SUBROUTINE SGEMLQ( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
280288
$ MB, C, LDC, WORK, LWORK, INFO )
281289
END IF
282290
*
283-
WORK( 1 ) = SROUNDUP_LWORK( LW )
291+
WORK( 1 ) = SROUNDUP_LWORK( LWMIN )
284292
*
285293
RETURN
286294
*

‎SRC/sgemqr.f‎

Lines changed: 15 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -189,12 +189,13 @@ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
189189
* ..
190190
* .. Local Scalars ..
191191
LOGICAL LEFT, RIGHT, TRAN, NOTRAN, LQUERY
192-
INTEGER MB, NB, LW, NBLCKS, MN
192+
INTEGER MB, NB, LW, NBLCKS, MN, MINMNK, LWMIN
193193
* ..
194194
* .. External Functions ..
195195
LOGICAL LSAME
196+
EXTERNAL LSAME
196197
REAL SROUNDUP_LWORK
197-
EXTERNAL LSAME, SROUNDUP_LWORK
198+
EXTERNAL SROUNDUP_LWORK
198199
* ..
199200
* .. External Subroutines ..
200201
EXTERNAL SGEMQRT, SLAMTSQR, XERBLA
@@ -206,7 +207,7 @@ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
206207
*
207208
* Test the input arguments
208209
*
209-
LQUERY = LWORK.EQ.-1
210+
LQUERY = ( LWORK.EQ.-1 )
210211
NOTRAN = LSAME( TRANS, 'N' )
211212
TRAN = LSAME( TRANS, 'T' )
212213
LEFT = LSAME( SIDE, 'L' )
@@ -221,6 +222,13 @@ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
221222
LW = MB * NB
222223
MN = N
223224
END IF
225+
*
226+
MINMNK = MIN( M, N, K )
227+
IF( MINMNK.EQ.0 ) THEN
228+
LWMIN = 1
229+
ELSE
230+
LWMIN = MAX( 1, LW )
231+
END IF
224232
*
225233
IF( ( MB.GT.K ) .AND. ( MN.GT.K ) ) THEN
226234
IF( MOD( MN - K, MB - K ).EQ.0 ) THEN
@@ -249,12 +257,12 @@ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
249257
INFO = -9
250258
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
251259
INFO = -11
252-
ELSE IF( ( LWORK.LT.MAX( 1, LW ) ) .AND. ( .NOT.LQUERY ) ) THEN
260+
ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
253261
INFO = -13
254262
END IF
255263
*
256264
IF( INFO.EQ.0 ) THEN
257-
WORK( 1 ) = SROUNDUP_LWORK(LW)
265+
WORK( 1 ) = SROUNDUP_LWORK( LWMIN )
258266
END IF
259267
*
260268
IF( INFO.NE.0 ) THEN
@@ -266,7 +274,7 @@ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
266274
*
267275
* Quick return if possible
268276
*
269-
IF( MIN( M, N, K ).EQ.0 ) THEN
277+
IF( MINMNK.EQ.0 ) THEN
270278
RETURN
271279
END IF
272280
*
@@ -279,7 +287,7 @@ SUBROUTINE SGEMQR( SIDE, TRANS, M, N, K, A, LDA, T, TSIZE,
279287
$ NB, C, LDC, WORK, LWORK, INFO )
280288
END IF
281289
*
282-
WORK( 1 ) = SROUNDUP_LWORK(LW)
290+
WORK( 1 ) = SROUNDUP_LWORK( LWMIN )
283291
*
284292
RETURN
285293
*

‎SRC/sgeqlf.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
*>
@@ -189,8 +190,9 @@ SUBROUTINE SGEQLF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
189190
END IF
190191
WORK( 1 ) = SROUNDUP_LWORK(LWKOPT)
191192
*
192-
IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
193-
INFO = -7
193+
IF( .NOT.LQUERY ) THEN
194+
IF( LWORK.LE.0 .OR. ( M.GT.0 .AND. LWORK.LT.MAX( 1, N ) ) )
195+
$ INFO = -7
194196
END IF
195197
END IF
196198
*

0 commit comments

Comments
 (0)