216 SUBROUTINE sgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
217 $ lwork, iwork, info )
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...
LOGICAL function lsame(CA, CB)
LSAME
REAL function slamch(CMACH)
SLAMCH
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
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.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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 sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
SGESDD
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