210
211
212
213
214
215
216 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
217 REAL RCOND
218
219
220 INTEGER IWORK( * )
221 REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
222
223
224
225
226
227 REAL ZERO, ONE, TWO
228 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
229
230
231 LOGICAL LQUERY
232 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
233 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
234 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
235 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
236
237
240
241
242 INTEGER ILAENV
243 REAL SLAMCH, SLANGE
245
246
247 INTRINSIC int, log, max, min, real
248
249
250
251
252
253 info = 0
254 minmn = min( m, n )
255 maxmn = max( m, n )
256 lquery = ( lwork.EQ.-1 )
257 IF( m.LT.0 ) THEN
258 info = -1
259 ELSE IF( n.LT.0 ) THEN
260 info = -2
261 ELSE IF( nrhs.LT.0 ) THEN
262 info = -3
263 ELSE IF( lda.LT.max( 1, m ) ) THEN
264 info = -5
265 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
266 info = -7
267 END IF
268
269
270
271
272
273
274
275
276 IF( info.EQ.0 ) THEN
277 minwrk = 1
278 maxwrk = 1
279 liwork = 1
280 IF( minmn.GT.0 ) THEN
281 smlsiz =
ilaenv( 9,
'SGELSD',
' ', 0, 0, 0, 0 )
282 mnthr =
ilaenv( 6,
'SGELSD',
' ', m, n, nrhs, -1 )
283 nlvl = max( int( log( real( minmn ) / real( smlsiz + 1 ) ) /
284 $ log( two ) ) + 1, 0 )
285 liwork = 3*minmn*nlvl + 11*minmn
286 mm = m
287 IF( m.GE.n .AND. m.GE.mnthr ) THEN
288
289
290
291
292 mm = n
293 maxwrk = max( maxwrk, n + n*
ilaenv( 1,
'SGEQRF',
' ', m,
294 $ n, -1, -1 ) )
295 maxwrk = max( maxwrk, n + nrhs*
ilaenv( 1,
'SORMQR',
'LT',
296 $ m, nrhs, n, -1 ) )
297 END IF
298 IF( m.GE.n ) THEN
299
300
301
302 maxwrk = max( maxwrk, 3*n + ( mm + n )*
ilaenv( 1,
303 $ 'SGEBRD', ' ', mm, n, -1, -1 ) )
304 maxwrk = max( maxwrk, 3*n + nrhs*
ilaenv( 1,
'SORMBR',
305 $ 'QLT', mm, nrhs, n, -1 ) )
306 maxwrk = max( maxwrk, 3*n + ( n - 1 )*
ilaenv( 1,
307 $ 'SORMBR', 'PLN', n, nrhs, n, -1 ) )
308 wlalsd = 9*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs +
309 $ ( smlsiz + 1 )**2
310 maxwrk = max( maxwrk, 3*n + wlalsd )
311 minwrk = max( 3*n + mm, 3*n + nrhs, 3*n + wlalsd )
312 END IF
313 IF( n.GT.m ) THEN
314 wlalsd = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs +
315 $ ( smlsiz + 1 )**2
316 IF( n.GE.mnthr ) THEN
317
318
319
320
321 maxwrk = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1,
322 $ -1 )
323 maxwrk = max( maxwrk, m*m + 4*m + 2*m*
ilaenv( 1,
324 $ 'SGEBRD', ' ', m, m, -1, -1 ) )
325 maxwrk = max( maxwrk, m*m + 4*m + nrhs*
ilaenv( 1,
326 $ 'SORMBR', 'QLT', m, nrhs, m, -1 ) )
327 maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*
ilaenv( 1,
328 $ 'SORMBR', 'PLN', m, nrhs, m, -1 ) )
329 IF( nrhs.GT.1 ) THEN
330 maxwrk = max( maxwrk, m*m + m + m*nrhs )
331 ELSE
332 maxwrk = max( maxwrk, m*m + 2*m )
333 END IF
334 maxwrk = max( maxwrk, m + nrhs*
ilaenv( 1,
'SORMLQ',
335 $ 'LT', n, nrhs, m, -1 ) )
336 maxwrk = max( maxwrk, m*m + 4*m + wlalsd )
337
338
339 maxwrk = max( maxwrk,
340 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
341 ELSE
342
343
344
345 maxwrk = 3*m + ( n + m )*
ilaenv( 1,
'SGEBRD',
' ', m,
346 $ n, -1, -1 )
347 maxwrk = max( maxwrk, 3*m + nrhs*
ilaenv( 1,
'SORMBR',
348 $ 'QLT', m, nrhs, n, -1 ) )
349 maxwrk = max( maxwrk, 3*m + m*
ilaenv( 1,
'SORMBR',
350 $ 'PLN', n, nrhs, m, -1 ) )
351 maxwrk = max( maxwrk, 3*m + wlalsd )
352 END IF
353 minwrk = max( 3*m + nrhs, 3*m + m, 3*m + wlalsd )
354 END IF
355 END IF
356 minwrk = min( minwrk, maxwrk )
357 work( 1 ) = maxwrk
358 iwork( 1 ) = liwork
359
360 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
361 info = -12
362 END IF
363 END IF
364
365 IF( info.NE.0 ) THEN
366 CALL xerbla(
'SGELSD', -info )
367 RETURN
368 ELSE IF( lquery ) THEN
369 RETURN
370 END IF
371
372
373
374 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
375 rank = 0
376 RETURN
377 END IF
378
379
380
383 smlnum = sfmin / eps
384 bignum = one / smlnum
385 CALL slabad( smlnum, bignum )
386
387
388
389 anrm =
slange(
'M', m, n, a, lda, work )
390 iascl = 0
391 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
392
393
394
395 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
396 iascl = 1
397 ELSE IF( anrm.GT.bignum ) THEN
398
399
400
401 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
402 iascl = 2
403 ELSE IF( anrm.EQ.zero ) THEN
404
405
406
407 CALL slaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
408 CALL slaset(
'F', minmn, 1, zero, zero, s, 1 )
409 rank = 0
410 GO TO 10
411 END IF
412
413
414
415 bnrm =
slange(
'M', m, nrhs, b, ldb, work )
416 ibscl = 0
417 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
418
419
420
421 CALL slascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
422 ibscl = 1
423 ELSE IF( bnrm.GT.bignum ) THEN
424
425
426
427 CALL slascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
428 ibscl = 2
429 END IF
430
431
432
433 IF( m.LT.n )
434 $
CALL slaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
435
436
437
438 IF( m.GE.n ) THEN
439
440
441
442 mm = m
443 IF( m.GE.mnthr ) THEN
444
445
446
447 mm = n
448 itau = 1
449 nwork = itau + n
450
451
452
453
454 CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
455 $ lwork-nwork+1, info )
456
457
458
459
460 CALL sormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ), b,
461 $ ldb, work( nwork ), lwork-nwork+1, info )
462
463
464
465 IF( n.GT.1 ) THEN
466 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
467 END IF
468 END IF
469
470 ie = 1
471 itauq = ie + n
472 itaup = itauq + n
473 nwork = itaup + n
474
475
476
477
478 CALL sgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
479 $ work( itaup ), work( nwork ), lwork-nwork+1,
480 $ info )
481
482
483
484
485 CALL sormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
486 $ b, ldb, work( nwork ), lwork-nwork+1, info )
487
488
489
490 CALL slalsd(
'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
491 $ rcond, rank, work( nwork ), iwork, info )
492 IF( info.NE.0 ) THEN
493 GO TO 10
494 END IF
495
496
497
498 CALL sormbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
499 $ b, ldb, work( nwork ), lwork-nwork+1, info )
500
501 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
502 $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
503
504
505
506
507 ldwork = m
508 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
509 $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
510 itau = 1
511 nwork = m + 1
512
513
514
515
516 CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
517 $ lwork-nwork+1, info )
518 il = nwork
519
520
521
522 CALL slacpy(
'L', m, m, a, lda, work( il ), ldwork )
523 CALL slaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
524 $ ldwork )
525 ie = il + ldwork*m
526 itauq = ie + m
527 itaup = itauq + m
528 nwork = itaup + m
529
530
531
532
533 CALL sgebrd( m, m, work( il ), ldwork, s, work( ie ),
534 $ work( itauq ), work( itaup ), work( nwork ),
535 $ lwork-nwork+1, info )
536
537
538
539
540 CALL sormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
541 $ work( itauq ), b, ldb, work( nwork ),
542 $ lwork-nwork+1, info )
543
544
545
546 CALL slalsd(
'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
547 $ rcond, rank, work( nwork ), iwork, info )
548 IF( info.NE.0 ) THEN
549 GO TO 10
550 END IF
551
552
553
554 CALL sormbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
555 $ work( itaup ), b, ldb, work( nwork ),
556 $ lwork-nwork+1, info )
557
558
559
560 CALL slaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
561 nwork = itau + m
562
563
564
565
566 CALL sormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ), b,
567 $ ldb, work( nwork ), lwork-nwork+1, info )
568
569 ELSE
570
571
572
573 ie = 1
574 itauq = ie + m
575 itaup = itauq + m
576 nwork = itaup + m
577
578
579
580
581 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
582 $ work( itaup ), work( nwork ), lwork-nwork+1,
583 $ info )
584
585
586
587
588 CALL sormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
589 $ b, ldb, work( nwork ), lwork-nwork+1, info )
590
591
592
593 CALL slalsd(
'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
594 $ rcond, rank, work( nwork ), iwork, info )
595 IF( info.NE.0 ) THEN
596 GO TO 10
597 END IF
598
599
600
601 CALL sormbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
602 $ b, ldb, work( nwork ), lwork-nwork+1, info )
603
604 END IF
605
606
607
608 IF( iascl.EQ.1 ) THEN
609 CALL slascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
610 CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
611 $ info )
612 ELSE IF( iascl.EQ.2 ) THEN
613 CALL slascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
614 CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
615 $ info )
616 END IF
617 IF( ibscl.EQ.1 ) THEN
618 CALL slascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
619 ELSE IF( ibscl.EQ.2 ) THEN
620 CALL slascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
621 END IF
622
623 10 CONTINUE
624 work( 1 ) = maxwrk
625 iwork( 1 ) = liwork
626 RETURN
627
628
629
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
subroutine xerbla(SRNAME, INFO)
XERBLA
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
subroutine slalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
SLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
real function slamch(CMACH)
SLAMCH