209 SUBROUTINE sgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
210 $ WORK, LWORK, IWORK, INFO )
219 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
223 REAL A( LDA, * ), S( * ), U( LDU, * ),
224 $ vt( ldvt, * ), work( * )
231 parameter( zero = 0.0e0, one = 1.0e0 )
234 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
235 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
236 $ ir, iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
237 $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
238 $ mnthr, nwork, wrkbl
239 INTEGER LWORK_SGEBRD_MN, LWORK_SGEBRD_MM,
240 $ lwork_sgebrd_nn, lwork_sgelqf_mn,
242 $ lwork_sorgbr_p_mm, lwork_sorgbr_q_nn,
243 $ lwork_sorglq_mn, lwork_sorglq_nn,
244 $ lwork_sorgqr_mm, lwork_sorgqr_mn,
245 $ lwork_sormbr_prt_mm, lwork_sormbr_qln_mm,
246 $ lwork_sormbr_prt_mn, lwork_sormbr_qln_mn,
247 $ lwork_sormbr_prt_nn, lwork_sormbr_qln_nn
248 REAL ANRM, BIGNUM, EPS, SMLNUM
261 LOGICAL LSAME, SISNAN
262 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
263 EXTERNAL slamch, slange, lsame, sisnan,
267 INTRINSIC int, max, min, sqrt
275 wntqa = lsame( jobz,
'A' )
276 wntqs = lsame( jobz,
'S' )
277 wntqas = wntqa .OR. wntqs
278 wntqo = lsame( jobz,
'O' )
279 wntqn = lsame( jobz,
'N' )
280 lquery = ( lwork.EQ.-1 )
282 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) )
THEN
284 ELSE IF( m.LT.0 )
THEN
286 ELSE IF( n.LT.0 )
THEN
288 ELSE IF( lda.LT.max( 1, m ) )
THEN
290 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
291 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) )
THEN
293 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
294 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
295 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) )
THEN
310 mnthr = int( real( minmn )*11.0e0 / 6.0e0 )
311 IF( m.GE.n .AND. minmn.GT.0 )
THEN
324 CALL sgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
325 $ dum(1), dum(1), -1, ierr )
326 lwork_sgebrd_mn = int( dum(1) )
328 CALL sgebrd( n, n, dum(1), n, dum(1), dum(1), dum(1),
329 $ dum(1), dum(1), -1, ierr )
330 lwork_sgebrd_nn = int( dum(1) )
332 CALL sgeqrf( m, n, dum(1), m, dum(1), dum(1), -1, ierr )
333 lwork_sgeqrf_mn = int( dum(1) )
335 CALL sorgbr(
'Q', n, n, n, dum(1), n, dum(1), dum(1), -1,
337 lwork_sorgbr_q_nn = int( dum(1) )
339 CALL sorgqr( m, m, n, dum(1), m, dum(1), dum(1), -1,
341 lwork_sorgqr_mm = int( dum(1) )
343 CALL sorgqr( m, n, n, dum(1), m, dum(1), dum(1), -1,
345 lwork_sorgqr_mn = int( dum(1) )
347 CALL sormbr(
'P',
'R',
'T', n, n, n, dum(1), n,
348 $ dum(1), dum(1), n, dum(1), -1, ierr )
349 lwork_sormbr_prt_nn = int( dum(1) )
351 CALL sormbr(
'Q',
'L',
'N', n, n, n, dum(1), n,
352 $ dum(1), dum(1), n, dum(1), -1, ierr )
353 lwork_sormbr_qln_nn = int( dum(1) )
355 CALL sormbr(
'Q',
'L',
'N', m, n, n, dum(1), m,
356 $ dum(1), dum(1), m, dum(1), -1, ierr )
357 lwork_sormbr_qln_mn = int( dum(1) )
359 CALL sormbr(
'Q',
'L',
'N', m, m, n, dum(1), m,
360 $ dum(1), dum(1), m, dum(1), -1, ierr )
361 lwork_sormbr_qln_mm = int( dum(1) )
363 IF( m.GE.mnthr )
THEN
368 wrkbl = n + lwork_sgeqrf_mn
369 wrkbl = max( wrkbl, 3*n + lwork_sgebrd_nn )
370 maxwrk = max( wrkbl, bdspac + n )
372 ELSE IF( wntqo )
THEN
376 wrkbl = n + lwork_sgeqrf_mn
377 wrkbl = max( wrkbl, n + lwork_sorgqr_mn )
378 wrkbl = max( wrkbl, 3*n + lwork_sgebrd_nn )
379 wrkbl = max( wrkbl, 3*n + lwork_sormbr_qln_nn )
380 wrkbl = max( wrkbl, 3*n + lwork_sormbr_prt_nn )
381 wrkbl = max( wrkbl, 3*n + bdspac )
382 maxwrk = wrkbl + 2*n*n
383 minwrk = bdspac + 2*n*n + 3*n
384 ELSE IF( wntqs )
THEN
388 wrkbl = n + lwork_sgeqrf_mn
389 wrkbl = max( wrkbl, n + lwork_sorgqr_mn )
390 wrkbl = max( wrkbl, 3*n + lwork_sgebrd_nn )
391 wrkbl = max( wrkbl, 3*n + lwork_sormbr_qln_nn )
392 wrkbl = max( wrkbl, 3*n + lwork_sormbr_prt_nn )
393 wrkbl = max( wrkbl, 3*n + bdspac )
395 minwrk = bdspac + n*n + 3*n
396 ELSE IF( wntqa )
THEN
400 wrkbl = n + lwork_sgeqrf_mn
401 wrkbl = max( wrkbl, n + lwork_sorgqr_mm )
402 wrkbl = max( wrkbl, 3*n + lwork_sgebrd_nn )
403 wrkbl = max( wrkbl, 3*n + lwork_sormbr_qln_nn )
404 wrkbl = max( wrkbl, 3*n + lwork_sormbr_prt_nn )
405 wrkbl = max( wrkbl, 3*n + bdspac )
407 minwrk = n*n + max( 3*n + bdspac, n + m )
413 wrkbl = 3*n + lwork_sgebrd_mn
416 maxwrk = max( wrkbl, 3*n + bdspac )
417 minwrk = 3*n + max( m, bdspac )
418 ELSE IF( wntqo )
THEN
420 wrkbl = max( wrkbl, 3*n + lwork_sormbr_prt_nn )
421 wrkbl = max( wrkbl, 3*n + lwork_sormbr_qln_mn )
422 wrkbl = max( wrkbl, 3*n + bdspac )
424 minwrk = 3*n + max( m, n*n + bdspac )
425 ELSE IF( wntqs )
THEN
427 wrkbl = max( wrkbl, 3*n + lwork_sormbr_qln_mn )
428 wrkbl = max( wrkbl, 3*n + lwork_sormbr_prt_nn )
429 maxwrk = max( wrkbl, 3*n + bdspac )
430 minwrk = 3*n + max( m, bdspac )
431 ELSE IF( wntqa )
THEN
433 wrkbl = max( wrkbl, 3*n + lwork_sormbr_qln_mm )
434 wrkbl = max( wrkbl, 3*n + lwork_sormbr_prt_nn )
435 maxwrk = max( wrkbl, 3*n + bdspac )
436 minwrk = 3*n + max( m, bdspac )
439 ELSE IF( minmn.GT.0 )
THEN
452 CALL sgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
453 $ dum(1), dum(1), -1, ierr )
454 lwork_sgebrd_mn = int( dum(1) )
456 CALL sgebrd( m, m, a, m, s, dum(1), dum(1),
457 $ dum(1), dum(1), -1, ierr )
458 lwork_sgebrd_mm = int( dum(1) )
460 CALL sgelqf( m, n, a, m, dum(1), dum(1), -1, ierr )
461 lwork_sgelqf_mn = int( dum(1) )
463 CALL sorglq( n, n, m, dum(1), n, dum(1), dum(1), -1,
465 lwork_sorglq_nn = int( dum(1) )
467 CALL sorglq( m, n, m, a, m, dum(1), dum(1), -1, ierr )
468 lwork_sorglq_mn = int( dum(1) )
470 CALL sorgbr(
'P', m, m, m, a, n, dum(1), dum(1), -1,
472 lwork_sorgbr_p_mm = int( dum(1) )
474 CALL sormbr(
'P',
'R',
'T', m, m, m, dum(1), m,
475 $ dum(1), dum(1), m, dum(1), -1, ierr )
476 lwork_sormbr_prt_mm = int( dum(1) )
478 CALL sormbr(
'P',
'R',
'T', m, n, m, dum(1), m,
479 $ dum(1), dum(1), m, dum(1), -1, ierr )
480 lwork_sormbr_prt_mn = int( dum(1) )
482 CALL sormbr(
'P',
'R',
'T', n, n, m, dum(1), n,
483 $ dum(1), dum(1), n, dum(1), -1, ierr )
484 lwork_sormbr_prt_nn = int( dum(1) )
486 CALL sormbr(
'Q',
'L',
'N', m, m, m, dum(1), m,
487 $ dum(1), dum(1), m, dum(1), -1, ierr )
488 lwork_sormbr_qln_mm = int( dum(1) )
490 IF( n.GE.mnthr )
THEN
495 wrkbl = m + lwork_sgelqf_mn
496 wrkbl = max( wrkbl, 3*m + lwork_sgebrd_mm )
497 maxwrk = max( wrkbl, bdspac + m )
499 ELSE IF( wntqo )
THEN
503 wrkbl = m + lwork_sgelqf_mn
504 wrkbl = max( wrkbl, m + lwork_sorglq_mn )
505 wrkbl = max( wrkbl, 3*m + lwork_sgebrd_mm )
506 wrkbl = max( wrkbl, 3*m + lwork_sormbr_qln_mm )
507 wrkbl = max( wrkbl, 3*m + lwork_sormbr_prt_mm )
508 wrkbl = max( wrkbl, 3*m + bdspac )
509 maxwrk = wrkbl + 2*m*m
510 minwrk = bdspac + 2*m*m + 3*m
511 ELSE IF( wntqs )
THEN
515 wrkbl = m + lwork_sgelqf_mn
516 wrkbl = max( wrkbl, m + lwork_sorglq_mn )
517 wrkbl = max( wrkbl, 3*m + lwork_sgebrd_mm )
518 wrkbl = max( wrkbl, 3*m + lwork_sormbr_qln_mm )
519 wrkbl = max( wrkbl, 3*m + lwork_sormbr_prt_mm )
520 wrkbl = max( wrkbl, 3*m + bdspac )
522 minwrk = bdspac + m*m + 3*m
523 ELSE IF( wntqa )
THEN
527 wrkbl = m + lwork_sgelqf_mn
528 wrkbl = max( wrkbl, m + lwork_sorglq_nn )
529 wrkbl = max( wrkbl, 3*m + lwork_sgebrd_mm )
530 wrkbl = max( wrkbl, 3*m + lwork_sormbr_qln_mm )
531 wrkbl = max( wrkbl, 3*m + lwork_sormbr_prt_mm )
532 wrkbl = max( wrkbl, 3*m + bdspac )
534 minwrk = m*m + max( 3*m + bdspac, m + n )
540 wrkbl = 3*m + lwork_sgebrd_mn
543 maxwrk = max( wrkbl, 3*m + bdspac )
544 minwrk = 3*m + max( n, bdspac )
545 ELSE IF( wntqo )
THEN
547 wrkbl = max( wrkbl, 3*m + lwork_sormbr_qln_mm )
548 wrkbl = max( wrkbl, 3*m + lwork_sormbr_prt_mn )
549 wrkbl = max( wrkbl, 3*m + bdspac )
551 minwrk = 3*m + max( n, m*m + bdspac )
552 ELSE IF( wntqs )
THEN
554 wrkbl = max( wrkbl, 3*m + lwork_sormbr_qln_mm )
555 wrkbl = max( wrkbl, 3*m + lwork_sormbr_prt_mn )
556 maxwrk = max( wrkbl, 3*m + bdspac )
557 minwrk = 3*m + max( n, bdspac )
558 ELSE IF( wntqa )
THEN
560 wrkbl = max( wrkbl, 3*m + lwork_sormbr_qln_mm )
561 wrkbl = max( wrkbl, 3*m + lwork_sormbr_prt_nn )
562 maxwrk = max( wrkbl, 3*m + bdspac )
563 minwrk = 3*m + max( n, bdspac )
568 maxwrk = max( maxwrk, minwrk )
569 work( 1 ) = sroundup_lwork( maxwrk )
571 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
577 CALL xerbla(
'SGESDD', -info )
579 ELSE IF( lquery )
THEN
585 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
592 smlnum = sqrt( slamch(
'S' ) ) / eps
593 bignum = one / smlnum
597 anrm = slange(
'M', m, n, a, lda, dum )
598 IF( sisnan( anrm ) )
THEN
603 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
605 CALL slascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
606 ELSE IF( anrm.GT.bignum )
THEN
608 CALL slascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
617 IF( m.GE.mnthr )
THEN
631 CALL sgeqrf( m, n, a, lda, work( itau ),
633 $ lwork - nwork + 1, ierr )
637 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
648 CALL sgebrd( n, n, a, lda, s, work( ie ),
650 $ work( itaup ), work( nwork ), lwork-nwork+1,
657 CALL sbdsdc(
'U',
'N', n, s, work( ie ), dum, 1, dum,
659 $ dum, idum, work( nwork ), iwork, info )
661 ELSE IF( wntqo )
THEN
671 IF( lwork .GE. lda*n + n*n + 3*n + bdspac )
THEN
674 ldwrkr = ( lwork - n*n - 3*n - bdspac ) / n
683 CALL sgeqrf( m, n, a, lda, work( itau ),
685 $ lwork - nwork + 1, ierr )
689 CALL slacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
690 CALL slaset(
'L', n - 1, n - 1, zero, zero,
698 CALL sorgqr( m, n, n, a, lda, work( itau ),
699 $ work( nwork ), lwork - nwork + 1, ierr )
709 CALL sgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
710 $ work( itauq ), work( itaup ), work( nwork ),
711 $ lwork - nwork + 1, ierr )
723 CALL sbdsdc(
'U',
'I', n, s, work( ie ), work( iu ),
725 $ vt, ldvt, dum, idum, work( nwork ), iwork,
733 CALL sormbr(
'Q',
'L',
'N', n, n, n, work( ir ),
735 $ work( itauq ), work( iu ), n, work( nwork ),
736 $ lwork - nwork + 1, ierr )
737 CALL sormbr(
'P',
'R',
'T', n, n, n, work( ir ),
739 $ work( itaup ), vt, ldvt, work( nwork ),
740 $ lwork - nwork + 1, ierr )
747 DO 10 i = 1, m, ldwrkr
748 chunk = min( m - i + 1, ldwrkr )
749 CALL sgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
750 $ lda, work( iu ), n, zero, work( ir ),
752 CALL slacpy(
'F', chunk, n, work( ir ), ldwrkr,
756 ELSE IF( wntqs )
THEN
774 CALL sgeqrf( m, n, a, lda, work( itau ),
776 $ lwork - nwork + 1, ierr )
780 CALL slacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
781 CALL slaset(
'L', n - 1, n - 1, zero, zero,
789 CALL sorgqr( m, n, n, a, lda, work( itau ),
790 $ work( nwork ), lwork - nwork + 1, ierr )
800 CALL sgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
801 $ work( itauq ), work( itaup ), work( nwork ),
802 $ lwork - nwork + 1, ierr )
809 CALL sbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
810 $ ldvt, dum, idum, work( nwork ), iwork,
818 CALL sormbr(
'Q',
'L',
'N', n, n, n, work( ir ),
820 $ work( itauq ), u, ldu, work( nwork ),
821 $ lwork - nwork + 1, ierr )
823 CALL sormbr(
'P',
'R',
'T', n, n, n, work( ir ),
825 $ work( itaup ), vt, ldvt, work( nwork ),
826 $ lwork - nwork + 1, ierr )
832 CALL slacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
833 CALL sgemm(
'N',
'N', m, n, n, one, a, lda,
835 $ ldwrkr, zero, u, ldu )
837 ELSE IF( wntqa )
THEN
855 CALL sgeqrf( m, n, a, lda, work( itau ),
857 $ lwork - nwork + 1, ierr )
858 CALL slacpy(
'L', m, n, a, lda, u, ldu )
863 CALL sorgqr( m, m, n, u, ldu, work( itau ),
864 $ work( nwork ), lwork - nwork + 1, ierr )
868 CALL slaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ),
879 CALL sgebrd( n, n, a, lda, s, work( ie ),
881 $ work( itaup ), work( nwork ), lwork-nwork+1,
889 CALL sbdsdc(
'U',
'I', n, s, work( ie ), work( iu ),
891 $ vt, ldvt, dum, idum, work( nwork ), iwork,
899 CALL sormbr(
'Q',
'L',
'N', n, n, n, a, lda,
900 $ work( itauq ), work( iu ), ldwrku,
901 $ work( nwork ), lwork - nwork + 1, ierr )
902 CALL sormbr(
'P',
'R',
'T', n, n, n, a, lda,
903 $ work( itaup ), vt, ldvt, work( nwork ),
904 $ lwork - nwork + 1, ierr )
910 CALL sgemm(
'N',
'N', m, n, n, one, u, ldu,
912 $ ldwrku, zero, a, lda )
916 CALL slacpy(
'F', m, n, a, lda, u, ldu )
936 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
937 $ work( itaup ), work( nwork ), lwork-nwork+1,
945 CALL sbdsdc(
'U',
'N', n, s, work( ie ), dum, 1, dum,
947 $ dum, idum, work( nwork ), iwork, info )
948 ELSE IF( wntqo )
THEN
951 IF( lwork .GE. m*n + 3*n + bdspac )
THEN
956 nwork = iu + ldwrku*n
957 CALL slaset(
'F', m, n, zero, zero, work( iu ),
966 nwork = iu + ldwrku*n
971 ldwrkr = ( lwork - n*n - 3*n ) / n
973 nwork = iu + ldwrku*n
980 CALL sbdsdc(
'U',
'I', n, s, work( ie ), work( iu ),
981 $ ldwrku, vt, ldvt, dum, idum, work( nwork ),
988 CALL sormbr(
'P',
'R',
'T', n, n, n, a, lda,
989 $ work( itaup ), vt, ldvt, work( nwork ),
990 $ lwork - nwork + 1, ierr )
992 IF( lwork .GE. m*n + 3*n + bdspac )
THEN
999 CALL sormbr(
'Q',
'L',
'N', m, n, n, a, lda,
1000 $ work( itauq ), work( iu ), ldwrku,
1001 $ work( nwork ), lwork - nwork + 1, ierr )
1005 CALL slacpy(
'F', m, n, work( iu ), ldwrku, a,
1014 CALL sorgbr(
'Q', m, n, n, a, lda, work( itauq ),
1015 $ work( nwork ), lwork - nwork + 1, ierr )
1023 DO 20 i = 1, m, ldwrkr
1024 chunk = min( m - i + 1, ldwrkr )
1025 CALL sgemm(
'N',
'N', chunk, n, n, one, a( i,
1027 $ lda, work( iu ), ldwrku, zero,
1028 $ work( ir ), ldwrkr )
1029 CALL slacpy(
'F', chunk, n, work( ir ), ldwrkr,
1034 ELSE IF( wntqs )
THEN
1042 CALL slaset(
'F', m, n, zero, zero, u, ldu )
1043 CALL sbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
1044 $ ldvt, dum, idum, work( nwork ), iwork,
1052 CALL sormbr(
'Q',
'L',
'N', m, n, n, a, lda,
1053 $ work( itauq ), u, ldu, work( nwork ),
1054 $ lwork - nwork + 1, ierr )
1055 CALL sormbr(
'P',
'R',
'T', n, n, n, a, lda,
1056 $ work( itaup ), vt, ldvt, work( nwork ),
1057 $ lwork - nwork + 1, ierr )
1058 ELSE IF( wntqa )
THEN
1066 CALL slaset(
'F', m, m, zero, zero, u, ldu )
1067 CALL sbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
1068 $ ldvt, dum, idum, work( nwork ), iwork,
1074 CALL slaset(
'F', m - n, m - n, zero, one, u(n+1,
1084 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1085 $ work( itauq ), u, ldu, work( nwork ),
1086 $ lwork - nwork + 1, ierr )
1087 CALL sormbr(
'P',
'R',
'T', n, n, m, a, lda,
1088 $ work( itaup ), vt, ldvt, work( nwork ),
1089 $ lwork - nwork + 1, ierr )
1100 IF( n.GE.mnthr )
THEN
1114 CALL sgelqf( m, n, a, lda, work( itau ),
1116 $ lwork - nwork + 1, ierr )
1120 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
1131 CALL sgebrd( m, m, a, lda, s, work( ie ),
1133 $ work( itaup ), work( nwork ), lwork-nwork+1,
1140 CALL sbdsdc(
'U',
'N', m, s, work( ie ), dum, 1, dum,
1142 $ dum, idum, work( nwork ), iwork, info )
1144 ELSE IF( wntqo )
THEN
1156 IF( lwork .GE. m*n + m*m + 3*m + bdspac )
THEN
1161 chunk = ( lwork - m*m ) / m
1163 itau = il + ldwrkl*m
1170 CALL sgelqf( m, n, a, lda, work( itau ),
1172 $ lwork - nwork + 1, ierr )
1176 CALL slacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1177 CALL slaset(
'U', m - 1, m - 1, zero, zero,
1178 $ work( il + ldwrkl ), ldwrkl )
1184 CALL sorglq( m, n, m, a, lda, work( itau ),
1185 $ work( nwork ), lwork - nwork + 1, ierr )
1195 CALL sgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1196 $ work( itauq ), work( itaup ), work( nwork ),
1197 $ lwork - nwork + 1, ierr )
1204 CALL sbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1205 $ work( ivt ), m, dum, idum, work( nwork ),
1213 CALL sormbr(
'Q',
'L',
'N', m, m, m, work( il ),
1215 $ work( itauq ), u, ldu, work( nwork ),
1216 $ lwork - nwork + 1, ierr )
1217 CALL sormbr(
'P',
'R',
'T', m, m, m, work( il ),
1219 $ work( itaup ), work( ivt ), m,
1220 $ work( nwork ), lwork - nwork + 1, ierr )
1228 DO 30 i = 1, n, chunk
1229 blk = min( n - i + 1, chunk )
1230 CALL sgemm(
'N',
'N', m, blk, m, one, work( ivt ),
1232 $ a( 1, i ), lda, zero, work( il ), ldwrkl )
1233 CALL slacpy(
'F', m, blk, work( il ), ldwrkl,
1237 ELSE IF( wntqs )
THEN
1248 itau = il + ldwrkl*m
1255 CALL sgelqf( m, n, a, lda, work( itau ),
1257 $ lwork - nwork + 1, ierr )
1261 CALL slacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1262 CALL slaset(
'U', m - 1, m - 1, zero, zero,
1263 $ work( il + ldwrkl ), ldwrkl )
1269 CALL sorglq( m, n, m, a, lda, work( itau ),
1270 $ work( nwork ), lwork - nwork + 1, ierr )
1280 CALL sgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1281 $ work( itauq ), work( itaup ), work( nwork ),
1282 $ lwork - nwork + 1, ierr )
1289 CALL sbdsdc(
'U',
'I', m, s, work( ie ), u, ldu, vt,
1290 $ ldvt, dum, idum, work( nwork ), iwork,
1298 CALL sormbr(
'Q',
'L',
'N', m, m, m, work( il ),
1300 $ work( itauq ), u, ldu, work( nwork ),
1301 $ lwork - nwork + 1, ierr )
1302 CALL sormbr(
'P',
'R',
'T', m, m, m, work( il ),
1304 $ work( itaup ), vt, ldvt, work( nwork ),
1305 $ lwork - nwork + 1, ierr )
1311 CALL slacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1312 CALL sgemm(
'N',
'N', m, n, m, one, work( il ),
1314 $ a, lda, zero, vt, ldvt )
1316 ELSE IF( wntqa )
THEN
1327 itau = ivt + ldwkvt*m
1334 CALL sgelqf( m, n, a, lda, work( itau ),
1336 $ lwork - nwork + 1, ierr )
1337 CALL slacpy(
'U', m, n, a, lda, vt, ldvt )
1343 CALL sorglq( n, n, m, vt, ldvt, work( itau ),
1344 $ work( nwork ), lwork - nwork + 1, ierr )
1348 CALL slaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ),
1359 CALL sgebrd( m, m, a, lda, s, work( ie ),
1361 $ work( itaup ), work( nwork ), lwork-nwork+1,
1369 CALL sbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1370 $ work( ivt ), ldwkvt, dum, idum,
1371 $ work( nwork ), iwork, info )
1378 CALL sormbr(
'Q',
'L',
'N', m, m, m, a, lda,
1379 $ work( itauq ), u, ldu, work( nwork ),
1380 $ lwork - nwork + 1, ierr )
1381 CALL sormbr(
'P',
'R',
'T', m, m, m, a, lda,
1382 $ work( itaup ), work( ivt ), ldwkvt,
1383 $ work( nwork ), lwork - nwork + 1, ierr )
1389 CALL sgemm(
'N',
'N', m, n, m, one, work( ivt ),
1391 $ vt, ldvt, zero, a, lda )
1395 CALL slacpy(
'F', m, n, a, lda, vt, ldvt )
1415 CALL sgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1416 $ work( itaup ), work( nwork ), lwork-nwork+1,
1424 CALL sbdsdc(
'L',
'N', m, s, work( ie ), dum, 1, dum,
1426 $ dum, idum, work( nwork ), iwork, info )
1427 ELSE IF( wntqo )
THEN
1431 IF( lwork .GE. m*n + 3*m + bdspac )
THEN
1435 CALL slaset(
'F', m, n, zero, zero, work( ivt ),
1437 nwork = ivt + ldwkvt*n
1444 nwork = ivt + ldwkvt*m
1449 chunk = ( lwork - m*m - 3*m ) / m
1457 CALL sbdsdc(
'L',
'I', m, s, work( ie ), u, ldu,
1458 $ work( ivt ), ldwkvt, dum, idum,
1459 $ work( nwork ), iwork, info )
1465 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1466 $ work( itauq ), u, ldu, work( nwork ),
1467 $ lwork - nwork + 1, ierr )
1469 IF( lwork .GE. m*n + 3*m + bdspac )
THEN
1476 CALL sormbr(
'P',
'R',
'T', m, n, m, a, lda,
1477 $ work( itaup ), work( ivt ), ldwkvt,
1478 $ work( nwork ), lwork - nwork + 1, ierr )
1482 CALL slacpy(
'F', m, n, work( ivt ), ldwkvt, a,
1491 CALL sorgbr(
'P', m, n, m, a, lda, work( itaup ),
1492 $ work( nwork ), lwork - nwork + 1, ierr )
1500 DO 40 i = 1, n, chunk
1501 blk = min( n - i + 1, chunk )
1502 CALL sgemm(
'N',
'N', m, blk, m, one,
1504 $ ldwkvt, a( 1, i ), lda, zero,
1506 CALL slacpy(
'F', m, blk, work( il ), m, a( 1,
1511 ELSE IF( wntqs )
THEN
1519 CALL slaset(
'F', m, n, zero, zero, vt, ldvt )
1520 CALL sbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1521 $ ldvt, dum, idum, work( nwork ), iwork,
1529 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1530 $ work( itauq ), u, ldu, work( nwork ),
1531 $ lwork - nwork + 1, ierr )
1532 CALL sormbr(
'P',
'R',
'T', m, n, m, a, lda,
1533 $ work( itaup ), vt, ldvt, work( nwork ),
1534 $ lwork - nwork + 1, ierr )
1535 ELSE IF( wntqa )
THEN
1543 CALL slaset(
'F', n, n, zero, zero, vt, ldvt )
1544 CALL sbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1545 $ ldvt, dum, idum, work( nwork ), iwork,
1551 CALL slaset(
'F', n-m, n-m, zero, one, vt(m+1,m+1),
1560 CALL sormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1561 $ work( itauq ), u, ldu, work( nwork ),
1562 $ lwork - nwork + 1, ierr )
1563 CALL sormbr(
'P',
'R',
'T', n, n, m, a, lda,
1564 $ work( itaup ), vt, ldvt, work( nwork ),
1565 $ lwork - nwork + 1, ierr )
1574 IF( iscl.EQ.1 )
THEN
1575 IF( anrm.GT.bignum )
1576 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
1578 IF( anrm.LT.smlnum )
1579 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
1585 work( 1 ) = sroundup_lwork( maxwrk )