162
163
164
165
166
167
168 CHARACTER TRANS
169 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
170
171
172 REAL A( LDA, * ), B( LDB, * ), WORK( * )
173
174
175
176
177
178
179 REAL ZERO, ONE
180 parameter( zero = 0.0e0, one = 1.0e0 )
181
182
183 LOGICAL LQUERY, TRAN
184 INTEGER I, IASCL, IBSCL, J, MAXMN, BROW,
185 $ SCLLEN, TSZO, TSZM, LWO, LWM, LW1, LW2,
186 $ WSIZEO, WSIZEM, INFO2
187 REAL ANRM, BIGNUM, BNRM, SMLNUM, TQ( 5 ), WORKQ( 1 )
188
189
190 LOGICAL LSAME
191 REAL SLAMCH, SLANGE
193
194
197
198
199 INTRINSIC real, max, min, int
200
201
202
203
204
205 info = 0
206 maxmn = max( m, n )
207 tran =
lsame( trans,
'T' )
208
209 lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
210 IF( .NOT.(
lsame( trans,
'N' ) .OR.
211 $
lsame( trans,
'T' ) ) )
THEN
212 info = -1
213 ELSE IF( m.LT.0 ) THEN
214 info = -2
215 ELSE IF( n.LT.0 ) THEN
216 info = -3
217 ELSE IF( nrhs.LT.0 ) THEN
218 info = -4
219 ELSE IF( lda.LT.max( 1, m ) ) THEN
220 info = -6
221 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
222 info = -8
223 END IF
224
225 IF( info.EQ.0 ) THEN
226
227
228
229 IF( m.GE.n ) THEN
230 CALL sgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
231 tszo = int( tq( 1 ) )
232 lwo = int( workq( 1 ) )
233 CALL sgemqr(
'L', trans, m, nrhs, n, a, lda, tq,
234 $ tszo, b, ldb, workq, -1, info2 )
235 lwo = max( lwo, int( workq( 1 ) ) )
236 CALL sgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
237 tszm = int( tq( 1 ) )
238 lwm = int( workq( 1 ) )
239 CALL sgemqr(
'L', trans, m, nrhs, n, a, lda, tq,
240 $ tszm, b, ldb, workq, -1, info2 )
241 lwm = max( lwm, int( workq( 1 ) ) )
242 wsizeo = tszo + lwo
243 wsizem = tszm + lwm
244 ELSE
245 CALL sgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
246 tszo = int( tq( 1 ) )
247 lwo = int( workq( 1 ) )
248 CALL sgemlq(
'L', trans, n, nrhs, m, a, lda, tq,
249 $ tszo, b, ldb, workq, -1, info2 )
250 lwo = max( lwo, int( workq( 1 ) ) )
251 CALL sgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
252 tszm = int( tq( 1 ) )
253 lwm = int( workq( 1 ) )
254 CALL sgemlq(
'L', trans, n, nrhs, m, a, lda, tq,
255 $ tszm, b, ldb, workq, -1, info2 )
256 lwm = max( lwm, int( workq( 1 ) ) )
257 wsizeo = tszo + lwo
258 wsizem = tszm + lwm
259 END IF
260
261 IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
262 info = -10
263 END IF
264
265 work( 1 ) = real( wsizeo )
266
267 END IF
268
269 IF( info.NE.0 ) THEN
270 CALL xerbla(
'SGETSLS', -info )
271 RETURN
272 END IF
273 IF( lquery ) THEN
274 IF( lwork.EQ.-2 ) work( 1 ) = real( wsizem )
275 RETURN
276 END IF
277 IF( lwork.LT.wsizeo ) THEN
278 lw1 = tszm
279 lw2 = lwm
280 ELSE
281 lw1 = tszo
282 lw2 = lwo
283 END IF
284
285
286
287 IF( min( m, n, nrhs ).EQ.0 ) THEN
288 CALL slaset(
'FULL', max( m, n ), nrhs, zero, zero,
289 $ b, ldb )
290 RETURN
291 END IF
292
293
294
296 bignum = one / smlnum
297 CALL slabad( smlnum, bignum )
298
299
300
301 anrm =
slange(
'M', m, n, a, lda, work )
302 iascl = 0
303 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
304
305
306
307 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
308 iascl = 1
309 ELSE IF( anrm.GT.bignum ) THEN
310
311
312
313 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
314 iascl = 2
315 ELSE IF( anrm.EQ.zero ) THEN
316
317
318
319 CALL slaset(
'F', maxmn, nrhs, zero, zero, b, ldb )
320 GO TO 50
321 END IF
322
323 brow = m
324 IF ( tran ) THEN
325 brow = n
326 END IF
327 bnrm =
slange(
'M', brow, nrhs, b, ldb, work )
328 ibscl = 0
329 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
330
331
332
333 CALL slascl(
'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
334 $ info )
335 ibscl = 1
336 ELSE IF( bnrm.GT.bignum ) THEN
337
338
339
340 CALL slascl(
'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
341 $ info )
342 ibscl = 2
343 END IF
344
345 IF ( m.GE.n ) THEN
346
347
348
349 CALL sgeqr( m, n, a, lda, work( lw2+1 ), lw1,
350 $ work( 1 ), lw2, info )
351 IF ( .NOT.tran ) THEN
352
353
354
355
356
357 CALL sgemqr(
'L' ,
'T', m, nrhs, n, a, lda,
358 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
359 $ info )
360
361
362
363 CALL strtrs(
'U',
'N',
'N', n, nrhs,
364 $ a, lda, b, ldb, info )
365 IF( info.GT.0 ) THEN
366 RETURN
367 END IF
368 scllen = n
369 ELSE
370
371
372
373
374
375 CALL strtrs(
'U',
'T',
'N', n, nrhs,
376 $ a, lda, b, ldb, info )
377
378 IF( info.GT.0 ) THEN
379 RETURN
380 END IF
381
382
383
384 DO 20 j = 1, nrhs
385 DO 10 i = n + 1, m
386 b( i, j ) = zero
387 10 CONTINUE
388 20 CONTINUE
389
390
391
392 CALL sgemqr(
'L',
'N', m, nrhs, n, a, lda,
393 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
394 $ info )
395
396 scllen = m
397
398 END IF
399
400 ELSE
401
402
403
404 CALL sgelq( m, n, a, lda, work( lw2+1 ), lw1,
405 $ work( 1 ), lw2, info )
406
407
408
409 IF( .NOT.tran ) THEN
410
411
412
413
414
415 CALL strtrs(
'L',
'N',
'N', m, nrhs,
416 $ a, lda, b, ldb, info )
417
418 IF( info.GT.0 ) THEN
419 RETURN
420 END IF
421
422
423
424 DO 40 j = 1, nrhs
425 DO 30 i = m + 1, n
426 b( i, j ) = zero
427 30 CONTINUE
428 40 CONTINUE
429
430
431
432 CALL sgemlq(
'L',
'T', n, nrhs, m, a, lda,
433 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
434 $ info )
435
436
437
438 scllen = n
439
440 ELSE
441
442
443
444
445
446 CALL sgemlq(
'L',
'N', n, nrhs, m, a, lda,
447 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
448 $ info )
449
450
451
452
453
454 CALL strtrs(
'Lower',
'Transpose',
'Non-unit', m, nrhs,
455 $ a, lda, b, ldb, info )
456
457 IF( info.GT.0 ) THEN
458 RETURN
459 END IF
460
461 scllen = m
462
463 END IF
464
465 END IF
466
467
468
469 IF( iascl.EQ.1 ) THEN
470 CALL slascl(
'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
471 $ info )
472 ELSE IF( iascl.EQ.2 ) THEN
473 CALL slascl(
'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
474 $ info )
475 END IF
476 IF( ibscl.EQ.1 ) THEN
477 CALL slascl(
'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
478 $ info )
479 ELSE IF( ibscl.EQ.2 ) THEN
480 CALL slascl(
'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
481 $ info )
482 END IF
483
484 50 CONTINUE
485 work( 1 ) = real( tszo + lwo )
486 RETURN
487
488
489
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine strtrs(UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, INFO)
STRTRS
subroutine sgelq(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
SGELQ
subroutine sgemlq(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
SGEMLQ
subroutine sgemqr(SIDE, TRANS, M, N, K, A, LDA, T, TSIZE, C, LDC, WORK, LWORK, INFO)
SGEMQR
subroutine sgeqr(M, N, A, LDA, T, TSIZE, WORK, LWORK, INFO)
SGEQR
real function slamch(CMACH)
SLAMCH