226 INTEGER info, lda, ldu, ldvt, lwork, m, n
230 REAL a( lda, * ), s( * ), u( ldu, * ),
231 $ vt( ldvt, * ), work( * )
238 parameter( zero = 0.0e0, one = 1.0e0 )
241 LOGICAL lquery, wntqa, wntqas, wntqn, wntqo, wntqs
242 INTEGER bdspac, blk, chunk, i, ie, ierr, il,
243 $ ir, iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
244 $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
245 $ mnthr, nwork, wrkbl
246 REAL anrm, bignum, eps, smlnum
264 INTRINSIC int, max, min, sqrt
272 wntqa =
lsame( jobz,
'A' )
273 wntqs =
lsame( jobz,
'S' )
274 wntqas = wntqa .OR. wntqs
275 wntqo =
lsame( jobz,
'O' )
276 wntqn =
lsame( jobz,
'N' )
277 lquery = ( lwork.EQ.-1 )
279 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) )
THEN
281 ELSE IF( m.LT.0 )
THEN
283 ELSE IF( n.LT.0 )
THEN
285 ELSE IF( lda.LT.max( 1, m ) )
THEN
287 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
288 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) )
THEN
290 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
291 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
292 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) )
THEN
306 IF( m.GE.n .AND. minmn.GT.0 )
THEN
310 mnthr = int( minmn*11.0e0 / 6.0e0 )
316 IF( m.GE.mnthr )
THEN
321 wrkbl = n + n*
ilaenv( 1,
'SGEQRF',
' ', m, n, -1,
323 wrkbl = max( wrkbl, 3*n+2*n*
324 $
ilaenv( 1,
'SGEBRD',
' ', n, n, -1, -1 ) )
325 maxwrk = max( wrkbl, bdspac+n )
327 ELSE IF( wntqo )
THEN
331 wrkbl = n + n*
ilaenv( 1,
'SGEQRF',
' ', m, n, -1, -1 )
332 wrkbl = max( wrkbl, n+n*
ilaenv( 1,
'SORGQR',
' ', m,
334 wrkbl = max( wrkbl, 3*n+2*n*
335 $
ilaenv( 1,
'SGEBRD',
' ', n, n, -1, -1 ) )
336 wrkbl = max( wrkbl, 3*n+n*
337 $
ilaenv( 1,
'SORMBR',
'QLN', n, n, n, -1 ) )
338 wrkbl = max( wrkbl, 3*n+n*
339 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, n, -1 ) )
340 wrkbl = max( wrkbl, bdspac+3*n )
341 maxwrk = wrkbl + 2*n*n
342 minwrk = bdspac + 2*n*n + 3*n
343 ELSE IF( wntqs )
THEN
347 wrkbl = n + n*
ilaenv( 1,
'SGEQRF',
' ', m, n, -1, -1 )
348 wrkbl = max( wrkbl, n+n*
ilaenv( 1,
'SORGQR',
' ', m,
350 wrkbl = max( wrkbl, 3*n+2*n*
351 $
ilaenv( 1,
'SGEBRD',
' ', n, n, -1, -1 ) )
352 wrkbl = max( wrkbl, 3*n+n*
353 $
ilaenv( 1,
'SORMBR',
'QLN', n, n, n, -1 ) )
354 wrkbl = max( wrkbl, 3*n+n*
355 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, n, -1 ) )
356 wrkbl = max( wrkbl, bdspac+3*n )
358 minwrk = bdspac + n*n + 3*n
359 ELSE IF( wntqa )
THEN
363 wrkbl = n + n*
ilaenv( 1,
'SGEQRF',
' ', m, n, -1, -1 )
364 wrkbl = max( wrkbl, n+m*
ilaenv( 1,
'SORGQR',
' ', m,
366 wrkbl = max( wrkbl, 3*n+2*n*
367 $
ilaenv( 1,
'SGEBRD',
' ', n, n, -1, -1 ) )
368 wrkbl = max( wrkbl, 3*n+n*
369 $
ilaenv( 1,
'SORMBR',
'QLN', n, n, n, -1 ) )
370 wrkbl = max( wrkbl, 3*n+n*
371 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, n, -1 ) )
372 wrkbl = max( wrkbl, bdspac+3*n )
374 minwrk = bdspac + n*n + 2*n + m
380 wrkbl = 3*n + ( m+n )*
ilaenv( 1,
'SGEBRD',
' ', m, n, -1,
383 maxwrk = max( wrkbl, bdspac+3*n )
384 minwrk = 3*n + max( m, bdspac )
385 ELSE IF( wntqo )
THEN
386 wrkbl = max( wrkbl, 3*n+n*
387 $
ilaenv( 1,
'SORMBR',
'QLN', m, n, n, -1 ) )
388 wrkbl = max( wrkbl, 3*n+n*
389 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, n, -1 ) )
390 wrkbl = max( wrkbl, bdspac+3*n )
392 minwrk = 3*n + max( m, n*n+bdspac )
393 ELSE IF( wntqs )
THEN
394 wrkbl = max( wrkbl, 3*n+n*
395 $
ilaenv( 1,
'SORMBR',
'QLN', m, n, n, -1 ) )
396 wrkbl = max( wrkbl, 3*n+n*
397 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, n, -1 ) )
398 maxwrk = max( wrkbl, bdspac+3*n )
399 minwrk = 3*n + max( m, bdspac )
400 ELSE IF( wntqa )
THEN
401 wrkbl = max( wrkbl, 3*n+m*
402 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, n, -1 ) )
403 wrkbl = max( wrkbl, 3*n+n*
404 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, n, -1 ) )
405 maxwrk = max( maxwrk, bdspac+3*n )
406 minwrk = 3*n + max( m, bdspac )
409 ELSE IF ( minmn.GT.0 )
THEN
413 mnthr = int( minmn*11.0e0 / 6.0e0 )
419 IF( n.GE.mnthr )
THEN
424 wrkbl = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1,
426 wrkbl = max( wrkbl, 3*m+2*m*
427 $
ilaenv( 1,
'SGEBRD',
' ', m, m, -1, -1 ) )
428 maxwrk = max( wrkbl, bdspac+m )
430 ELSE IF( wntqo )
THEN
434 wrkbl = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1, -1 )
435 wrkbl = max( wrkbl, m+m*
ilaenv( 1,
'SORGLQ',
' ', m,
437 wrkbl = max( wrkbl, 3*m+2*m*
438 $
ilaenv( 1,
'SGEBRD',
' ', m, m, -1, -1 ) )
439 wrkbl = max( wrkbl, 3*m+m*
440 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, m, -1 ) )
441 wrkbl = max( wrkbl, 3*m+m*
442 $
ilaenv( 1,
'SORMBR',
'PRT', m, m, m, -1 ) )
443 wrkbl = max( wrkbl, bdspac+3*m )
444 maxwrk = wrkbl + 2*m*m
445 minwrk = bdspac + 2*m*m + 3*m
446 ELSE IF( wntqs )
THEN
450 wrkbl = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1, -1 )
451 wrkbl = max( wrkbl, m+m*
ilaenv( 1,
'SORGLQ',
' ', m,
453 wrkbl = max( wrkbl, 3*m+2*m*
454 $
ilaenv( 1,
'SGEBRD',
' ', m, m, -1, -1 ) )
455 wrkbl = max( wrkbl, 3*m+m*
456 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, m, -1 ) )
457 wrkbl = max( wrkbl, 3*m+m*
458 $
ilaenv( 1,
'SORMBR',
'PRT', m, m, m, -1 ) )
459 wrkbl = max( wrkbl, bdspac+3*m )
461 minwrk = bdspac + m*m + 3*m
462 ELSE IF( wntqa )
THEN
466 wrkbl = m + m*
ilaenv( 1,
'SGELQF',
' ', m, n, -1, -1 )
467 wrkbl = max( wrkbl, m+n*
ilaenv( 1,
'SORGLQ',
' ', n,
469 wrkbl = max( wrkbl, 3*m+2*m*
470 $
ilaenv( 1,
'SGEBRD',
' ', m, m, -1, -1 ) )
471 wrkbl = max( wrkbl, 3*m+m*
472 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, m, -1 ) )
473 wrkbl = max( wrkbl, 3*m+m*
474 $
ilaenv( 1,
'SORMBR',
'PRT', m, m, m, -1 ) )
475 wrkbl = max( wrkbl, bdspac+3*m )
477 minwrk = bdspac + m*m + 3*m
483 wrkbl = 3*m + ( m+n )*
ilaenv( 1,
'SGEBRD',
' ', m, n, -1,
486 maxwrk = max( wrkbl, bdspac+3*m )
487 minwrk = 3*m + max( n, bdspac )
488 ELSE IF( wntqo )
THEN
489 wrkbl = max( wrkbl, 3*m+m*
490 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, n, -1 ) )
491 wrkbl = max( wrkbl, 3*m+m*
492 $
ilaenv( 1,
'SORMBR',
'PRT', m, n, m, -1 ) )
493 wrkbl = max( wrkbl, bdspac+3*m )
495 minwrk = 3*m + max( n, m*m+bdspac )
496 ELSE IF( wntqs )
THEN
497 wrkbl = max( wrkbl, 3*m+m*
498 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, n, -1 ) )
499 wrkbl = max( wrkbl, 3*m+m*
500 $
ilaenv( 1,
'SORMBR',
'PRT', m, n, m, -1 ) )
501 maxwrk = max( wrkbl, bdspac+3*m )
502 minwrk = 3*m + max( n, bdspac )
503 ELSE IF( wntqa )
THEN
504 wrkbl = max( wrkbl, 3*m+m*
505 $
ilaenv( 1,
'SORMBR',
'QLN', m, m, n, -1 ) )
506 wrkbl = max( wrkbl, 3*m+m*
507 $
ilaenv( 1,
'SORMBR',
'PRT', n, n, m, -1 ) )
508 maxwrk = max( wrkbl, bdspac+3*m )
509 minwrk = 3*m + max( n, bdspac )
513 maxwrk = max( maxwrk, minwrk )
516 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
522 CALL xerbla(
'SGESDD', -info )
524 ELSE IF( lquery )
THEN
530 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
537 smlnum = sqrt(
slamch(
'S' ) ) / eps
538 bignum = one / smlnum
542 anrm =
slange(
'M', m, n, a, lda, dum )
544 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
546 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
547 ELSE IF( anrm.GT.bignum )
THEN
549 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
558 IF( m.GE.mnthr )
THEN
571 CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
572 $ lwork-nwork+1, ierr )
576 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
585 CALL sgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
586 $ work( itaup ), work( nwork ), lwork-nwork+1,
593 CALL sbdsdc(
'U',
'N', n, s, work( ie ), dum, 1, dum, 1,
594 $ dum, idum, work( nwork ), iwork, info )
596 ELSE IF( wntqo )
THEN
606 IF( lwork.GE.lda*n+n*n+3*n+bdspac )
THEN
609 ldwrkr = ( lwork-n*n-3*n-bdspac ) / n
617 CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
618 $ lwork-nwork+1, ierr )
622 CALL slacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
623 CALL slaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
629 CALL sorgqr( m, n, n, a, lda, work( itau ),
630 $ work( nwork ), lwork-nwork+1, ierr )
639 CALL sgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
640 $ work( itauq ), work( itaup ), work( nwork ),
641 $ lwork-nwork+1, ierr )
653 CALL sbdsdc(
'U',
'I', n, s, work( ie ), work( iu ), n,
654 $ vt, ldvt, dum, idum, work( nwork ), iwork,
661 CALL sormbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
662 $ work( itauq ), work( iu ), n, work( nwork ),
663 $ lwork-nwork+1, ierr )
664 CALL sormbr(
'P',
'R',
'T', n, n, n, work( ir ), ldwrkr,
665 $ work( itaup ), vt, ldvt, work( nwork ),
666 $ lwork-nwork+1, ierr )
672 DO 10 i = 1, m, ldwrkr
673 chunk = min( m-i+1, ldwrkr )
674 CALL sgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
675 $ lda, work( iu ), n, zero, work( ir ),
677 CALL slacpy(
'F', chunk, n, work( ir ), ldwrkr,
681 ELSE IF( wntqs )
THEN
698 CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
699 $ lwork-nwork+1, ierr )
703 CALL slacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
704 CALL slaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
710 CALL sorgqr( m, n, n, a, lda, work( itau ),
711 $ work( nwork ), lwork-nwork+1, ierr )
720 CALL sgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
721 $ work( itauq ), work( itaup ), work( nwork ),
722 $ lwork-nwork+1, ierr )
729 CALL sbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
730 $ ldvt, dum, idum, work( nwork ), iwork,
737 CALL sormbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
738 $ work( itauq ), u, ldu, work( nwork ),
739 $ lwork-nwork+1, ierr )
741 CALL sormbr(
'P',
'R',
'T', n, n, n, work( ir ), ldwrkr,
742 $ work( itaup ), vt, ldvt, work( nwork ),
743 $ lwork-nwork+1, ierr )
749 CALL slacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
750 CALL sgemm(
'N',
'N', m, n, n, one, a, lda, work( ir ),
751 $ ldwrkr, zero, u, ldu )
753 ELSE IF( wntqa )
THEN
770 CALL sgeqrf( m, n, a, lda, work( itau ), work( nwork ),
771 $ lwork-nwork+1, ierr )
772 CALL slacpy(
'L', m, n, a, lda, u, ldu )
776 CALL sorgqr( m, m, n, u, ldu, work( itau ),
777 $ work( nwork ), lwork-nwork+1, ierr )
781 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
790 CALL sgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
791 $ work( itaup ), work( nwork ), lwork-nwork+1,
799 CALL sbdsdc(
'U',
'I', n, s, work( ie ), work( iu ), n,
800 $ vt, ldvt, dum, idum, work( nwork ), iwork,
807 CALL sormbr(
'Q',
'L',
'N', n, n, n, a, lda,
808 $ work( itauq ), work( iu ), ldwrku,
809 $ work( nwork ), lwork-nwork+1, ierr )
810 CALL sormbr(
'P',
'R',
'T', n, n, n, a, lda,
811 $ work( itaup ), vt, ldvt, work( nwork ),
812 $ lwork-nwork+1, ierr )
818 CALL sgemm(
'N',
'N', m, n, n, one, u, ldu, work( iu ),
819 $ ldwrku, zero, a, lda )
823 CALL slacpy(
'F', m, n, a, lda, u, ldu )
842 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
843 $ work( itaup ), work( nwork ), lwork-nwork+1,
850 CALL sbdsdc(
'U',
'N', n, s, work( ie ), dum, 1, dum, 1,
851 $ dum, idum, work( nwork ), iwork, info )
852 ELSE IF( wntqo )
THEN
854 IF( lwork.GE.m*n+3*n+bdspac )
THEN
859 nwork = iu + ldwrku*n
860 CALL slaset(
'F', m, n, zero, zero, work( iu ),
867 nwork = iu + ldwrku*n
872 ldwrkr = ( lwork-n*n-3*n ) / n
874 nwork = iu + ldwrku*n
881 CALL sbdsdc(
'U',
'I', n, s, work( ie ), work( iu ),
882 $ ldwrku, vt, ldvt, dum, idum, work( nwork ),
888 CALL sormbr(
'P',
'R',
'T', n, n, n, a, lda,
889 $ work( itaup ), vt, ldvt, work( nwork ),
890 $ lwork-nwork+1, ierr )
892 IF( lwork.GE.m*n+3*n+bdspac )
THEN
897 CALL sormbr(
'Q',
'L',
'N', m, n, n, a, lda,
898 $ work( itauq ), work( iu ), ldwrku,
899 $ work( nwork ), lwork-nwork+1, ierr )
903 CALL slacpy(
'F', m, n, work( iu ), ldwrku, a, lda )
909 CALL sorgbr(
'Q', m, n, n, a, lda, work( itauq ),
910 $ work( nwork ), lwork-nwork+1, ierr )
917 DO 20 i = 1, m, ldwrkr
918 chunk = min( m-i+1, ldwrkr )
919 CALL sgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
920 $ lda, work( iu ), ldwrku, zero,
921 $ work( ir ), ldwrkr )
922 CALL slacpy(
'F', chunk, n, work( ir ), ldwrkr,
927 ELSE IF( wntqs )
THEN
934 CALL slaset(
'F', m, n, zero, zero, u, ldu )
935 CALL sbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
936 $ ldvt, dum, idum, work( nwork ), iwork,
943 CALL sormbr(
'Q',
'L',
'N', m, n, n, a, lda,
944 $ work( itauq ), u, ldu, work( nwork ),
945 $ lwork-nwork+1, ierr )
946 CALL sormbr(
'P',
'R',
'T', n, n, n, a, lda,
947 $ work( itaup ), vt, ldvt, work( nwork ),
948 $ lwork-nwork+1, ierr )
949 ELSE IF( wntqa )
THEN
956 CALL slaset(
'F', m, m, zero, zero, u, ldu )
957 CALL sbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
958 $ ldvt, dum, idum, work( nwork ), iwork,
964 CALL slaset(
'F', m-n, m-n, zero, one, u( n+1, n+1 ),
972 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
973 $ work( itauq ), u, ldu, work( nwork ),
974 $ lwork-nwork+1, ierr )
975 CALL sormbr(
'P',
'R',
'T', n, n, m, a, lda,
976 $ work( itaup ), vt, ldvt, work( nwork ),
977 $ lwork-nwork+1, ierr )
988 IF( n.GE.mnthr )
THEN
1001 CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
1002 $ lwork-nwork+1, ierr )
1006 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1015 CALL sgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1016 $ work( itaup ), work( nwork ), lwork-nwork+1,
1023 CALL sbdsdc(
'U',
'N', m, s, work( ie ), dum, 1, dum, 1,
1024 $ dum, idum, work( nwork ), iwork, info )
1026 ELSE IF( wntqo )
THEN
1037 IF( lwork.GE.m*n+m*m+3*m+bdspac )
THEN
1045 chunk = ( lwork-m*m ) / m
1047 itau = il + ldwrkl*m
1053 CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
1054 $ lwork-nwork+1, ierr )
1058 CALL slacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1059 CALL slaset(
'U', m-1, m-1, zero, zero,
1060 $ work( il+ldwrkl ), ldwrkl )
1065 CALL sorglq( m, n, m, a, lda, work( itau ),
1066 $ work( nwork ), lwork-nwork+1, ierr )
1075 CALL sgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1076 $ work( itauq ), work( itaup ), work( nwork ),
1077 $ lwork-nwork+1, ierr )
1084 CALL sbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1085 $ work( ivt ), m, dum, idum, work( nwork ),
1092 CALL sormbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1093 $ work( itauq ), u, ldu, work( nwork ),
1094 $ lwork-nwork+1, ierr )
1095 CALL sormbr(
'P',
'R',
'T', m, m, m, work( il ), ldwrkl,
1096 $ work( itaup ), work( ivt ), m,
1097 $ work( nwork ), lwork-nwork+1, ierr )
1103 DO 30 i = 1, n, chunk
1104 blk = min( n-i+1, chunk )
1105 CALL sgemm(
'N',
'N', m, blk, m, one, work( ivt ), m,
1106 $ a( 1, i ), lda, zero, work( il ), ldwrkl )
1107 CALL slacpy(
'F', m, blk, work( il ), ldwrkl,
1111 ELSE IF( wntqs )
THEN
1122 itau = il + ldwrkl*m
1128 CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
1129 $ lwork-nwork+1, ierr )
1133 CALL slacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1134 CALL slaset(
'U', m-1, m-1, zero, zero,
1135 $ work( il+ldwrkl ), ldwrkl )
1140 CALL sorglq( m, n, m, a, lda, work( itau ),
1141 $ work( nwork ), lwork-nwork+1, ierr )
1150 CALL sgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1151 $ work( itauq ), work( itaup ), work( nwork ),
1152 $ lwork-nwork+1, ierr )
1159 CALL sbdsdc(
'U',
'I', m, s, work( ie ), u, ldu, vt,
1160 $ ldvt, dum, idum, work( nwork ), iwork,
1167 CALL sormbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1168 $ work( itauq ), u, ldu, work( nwork ),
1169 $ lwork-nwork+1, ierr )
1170 CALL sormbr(
'P',
'R',
'T', m, m, m, work( il ), ldwrkl,
1171 $ work( itaup ), vt, ldvt, work( nwork ),
1172 $ lwork-nwork+1, ierr )
1178 CALL slacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1179 CALL sgemm(
'N',
'N', m, n, m, one, work( il ), ldwrkl,
1180 $ a, lda, zero, vt, ldvt )
1182 ELSE IF( wntqa )
THEN
1193 itau = ivt + ldwkvt*m
1199 CALL sgelqf( m, n, a, lda, work( itau ), work( nwork ),
1200 $ lwork-nwork+1, ierr )
1201 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
1206 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
1207 $ work( nwork ), lwork-nwork+1, ierr )
1211 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1220 CALL sgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1221 $ work( itaup ), work( nwork ), lwork-nwork+1,
1229 CALL sbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1230 $ work( ivt ), ldwkvt, dum, idum,
1231 $ work( nwork ), iwork, info )
1237 CALL sormbr(
'Q',
'L',
'N', m, m, m, a, lda,
1238 $ work( itauq ), u, ldu, work( nwork ),
1239 $ lwork-nwork+1, ierr )
1240 CALL sormbr(
'P',
'R',
'T', m, m, m, a, lda,
1241 $ work( itaup ), work( ivt ), ldwkvt,
1242 $ work( nwork ), lwork-nwork+1, ierr )
1248 CALL sgemm(
'N',
'N', m, n, m, one, work( ivt ), ldwkvt,
1249 $ vt, ldvt, zero, a, lda )
1253 CALL slacpy(
'F', m, n, a, lda, vt, ldvt )
1272 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1273 $ work( itaup ), work( nwork ), lwork-nwork+1,
1280 CALL sbdsdc(
'L',
'N', m, s, work( ie ), dum, 1, dum, 1,
1281 $ dum, idum, work( nwork ), iwork, info )
1282 ELSE IF( wntqo )
THEN
1285 IF( lwork.GE.m*n+3*m+bdspac )
THEN
1289 CALL slaset(
'F', m, n, zero, zero, work( ivt ),
1291 nwork = ivt + ldwkvt*n
1296 nwork = ivt + ldwkvt*m
1301 chunk = ( lwork-m*m-3*m ) / m
1309 CALL sbdsdc(
'L',
'I', m, s, work( ie ), u, ldu,
1310 $ work( ivt ), ldwkvt, dum, idum,
1311 $ work( nwork ), iwork, info )
1316 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1317 $ work( itauq ), u, ldu, work( nwork ),
1318 $ lwork-nwork+1, ierr )
1320 IF( lwork.GE.m*n+3*m+bdspac )
THEN
1325 CALL sormbr(
'P',
'R',
'T', m, n, m, a, lda,
1326 $ work( itaup ), work( ivt ), ldwkvt,
1327 $ work( nwork ), lwork-nwork+1, ierr )
1331 CALL slacpy(
'F', m, n, work( ivt ), ldwkvt, a, lda )
1337 CALL sorgbr(
'P', m, n, m, a, lda, work( itaup ),
1338 $ work( nwork ), lwork-nwork+1, ierr )
1345 DO 40 i = 1, n, chunk
1346 blk = min( n-i+1, chunk )
1347 CALL sgemm(
'N',
'N', m, blk, m, one, work( ivt ),
1348 $ ldwkvt, a( 1, i ), lda, zero,
1350 CALL slacpy(
'F', m, blk, work( il ), m, a( 1, i ),
1354 ELSE IF( wntqs )
THEN
1361 CALL slaset(
'F', m, n, zero, zero, vt, ldvt )
1362 CALL sbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1363 $ ldvt, dum, idum, work( nwork ), iwork,
1370 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1371 $ work( itauq ), u, ldu, work( nwork ),
1372 $ lwork-nwork+1, ierr )
1373 CALL sormbr(
'P',
'R',
'T', m, n, m, a, lda,
1374 $ work( itaup ), vt, ldvt, work( nwork ),
1375 $ lwork-nwork+1, ierr )
1376 ELSE IF( wntqa )
THEN
1383 CALL slaset(
'F', n, n, zero, zero, vt, ldvt )
1384 CALL sbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1385 $ ldvt, dum, idum, work( nwork ), iwork,
1391 CALL slaset(
'F', n-m, n-m, zero, one, vt( m+1, m+1 ),
1399 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1400 $ work( itauq ), u, ldu, work( nwork ),
1401 $ lwork-nwork+1, ierr )
1402 CALL sormbr(
'P',
'R',
'T', n, n, m, a, lda,
1403 $ work( itaup ), vt, ldvt, work( nwork ),
1404 $ lwork-nwork+1, ierr )
1413 IF( iscl.EQ.1 )
THEN
1414 IF( anrm.GT.bignum )
1415 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
1417 IF( anrm.LT.smlnum )
1418 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
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...
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 sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
real function slamch(CMACH)
SLAMCH
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMBR
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ