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, SROUNDUP_LWORK
193
194
197
198
199 INTRINSIC 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
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
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
298
299
300 anrm =
slange(
'M', m, n, a, lda, work )
301 iascl = 0
302 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
303
304
305
306 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
307 iascl = 1
308 ELSE IF( anrm.GT.bignum ) THEN
309
310
311
312 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
313 iascl = 2
314 ELSE IF( anrm.EQ.zero ) THEN
315
316
317
318 CALL slaset(
'F', maxmn, nrhs, zero, zero, b, ldb )
319 GO TO 50
320 END IF
321
322 brow = m
323 IF ( tran ) THEN
324 brow = n
325 END IF
326 bnrm =
slange(
'M', brow, nrhs, b, ldb, work )
327 ibscl = 0
328 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
329
330
331
332 CALL slascl(
'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
333 $ info )
334 ibscl = 1
335 ELSE IF( bnrm.GT.bignum ) THEN
336
337
338
339 CALL slascl(
'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
340 $ info )
341 ibscl = 2
342 END IF
343
344 IF ( m.GE.n ) THEN
345
346
347
348 CALL sgeqr( m, n, a, lda, work( lw2+1 ), lw1,
349 $ work( 1 ), lw2, info )
350 IF ( .NOT.tran ) THEN
351
352
353
354
355
356 CALL sgemqr(
'L' ,
'T', m, nrhs, n, a, lda,
357 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
358 $ info )
359
360
361
362 CALL strtrs(
'U',
'N',
'N', n, nrhs,
363 $ a, lda, b, ldb, info )
364 IF( info.GT.0 ) THEN
365 RETURN
366 END IF
367 scllen = n
368 ELSE
369
370
371
372
373
374 CALL strtrs(
'U',
'T',
'N', n, nrhs,
375 $ a, lda, b, ldb, info )
376
377 IF( info.GT.0 ) THEN
378 RETURN
379 END IF
380
381
382
383 DO 20 j = 1, nrhs
384 DO 10 i = n + 1, m
385 b( i, j ) = zero
386 10 CONTINUE
387 20 CONTINUE
388
389
390
391 CALL sgemqr(
'L',
'N', m, nrhs, n, a, lda,
392 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
393 $ info )
394
395 scllen = m
396
397 END IF
398
399 ELSE
400
401
402
403 CALL sgelq( m, n, a, lda, work( lw2+1 ), lw1,
404 $ work( 1 ), lw2, info )
405
406
407
408 IF( .NOT.tran ) THEN
409
410
411
412
413
414 CALL strtrs(
'L',
'N',
'N', m, nrhs,
415 $ a, lda, b, ldb, info )
416
417 IF( info.GT.0 ) THEN
418 RETURN
419 END IF
420
421
422
423 DO 40 j = 1, nrhs
424 DO 30 i = m + 1, n
425 b( i, j ) = zero
426 30 CONTINUE
427 40 CONTINUE
428
429
430
431 CALL sgemlq(
'L',
'T', n, nrhs, m, a, lda,
432 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
433 $ info )
434
435
436
437 scllen = n
438
439 ELSE
440
441
442
443
444
445 CALL sgemlq(
'L',
'N', n, nrhs, m, a, lda,
446 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
447 $ info )
448
449
450
451
452
453 CALL strtrs(
'Lower',
'Transpose',
'Non-unit', m, nrhs,
454 $ a, lda, b, ldb, info )
455
456 IF( info.GT.0 ) THEN
457 RETURN
458 END IF
459
460 scllen = m
461
462 END IF
463
464 END IF
465
466
467
468 IF( iascl.EQ.1 ) THEN
469 CALL slascl(
'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
470 $ info )
471 ELSE IF( iascl.EQ.2 ) THEN
472 CALL slascl(
'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
473 $ info )
474 END IF
475 IF( ibscl.EQ.1 ) THEN
476 CALL slascl(
'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
477 $ info )
478 ELSE IF( ibscl.EQ.2 ) THEN
479 CALL slascl(
'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
480 $ info )
481 END IF
482
483 50 CONTINUE
485 RETURN
486
487
488
subroutine xerbla(srname, info)
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
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 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.
logical function lsame(ca, cb)
LSAME
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine strtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
STRTRS