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 )
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.
logical function lsame(CA, CB)
LSAME
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
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
double precision function dlamch(CMACH)
DLAMCH
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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...