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,