216 SUBROUTINE dgesdd( 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 DOUBLE PRECISION A( lda, * ), S( * ), U( ldu, * ),
231 $ vt( ldvt, * ), work( * )
237 DOUBLE PRECISION ZERO, ONE
238 parameter( zero = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
250 DOUBLE PRECISION DUM( 1 )
260 DOUBLE PRECISION DLAMCH, DLANGE
261 EXTERNAL dlamch, dlange, ilaenv, lsame
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.0d0 / 6.0d0 )
316 IF( m.GE.mnthr )
THEN
321 wrkbl = n + n*ilaenv( 1,
'DGEQRF',
' ', m, n, -1,
323 wrkbl = max( wrkbl, 3*n+2*n*
324 $ ilaenv( 1,
'DGEBRD',
' ', n, n, -1, -1 ) )
325 maxwrk = max( wrkbl, bdspac+n )
327 ELSE IF( wntqo )
THEN
331 wrkbl = n + n*ilaenv( 1,
'DGEQRF',
' ', m, n, -1, -1 )
332 wrkbl = max( wrkbl, n+n*ilaenv( 1,
'DORGQR',
' ', m,
334 wrkbl = max( wrkbl, 3*n+2*n*
335 $ ilaenv( 1,
'DGEBRD',
' ', n, n, -1, -1 ) )
336 wrkbl = max( wrkbl, 3*n+n*
337 $ ilaenv( 1,
'DORMBR',
'QLN', n, n, n, -1 ) )
338 wrkbl = max( wrkbl, 3*n+n*
339 $ ilaenv( 1,
'DORMBR',
'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,
'DGEQRF',
' ', m, n, -1, -1 )
348 wrkbl = max( wrkbl, n+n*ilaenv( 1,
'DORGQR',
' ', m,
350 wrkbl = max( wrkbl, 3*n+2*n*
351 $ ilaenv( 1,
'DGEBRD',
' ', n, n, -1, -1 ) )
352 wrkbl = max( wrkbl, 3*n+n*
353 $ ilaenv( 1,
'DORMBR',
'QLN', n, n, n, -1 ) )
354 wrkbl = max( wrkbl, 3*n+n*
355 $ ilaenv( 1,
'DORMBR',
'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,
'DGEQRF',
' ', m, n, -1, -1 )
364 wrkbl = max( wrkbl, n+m*ilaenv( 1,
'DORGQR',
' ', m,
366 wrkbl = max( wrkbl, 3*n+2*n*
367 $ ilaenv( 1,
'DGEBRD',
' ', n, n, -1, -1 ) )
368 wrkbl = max( wrkbl, 3*n+n*
369 $ ilaenv( 1,
'DORMBR',
'QLN', n, n, n, -1 ) )
370 wrkbl = max( wrkbl, 3*n+n*
371 $ ilaenv( 1,
'DORMBR',
'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,
'DGEBRD',
' ', 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,
'DORMBR',
'QLN', m, n, n, -1 ) )
388 wrkbl = max( wrkbl, 3*n+n*
389 $ ilaenv( 1,
'DORMBR',
'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,
'DORMBR',
'QLN', m, n, n, -1 ) )
396 wrkbl = max( wrkbl, 3*n+n*
397 $ ilaenv( 1,
'DORMBR',
'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,
'DORMBR',
'QLN', m, m, n, -1 ) )
403 wrkbl = max( wrkbl, 3*n+n*
404 $ ilaenv( 1,
'DORMBR',
'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.0d0 / 6.0d0 )
419 IF( n.GE.mnthr )
THEN
424 wrkbl = m + m*ilaenv( 1,
'DGELQF',
' ', m, n, -1,
426 wrkbl = max( wrkbl, 3*m+2*m*
427 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
428 maxwrk = max( wrkbl, bdspac+m )
430 ELSE IF( wntqo )
THEN
434 wrkbl = m + m*ilaenv( 1,
'DGELQF',
' ', m, n, -1, -1 )
435 wrkbl = max( wrkbl, m+m*ilaenv( 1,
'DORGLQ',
' ', m,
437 wrkbl = max( wrkbl, 3*m+2*m*
438 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
439 wrkbl = max( wrkbl, 3*m+m*
440 $ ilaenv( 1,
'DORMBR',
'QLN', m, m, m, -1 ) )
441 wrkbl = max( wrkbl, 3*m+m*
442 $ ilaenv( 1,
'DORMBR',
'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,
'DGELQF',
' ', m, n, -1, -1 )
451 wrkbl = max( wrkbl, m+m*ilaenv( 1,
'DORGLQ',
' ', m,
453 wrkbl = max( wrkbl, 3*m+2*m*
454 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
455 wrkbl = max( wrkbl, 3*m+m*
456 $ ilaenv( 1,
'DORMBR',
'QLN', m, m, m, -1 ) )
457 wrkbl = max( wrkbl, 3*m+m*
458 $ ilaenv( 1,
'DORMBR',
'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,
'DGELQF',
' ', m, n, -1, -1 )
467 wrkbl = max( wrkbl, m+n*ilaenv( 1,
'DORGLQ',
' ', n,
469 wrkbl = max( wrkbl, 3*m+2*m*
470 $ ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
471 wrkbl = max( wrkbl, 3*m+m*
472 $ ilaenv( 1,
'DORMBR',
'QLN', m, m, m, -1 ) )
473 wrkbl = max( wrkbl, 3*m+m*
474 $ ilaenv( 1,
'DORMBR',
'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,
'DGEBRD',
' ', 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,
'DORMBR',
'QLN', m, m, n, -1 ) )
491 wrkbl = max( wrkbl, 3*m+m*
492 $ ilaenv( 1,
'DORMBR',
'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,
'DORMBR',
'QLN', m, m, n, -1 ) )
499 wrkbl = max( wrkbl, 3*m+m*
500 $ ilaenv( 1,
'DORMBR',
'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,
'DORMBR',
'QLN', m, m, n, -1 ) )
506 wrkbl = max( wrkbl, 3*m+m*
507 $ ilaenv( 1,
'DORMBR',
'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(
'DGESDD', -info )
524 ELSE IF( lquery )
THEN
530 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
537 smlnum = sqrt( dlamch(
'S' ) ) / eps
538 bignum = one / smlnum
542 anrm = dlange(
'M', m, n, a, lda, dum )
544 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
546 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
547 ELSE IF( anrm.GT.bignum )
THEN
549 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
558 IF( m.GE.mnthr )
THEN
571 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
572 $ lwork-nwork+1, ierr )
576 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
585 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
586 $ work( itaup ), work( nwork ), lwork-nwork+1,
593 CALL dbdsdc(
'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 dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
618 $ lwork-nwork+1, ierr )
622 CALL dlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
623 CALL dlaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
629 CALL dorgqr( m, n, n, a, lda, work( itau ),
630 $ work( nwork ), lwork-nwork+1, ierr )
639 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
640 $ work( itauq ), work( itaup ), work( nwork ),
641 $ lwork-nwork+1, ierr )
653 CALL dbdsdc(
'U',
'I', n, s, work( ie ), work( iu ), n,
654 $ vt, ldvt, dum, idum, work( nwork ), iwork,
661 CALL dormbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
662 $ work( itauq ), work( iu ), n, work( nwork ),
663 $ lwork-nwork+1, ierr )
664 CALL dormbr(
'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 dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
675 $ lda, work( iu ), n, zero, work( ir ),
677 CALL dlacpy(
'F', chunk, n, work( ir ), ldwrkr,
681 ELSE IF( wntqs )
THEN
698 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
699 $ lwork-nwork+1, ierr )
703 CALL dlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
704 CALL dlaset(
'L', n-1, n-1, zero, zero, work( ir+1 ),
710 CALL dorgqr( m, n, n, a, lda, work( itau ),
711 $ work( nwork ), lwork-nwork+1, ierr )
720 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
721 $ work( itauq ), work( itaup ), work( nwork ),
722 $ lwork-nwork+1, ierr )
729 CALL dbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
730 $ ldvt, dum, idum, work( nwork ), iwork,
737 CALL dormbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
738 $ work( itauq ), u, ldu, work( nwork ),
739 $ lwork-nwork+1, ierr )
741 CALL dormbr(
'P',
'R',
'T', n, n, n, work( ir ), ldwrkr,
742 $ work( itaup ), vt, ldvt, work( nwork ),
743 $ lwork-nwork+1, ierr )
749 CALL dlacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
750 CALL dgemm(
'N',
'N', m, n, n, one, a, lda, work( ir ),
751 $ ldwrkr, zero, u, ldu )
753 ELSE IF( wntqa )
THEN
770 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
771 $ lwork-nwork+1, ierr )
772 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
776 CALL dorgqr( m, m, n, u, ldu, work( itau ),
777 $ work( nwork ), lwork-nwork+1, ierr )
781 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
790 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
791 $ work( itaup ), work( nwork ), lwork-nwork+1,
799 CALL dbdsdc(
'U',
'I', n, s, work( ie ), work( iu ), n,
800 $ vt, ldvt, dum, idum, work( nwork ), iwork,
807 CALL dormbr(
'Q',
'L',
'N', n, n, n, a, lda,
808 $ work( itauq ), work( iu ), ldwrku,
809 $ work( nwork ), lwork-nwork+1, ierr )
810 CALL dormbr(
'P',
'R',
'T', n, n, n, a, lda,
811 $ work( itaup ), vt, ldvt, work( nwork ),
812 $ lwork-nwork+1, ierr )
818 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu, work( iu ),
819 $ ldwrku, zero, a, lda )
823 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
842 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
843 $ work( itaup ), work( nwork ), lwork-nwork+1,
850 CALL dbdsdc(
'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 dlaset(
'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 dbdsdc(
'U',
'I', n, s, work( ie ), work( iu ),
882 $ ldwrku, vt, ldvt, dum, idum, work( nwork ),
888 CALL dormbr(
'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 dormbr(
'Q',
'L',
'N', m, n, n, a, lda,
898 $ work( itauq ), work( iu ), ldwrku,
899 $ work( nwork ), lwork-nwork+1, ierr )
903 CALL dlacpy(
'F', m, n, work( iu ), ldwrku, a, lda )
909 CALL dorgbr(
'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 dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
920 $ lda, work( iu ), ldwrku, zero,
921 $ work( ir ), ldwrkr )
922 CALL dlacpy(
'F', chunk, n, work( ir ), ldwrkr,
927 ELSE IF( wntqs )
THEN
934 CALL dlaset(
'F', m, n, zero, zero, u, ldu )
935 CALL dbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
936 $ ldvt, dum, idum, work( nwork ), iwork,
943 CALL dormbr(
'Q',
'L',
'N', m, n, n, a, lda,
944 $ work( itauq ), u, ldu, work( nwork ),
945 $ lwork-nwork+1, ierr )
946 CALL dormbr(
'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 dlaset(
'F', m, m, zero, zero, u, ldu )
957 CALL dbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
958 $ ldvt, dum, idum, work( nwork ), iwork,
964 CALL dlaset(
'F', m-n, m-n, zero, one, u( n+1, n+1 ),
972 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
973 $ work( itauq ), u, ldu, work( nwork ),
974 $ lwork-nwork+1, ierr )
975 CALL dormbr(
'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 dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1002 $ lwork-nwork+1, ierr )
1006 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1015 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1016 $ work( itaup ), work( nwork ), lwork-nwork+1,
1023 CALL dbdsdc(
'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 dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1054 $ lwork-nwork+1, ierr )
1058 CALL dlacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1059 CALL dlaset(
'U', m-1, m-1, zero, zero,
1060 $ work( il+ldwrkl ), ldwrkl )
1065 CALL dorglq( m, n, m, a, lda, work( itau ),
1066 $ work( nwork ), lwork-nwork+1, ierr )
1075 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1076 $ work( itauq ), work( itaup ), work( nwork ),
1077 $ lwork-nwork+1, ierr )
1084 CALL dbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1085 $ work( ivt ), m, dum, idum, work( nwork ),
1092 CALL dormbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1093 $ work( itauq ), u, ldu, work( nwork ),
1094 $ lwork-nwork+1, ierr )
1095 CALL dormbr(
'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 dgemm(
'N',
'N', m, blk, m, one, work( ivt ), m,
1106 $ a( 1, i ), lda, zero, work( il ), ldwrkl )
1107 CALL dlacpy(
'F', m, blk, work( il ), ldwrkl,
1111 ELSE IF( wntqs )
THEN
1122 itau = il + ldwrkl*m
1128 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1129 $ lwork-nwork+1, ierr )
1133 CALL dlacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1134 CALL dlaset(
'U', m-1, m-1, zero, zero,
1135 $ work( il+ldwrkl ), ldwrkl )
1140 CALL dorglq( m, n, m, a, lda, work( itau ),
1141 $ work( nwork ), lwork-nwork+1, ierr )
1150 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1151 $ work( itauq ), work( itaup ), work( nwork ),
1152 $ lwork-nwork+1, ierr )
1159 CALL dbdsdc(
'U',
'I', m, s, work( ie ), u, ldu, vt,
1160 $ ldvt, dum, idum, work( nwork ), iwork,
1167 CALL dormbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1168 $ work( itauq ), u, ldu, work( nwork ),
1169 $ lwork-nwork+1, ierr )
1170 CALL dormbr(
'P',
'R',
'T', m, m, m, work( il ), ldwrkl,
1171 $ work( itaup ), vt, ldvt, work( nwork ),
1172 $ lwork-nwork+1, ierr )
1178 CALL dlacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1179 CALL dgemm(
'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 dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1200 $ lwork-nwork+1, ierr )
1201 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
1206 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
1207 $ work( nwork ), lwork-nwork+1, ierr )
1211 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1220 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1221 $ work( itaup ), work( nwork ), lwork-nwork+1,
1229 CALL dbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1230 $ work( ivt ), ldwkvt, dum, idum,
1231 $ work( nwork ), iwork, info )
1237 CALL dormbr(
'Q',
'L',
'N', m, m, m, a, lda,
1238 $ work( itauq ), u, ldu, work( nwork ),
1239 $ lwork-nwork+1, ierr )
1240 CALL dormbr(
'P',
'R',
'T', m, m, m, a, lda,
1241 $ work( itaup ), work( ivt ), ldwkvt,
1242 $ work( nwork ), lwork-nwork+1, ierr )
1248 CALL dgemm(
'N',
'N', m, n, m, one, work( ivt ), ldwkvt,
1249 $ vt, ldvt, zero, a, lda )
1253 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
1272 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1273 $ work( itaup ), work( nwork ), lwork-nwork+1,
1280 CALL dbdsdc(
'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 dlaset(
'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 dbdsdc(
'L',
'I', m, s, work( ie ), u, ldu,
1310 $ work( ivt ), ldwkvt, dum, idum,
1311 $ work( nwork ), iwork, info )
1316 CALL dormbr(
'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 dormbr(
'P',
'R',
'T', m, n, m, a, lda,
1326 $ work( itaup ), work( ivt ), ldwkvt,
1327 $ work( nwork ), lwork-nwork+1, ierr )
1331 CALL dlacpy(
'F', m, n, work( ivt ), ldwkvt, a, lda )
1337 CALL dorgbr(
'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 dgemm(
'N',
'N', m, blk, m, one, work( ivt ),
1348 $ ldwkvt, a( 1, i ), lda, zero,
1350 CALL dlacpy(
'F', m, blk, work( il ), m, a( 1, i ),
1354 ELSE IF( wntqs )
THEN
1361 CALL dlaset(
'F', m, n, zero, zero, vt, ldvt )
1362 CALL dbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1363 $ ldvt, dum, idum, work( nwork ), iwork,
1370 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1371 $ work( itauq ), u, ldu, work( nwork ),
1372 $ lwork-nwork+1, ierr )
1373 CALL dormbr(
'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 dlaset(
'F', n, n, zero, zero, vt, ldvt )
1384 CALL dbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1385 $ ldvt, dum, idum, work( nwork ), iwork,
1391 CALL dlaset(
'F', n-m, n-m, zero, one, vt( m+1, m+1 ),
1399 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1400 $ work( itauq ), u, ldu, work( nwork ),
1401 $ lwork-nwork+1, ierr )
1402 CALL dormbr(
'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 dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
1417 IF( anrm.LT.smlnum )
1418 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESDD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...