218 SUBROUTINE dgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
219 $ work, lwork, iwork, info )
229 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
233 DOUBLE PRECISION A( lda, * ), S( * ), U( ldu, * ),
234 $ vt( ldvt, * ), work( * )
240 DOUBLE PRECISION ZERO, ONE
241 parameter ( zero = 0.0d0, one = 1.0d0 )
244 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
245 INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
246 $ ir, iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
247 $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
248 $ mnthr, nwork, wrkbl
249 INTEGER LWORK_DGEBRD_MN, LWORK_DGEBRD_MM,
250 $ lwork_dgebrd_nn, lwork_dgelqf_mn,
252 $ lwork_dorgbr_p_mm, lwork_dorgbr_q_nn,
253 $ lwork_dorglq_mn, lwork_dorglq_nn,
254 $ lwork_dorgqr_mm, lwork_dorgqr_mn,
255 $ lwork_dormbr_prt_mm, lwork_dormbr_qln_mm,
256 $ lwork_dormbr_prt_mn, lwork_dormbr_qln_mn,
257 $ lwork_dormbr_prt_nn, lwork_dormbr_qln_nn
258 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
262 DOUBLE PRECISION DUM( 1 )
271 DOUBLE PRECISION DLAMCH, DLANGE
272 EXTERNAL dlamch, dlange, lsame
275 INTRINSIC int, max, min, sqrt
283 wntqa = lsame( jobz,
'A' )
284 wntqs = lsame( jobz,
'S' )
285 wntqas = wntqa .OR. wntqs
286 wntqo = lsame( jobz,
'O' )
287 wntqn = lsame( jobz,
'N' )
288 lquery = ( lwork.EQ.-1 )
290 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) )
THEN
292 ELSE IF( m.LT.0 )
THEN
294 ELSE IF( n.LT.0 )
THEN
296 ELSE IF( lda.LT.max( 1, m ) )
THEN
298 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
299 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) )
THEN
301 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
302 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
303 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) )
THEN
318 mnthr = int( minmn*11.0d0 / 6.0d0 )
319 IF( m.GE.n .AND. minmn.GT.0 )
THEN
332 CALL dgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
333 $ dum(1), dum(1), -1, ierr )
334 lwork_dgebrd_mn = int( dum(1) )
336 CALL dgebrd( n, n, dum(1), n, dum(1), dum(1), dum(1),
337 $ dum(1), dum(1), -1, ierr )
338 lwork_dgebrd_nn = int( dum(1) )
340 CALL dgeqrf( m, n, dum(1), m, dum(1), dum(1), -1, ierr )
341 lwork_dgeqrf_mn = int( dum(1) )
343 CALL dorgbr(
'Q', n, n, n, dum(1), n, dum(1), dum(1), -1,
345 lwork_dorgbr_q_nn = int( dum(1) )
347 CALL dorgqr( m, m, n, dum(1), m, dum(1), dum(1), -1, ierr )
348 lwork_dorgqr_mm = int( dum(1) )
350 CALL dorgqr( m, n, n, dum(1), m, dum(1), dum(1), -1, ierr )
351 lwork_dorgqr_mn = int( dum(1) )
353 CALL dormbr(
'P',
'R',
'T', n, n, n, dum(1), n,
354 $ dum(1), dum(1), n, dum(1), -1, ierr )
355 lwork_dormbr_prt_nn = int( dum(1) )
357 CALL dormbr(
'Q',
'L',
'N', n, n, n, dum(1), n,
358 $ dum(1), dum(1), n, dum(1), -1, ierr )
359 lwork_dormbr_qln_nn = int( dum(1) )
361 CALL dormbr(
'Q',
'L',
'N', m, n, n, dum(1), m,
362 $ dum(1), dum(1), m, dum(1), -1, ierr )
363 lwork_dormbr_qln_mn = int( dum(1) )
365 CALL dormbr(
'Q',
'L',
'N', m, m, n, dum(1), m,
366 $ dum(1), dum(1), m, dum(1), -1, ierr )
367 lwork_dormbr_qln_mm = int( dum(1) )
369 IF( m.GE.mnthr )
THEN
374 wrkbl = n + lwork_dgeqrf_mn
375 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
376 maxwrk = max( wrkbl, bdspac + n )
378 ELSE IF( wntqo )
THEN
382 wrkbl = n + lwork_dgeqrf_mn
383 wrkbl = max( wrkbl, n + lwork_dorgqr_mn )
384 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
385 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
386 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
387 wrkbl = max( wrkbl, 3*n + bdspac )
388 maxwrk = wrkbl + 2*n*n
389 minwrk = bdspac + 2*n*n + 3*n
390 ELSE IF( wntqs )
THEN
394 wrkbl = n + lwork_dgeqrf_mn
395 wrkbl = max( wrkbl, n + lwork_dorgqr_mn )
396 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
397 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
398 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
399 wrkbl = max( wrkbl, 3*n + bdspac )
401 minwrk = bdspac + n*n + 3*n
402 ELSE IF( wntqa )
THEN
406 wrkbl = n + lwork_dgeqrf_mn
407 wrkbl = max( wrkbl, n + lwork_dorgqr_mm )
408 wrkbl = max( wrkbl, 3*n + lwork_dgebrd_nn )
409 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_nn )
410 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
411 wrkbl = max( wrkbl, 3*n + bdspac )
413 minwrk = n*n + max( 3*n + bdspac, n + m )
419 wrkbl = 3*n + lwork_dgebrd_mn
422 maxwrk = max( wrkbl, 3*n + bdspac )
423 minwrk = 3*n + max( m, bdspac )
424 ELSE IF( wntqo )
THEN
426 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
427 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mn )
428 wrkbl = max( wrkbl, 3*n + bdspac )
430 minwrk = 3*n + max( m, n*n + bdspac )
431 ELSE IF( wntqs )
THEN
433 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mn )
434 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
435 maxwrk = max( wrkbl, 3*n + bdspac )
436 minwrk = 3*n + max( m, bdspac )
437 ELSE IF( wntqa )
THEN
439 wrkbl = max( wrkbl, 3*n + lwork_dormbr_qln_mm )
440 wrkbl = max( wrkbl, 3*n + lwork_dormbr_prt_nn )
441 maxwrk = max( wrkbl, 3*n + bdspac )
442 minwrk = 3*n + max( m, bdspac )
445 ELSE IF( minmn.GT.0 )
THEN
458 CALL dgebrd( m, n, dum(1), m, dum(1), dum(1), dum(1),
459 $ dum(1), dum(1), -1, ierr )
460 lwork_dgebrd_mn = int( dum(1) )
462 CALL dgebrd( m, m, a, m, s, dum(1), dum(1),
463 $ dum(1), dum(1), -1, ierr )
464 lwork_dgebrd_mm = int( dum(1) )
466 CALL dgelqf( m, n, a, m, dum(1), dum(1), -1, ierr )
467 lwork_dgelqf_mn = int( dum(1) )
469 CALL dorglq( n, n, m, dum(1), n, dum(1), dum(1), -1, ierr )
470 lwork_dorglq_nn = int( dum(1) )
472 CALL dorglq( m, n, m, a, m, dum(1), dum(1), -1, ierr )
473 lwork_dorglq_mn = int( dum(1) )
475 CALL dorgbr(
'P', m, m, m, a, n, dum(1), dum(1), -1, ierr )
476 lwork_dorgbr_p_mm = int( dum(1) )
478 CALL dormbr(
'P',
'R',
'T', m, m, m, dum(1), m,
479 $ dum(1), dum(1), m, dum(1), -1, ierr )
480 lwork_dormbr_prt_mm = int( dum(1) )
482 CALL dormbr(
'P',
'R',
'T', m, n, m, dum(1), m,
483 $ dum(1), dum(1), m, dum(1), -1, ierr )
484 lwork_dormbr_prt_mn = int( dum(1) )
486 CALL dormbr(
'P',
'R',
'T', n, n, m, dum(1), n,
487 $ dum(1), dum(1), n, dum(1), -1, ierr )
488 lwork_dormbr_prt_nn = int( dum(1) )
490 CALL dormbr(
'Q',
'L',
'N', m, m, m, dum(1), m,
491 $ dum(1), dum(1), m, dum(1), -1, ierr )
492 lwork_dormbr_qln_mm = int( dum(1) )
494 IF( n.GE.mnthr )
THEN
499 wrkbl = m + lwork_dgelqf_mn
500 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
501 maxwrk = max( wrkbl, bdspac + m )
503 ELSE IF( wntqo )
THEN
507 wrkbl = m + lwork_dgelqf_mn
508 wrkbl = max( wrkbl, m + lwork_dorglq_mn )
509 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
510 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
511 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
512 wrkbl = max( wrkbl, 3*m + bdspac )
513 maxwrk = wrkbl + 2*m*m
514 minwrk = bdspac + 2*m*m + 3*m
515 ELSE IF( wntqs )
THEN
519 wrkbl = m + lwork_dgelqf_mn
520 wrkbl = max( wrkbl, m + lwork_dorglq_mn )
521 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
522 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
523 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
524 wrkbl = max( wrkbl, 3*m + bdspac )
526 minwrk = bdspac + m*m + 3*m
527 ELSE IF( wntqa )
THEN
531 wrkbl = m + lwork_dgelqf_mn
532 wrkbl = max( wrkbl, m + lwork_dorglq_nn )
533 wrkbl = max( wrkbl, 3*m + lwork_dgebrd_mm )
534 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
535 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mm )
536 wrkbl = max( wrkbl, 3*m + bdspac )
538 minwrk = m*m + max( 3*m + bdspac, m + n )
544 wrkbl = 3*m + lwork_dgebrd_mn
547 maxwrk = max( wrkbl, 3*m + bdspac )
548 minwrk = 3*m + max( n, bdspac )
549 ELSE IF( wntqo )
THEN
551 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
552 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mn )
553 wrkbl = max( wrkbl, 3*m + bdspac )
555 minwrk = 3*m + max( n, m*m + bdspac )
556 ELSE IF( wntqs )
THEN
558 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
559 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_mn )
560 maxwrk = max( wrkbl, 3*m + bdspac )
561 minwrk = 3*m + max( n, bdspac )
562 ELSE IF( wntqa )
THEN
564 wrkbl = max( wrkbl, 3*m + lwork_dormbr_qln_mm )
565 wrkbl = max( wrkbl, 3*m + lwork_dormbr_prt_nn )
566 maxwrk = max( wrkbl, 3*m + bdspac )
567 minwrk = 3*m + max( n, bdspac )
572 maxwrk = max( maxwrk, minwrk )
575 IF( lwork.LT.minwrk .AND. .NOT.lquery )
THEN
581 CALL xerbla(
'DGESDD', -info )
583 ELSE IF( lquery )
THEN
589 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
596 smlnum = sqrt( dlamch(
'S' ) ) / eps
597 bignum = one / smlnum
601 anrm = dlange(
'M', m, n, a, lda, dum )
603 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
605 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
606 ELSE IF( anrm.GT.bignum )
THEN
608 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
617 IF( m.GE.mnthr )
THEN
631 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
632 $ lwork - nwork + 1, ierr )
636 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
646 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
647 $ work( itaup ), work( nwork ), lwork-nwork+1,
654 CALL dbdsdc(
'U',
'N', n, s, work( ie ), dum, 1, dum, 1,
655 $ dum, idum, work( nwork ), iwork, info )
657 ELSE IF( wntqo )
THEN
667 IF( lwork .GE. lda*n + n*n + 3*n + bdspac )
THEN
670 ldwrkr = ( lwork - n*n - 3*n - bdspac ) / n
679 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
680 $ lwork - nwork + 1, ierr )
684 CALL dlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
685 CALL dlaset(
'L', n - 1, n - 1, zero, zero, work(ir+1),
692 CALL dorgqr( m, n, n, a, lda, work( itau ),
693 $ work( nwork ), lwork - nwork + 1, ierr )
703 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
704 $ work( itauq ), work( itaup ), work( nwork ),
705 $ lwork - nwork + 1, ierr )
717 CALL dbdsdc(
'U',
'I', n, s, work( ie ), work( iu ), n,
718 $ vt, ldvt, dum, idum, work( nwork ), iwork,
726 CALL dormbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
727 $ work( itauq ), work( iu ), n, work( nwork ),
728 $ lwork - nwork + 1, ierr )
729 CALL dormbr(
'P',
'R',
'T', n, n, n, work( ir ), ldwrkr,
730 $ work( itaup ), vt, ldvt, work( nwork ),
731 $ lwork - nwork + 1, ierr )
738 DO 10 i = 1, m, ldwrkr
739 chunk = min( m - i + 1, ldwrkr )
740 CALL dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
741 $ lda, work( iu ), n, zero, work( ir ),
743 CALL dlacpy(
'F', chunk, n, work( ir ), ldwrkr,
747 ELSE IF( wntqs )
THEN
765 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
766 $ lwork - nwork + 1, ierr )
770 CALL dlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
771 CALL dlaset(
'L', n - 1, n - 1, zero, zero, work(ir+1),
778 CALL dorgqr( m, n, n, a, lda, work( itau ),
779 $ work( nwork ), lwork - nwork + 1, ierr )
789 CALL dgebrd( n, n, work( ir ), ldwrkr, s, work( ie ),
790 $ work( itauq ), work( itaup ), work( nwork ),
791 $ lwork - nwork + 1, ierr )
798 CALL dbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
799 $ ldvt, dum, idum, work( nwork ), iwork,
807 CALL dormbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
808 $ work( itauq ), u, ldu, work( nwork ),
809 $ lwork - nwork + 1, ierr )
811 CALL dormbr(
'P',
'R',
'T', n, n, n, work( ir ), ldwrkr,
812 $ work( itaup ), vt, ldvt, work( nwork ),
813 $ lwork - nwork + 1, ierr )
819 CALL dlacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
820 CALL dgemm(
'N',
'N', m, n, n, one, a, lda, work( ir ),
821 $ ldwrkr, zero, u, ldu )
823 ELSE IF( wntqa )
THEN
841 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
842 $ lwork - nwork + 1, ierr )
843 CALL dlacpy(
'L', m, n, a, lda, u, ldu )
848 CALL dorgqr( m, m, n, u, ldu, work( itau ),
849 $ work( nwork ), lwork - nwork + 1, ierr )
853 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
863 CALL dgebrd( n, n, a, lda, s, work( ie ), work( itauq ),
864 $ work( itaup ), work( nwork ), lwork-nwork+1,
872 CALL dbdsdc(
'U',
'I', n, s, work( ie ), work( iu ), n,
873 $ vt, ldvt, dum, idum, work( nwork ), iwork,
881 CALL dormbr(
'Q',
'L',
'N', n, n, n, a, lda,
882 $ work( itauq ), work( iu ), ldwrku,
883 $ work( nwork ), lwork - nwork + 1, ierr )
884 CALL dormbr(
'P',
'R',
'T', n, n, n, a, lda,
885 $ work( itaup ), vt, ldvt, work( nwork ),
886 $ lwork - nwork + 1, ierr )
892 CALL dgemm(
'N',
'N', m, n, n, one, u, ldu, work( iu ),
893 $ ldwrku, zero, a, lda )
897 CALL dlacpy(
'F', m, n, a, lda, u, ldu )
917 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
918 $ work( itaup ), work( nwork ), lwork-nwork+1,
926 CALL dbdsdc(
'U',
'N', n, s, work( ie ), dum, 1, dum, 1,
927 $ dum, idum, work( nwork ), iwork, info )
928 ELSE IF( wntqo )
THEN
931 IF( lwork .GE. m*n + 3*n + bdspac )
THEN
936 nwork = iu + ldwrku*n
937 CALL dlaset(
'F', m, n, zero, zero, work( iu ),
946 nwork = iu + ldwrku*n
951 ldwrkr = ( lwork - n*n - 3*n ) / n
953 nwork = iu + ldwrku*n
960 CALL dbdsdc(
'U',
'I', n, s, work( ie ), work( iu ),
961 $ ldwrku, vt, ldvt, dum, idum, work( nwork ),
968 CALL dormbr(
'P',
'R',
'T', n, n, n, a, lda,
969 $ work( itaup ), vt, ldvt, work( nwork ),
970 $ lwork - nwork + 1, ierr )
972 IF( lwork .GE. m*n + 3*n + bdspac )
THEN
979 CALL dormbr(
'Q',
'L',
'N', m, n, n, a, lda,
980 $ work( itauq ), work( iu ), ldwrku,
981 $ work( nwork ), lwork - nwork + 1, ierr )
985 CALL dlacpy(
'F', m, n, work( iu ), ldwrku, a, lda )
993 CALL dorgbr(
'Q', m, n, n, a, lda, work( itauq ),
994 $ work( nwork ), lwork - nwork + 1, ierr )
1002 DO 20 i = 1, m, ldwrkr
1003 chunk = min( m - i + 1, ldwrkr )
1004 CALL dgemm(
'N',
'N', chunk, n, n, one, a( i, 1 ),
1005 $ lda, work( iu ), ldwrku, zero,
1006 $ work( ir ), ldwrkr )
1007 CALL dlacpy(
'F', chunk, n, work( ir ), ldwrkr,
1012 ELSE IF( wntqs )
THEN
1020 CALL dlaset(
'F', m, n, zero, zero, u, ldu )
1021 CALL dbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
1022 $ ldvt, dum, idum, work( nwork ), iwork,
1030 CALL dormbr(
'Q',
'L',
'N', m, n, n, a, lda,
1031 $ work( itauq ), u, ldu, work( nwork ),
1032 $ lwork - nwork + 1, ierr )
1033 CALL dormbr(
'P',
'R',
'T', n, n, n, a, lda,
1034 $ work( itaup ), vt, ldvt, work( nwork ),
1035 $ lwork - nwork + 1, ierr )
1036 ELSE IF( wntqa )
THEN
1044 CALL dlaset(
'F', m, m, zero, zero, u, ldu )
1045 CALL dbdsdc(
'U',
'I', n, s, work( ie ), u, ldu, vt,
1046 $ ldvt, dum, idum, work( nwork ), iwork,
1052 CALL dlaset(
'F', m - n, m - n, zero, one, u(n+1,n+1),
1061 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1062 $ work( itauq ), u, ldu, work( nwork ),
1063 $ lwork - nwork + 1, ierr )
1064 CALL dormbr(
'P',
'R',
'T', n, n, m, a, lda,
1065 $ work( itaup ), vt, ldvt, work( nwork ),
1066 $ lwork - nwork + 1, ierr )
1077 IF( n.GE.mnthr )
THEN
1091 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1092 $ lwork - nwork + 1, ierr )
1096 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1106 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1107 $ work( itaup ), work( nwork ), lwork-nwork+1,
1114 CALL dbdsdc(
'U',
'N', m, s, work( ie ), dum, 1, dum, 1,
1115 $ dum, idum, work( nwork ), iwork, info )
1117 ELSE IF( wntqo )
THEN
1129 IF( lwork .GE. m*n + m*m + 3*m + bdspac )
THEN
1134 chunk = ( lwork - m*m ) / m
1136 itau = il + ldwrkl*m
1143 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1144 $ lwork - nwork + 1, ierr )
1148 CALL dlacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1149 CALL dlaset(
'U', m - 1, m - 1, zero, zero,
1150 $ work( il + ldwrkl ), ldwrkl )
1156 CALL dorglq( m, n, m, a, lda, work( itau ),
1157 $ work( nwork ), lwork - nwork + 1, ierr )
1167 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1168 $ work( itauq ), work( itaup ), work( nwork ),
1169 $ lwork - nwork + 1, ierr )
1176 CALL dbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1177 $ work( ivt ), m, dum, idum, work( nwork ),
1185 CALL dormbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1186 $ work( itauq ), u, ldu, work( nwork ),
1187 $ lwork - nwork + 1, ierr )
1188 CALL dormbr(
'P',
'R',
'T', m, m, m, work( il ), ldwrkl,
1189 $ work( itaup ), work( ivt ), m,
1190 $ work( nwork ), lwork - nwork + 1, ierr )
1198 DO 30 i = 1, n, chunk
1199 blk = min( n - i + 1, chunk )
1200 CALL dgemm(
'N',
'N', m, blk, m, one, work( ivt ), m,
1201 $ a( 1, i ), lda, zero, work( il ), ldwrkl )
1202 CALL dlacpy(
'F', m, blk, work( il ), ldwrkl,
1206 ELSE IF( wntqs )
THEN
1217 itau = il + ldwrkl*m
1224 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1225 $ lwork - nwork + 1, ierr )
1229 CALL dlacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1230 CALL dlaset(
'U', m - 1, m - 1, zero, zero,
1231 $ work( il + ldwrkl ), ldwrkl )
1237 CALL dorglq( m, n, m, a, lda, work( itau ),
1238 $ work( nwork ), lwork - nwork + 1, ierr )
1248 CALL dgebrd( m, m, work( il ), ldwrkl, s, work( ie ),
1249 $ work( itauq ), work( itaup ), work( nwork ),
1250 $ lwork - nwork + 1, ierr )
1257 CALL dbdsdc(
'U',
'I', m, s, work( ie ), u, ldu, vt,
1258 $ ldvt, dum, idum, work( nwork ), iwork,
1266 CALL dormbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1267 $ work( itauq ), u, ldu, work( nwork ),
1268 $ lwork - nwork + 1, ierr )
1269 CALL dormbr(
'P',
'R',
'T', m, m, m, work( il ), ldwrkl,
1270 $ work( itaup ), vt, ldvt, work( nwork ),
1271 $ lwork - nwork + 1, ierr )
1277 CALL dlacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1278 CALL dgemm(
'N',
'N', m, n, m, one, work( il ), ldwrkl,
1279 $ a, lda, zero, vt, ldvt )
1281 ELSE IF( wntqa )
THEN
1292 itau = ivt + ldwkvt*m
1299 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
1300 $ lwork - nwork + 1, ierr )
1301 CALL dlacpy(
'U', m, n, a, lda, vt, ldvt )
1307 CALL dorglq( n, n, m, vt, ldvt, work( itau ),
1308 $ work( nwork ), lwork - nwork + 1, ierr )
1312 CALL dlaset(
'U', m-1, m-1, zero, zero, a( 1, 2 ), lda )
1322 CALL dgebrd( m, m, a, lda, s, work( ie ), work( itauq ),
1323 $ work( itaup ), work( nwork ), lwork-nwork+1,
1331 CALL dbdsdc(
'U',
'I', m, s, work( ie ), u, ldu,
1332 $ work( ivt ), ldwkvt, dum, idum,
1333 $ work( nwork ), iwork, info )
1340 CALL dormbr(
'Q',
'L',
'N', m, m, m, a, lda,
1341 $ work( itauq ), u, ldu, work( nwork ),
1342 $ lwork - nwork + 1, ierr )
1343 CALL dormbr(
'P',
'R',
'T', m, m, m, a, lda,
1344 $ work( itaup ), work( ivt ), ldwkvt,
1345 $ work( nwork ), lwork - nwork + 1, ierr )
1351 CALL dgemm(
'N',
'N', m, n, m, one, work( ivt ), ldwkvt,
1352 $ vt, ldvt, zero, a, lda )
1356 CALL dlacpy(
'F', m, n, a, lda, vt, ldvt )
1376 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
1377 $ work( itaup ), work( nwork ), lwork-nwork+1,
1385 CALL dbdsdc(
'L',
'N', m, s, work( ie ), dum, 1, dum, 1,
1386 $ dum, idum, work( nwork ), iwork, info )
1387 ELSE IF( wntqo )
THEN
1391 IF( lwork .GE. m*n + 3*m + bdspac )
THEN
1395 CALL dlaset(
'F', m, n, zero, zero, work( ivt ),
1397 nwork = ivt + ldwkvt*n
1404 nwork = ivt + ldwkvt*m
1409 chunk = ( lwork - m*m - 3*m ) / m
1417 CALL dbdsdc(
'L',
'I', m, s, work( ie ), u, ldu,
1418 $ work( ivt ), ldwkvt, dum, idum,
1419 $ work( nwork ), iwork, info )
1425 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1426 $ work( itauq ), u, ldu, work( nwork ),
1427 $ lwork - nwork + 1, ierr )
1429 IF( lwork .GE. m*n + 3*m + bdspac )
THEN
1436 CALL dormbr(
'P',
'R',
'T', m, n, m, a, lda,
1437 $ work( itaup ), work( ivt ), ldwkvt,
1438 $ work( nwork ), lwork - nwork + 1, ierr )
1442 CALL dlacpy(
'F', m, n, work( ivt ), ldwkvt, a, lda )
1450 CALL dorgbr(
'P', m, n, m, a, lda, work( itaup ),
1451 $ work( nwork ), lwork - nwork + 1, ierr )
1459 DO 40 i = 1, n, chunk
1460 blk = min( n - i + 1, chunk )
1461 CALL dgemm(
'N',
'N', m, blk, m, one, work( ivt ),
1462 $ ldwkvt, a( 1, i ), lda, zero,
1464 CALL dlacpy(
'F', m, blk, work( il ), m, a( 1, i ),
1468 ELSE IF( wntqs )
THEN
1476 CALL dlaset(
'F', m, n, zero, zero, vt, ldvt )
1477 CALL dbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1478 $ ldvt, dum, idum, work( nwork ), iwork,
1486 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1487 $ work( itauq ), u, ldu, work( nwork ),
1488 $ lwork - nwork + 1, ierr )
1489 CALL dormbr(
'P',
'R',
'T', m, n, m, a, lda,
1490 $ work( itaup ), vt, ldvt, work( nwork ),
1491 $ lwork - nwork + 1, ierr )
1492 ELSE IF( wntqa )
THEN
1500 CALL dlaset(
'F', n, n, zero, zero, vt, ldvt )
1501 CALL dbdsdc(
'L',
'I', m, s, work( ie ), u, ldu, vt,
1502 $ ldvt, dum, idum, work( nwork ), iwork,
1508 CALL dlaset(
'F', n-m, n-m, zero, one, vt(m+1,m+1),
1517 CALL dormbr(
'Q',
'L',
'N', m, m, n, a, lda,
1518 $ work( itauq ), u, ldu, work( nwork ),
1519 $ lwork - nwork + 1, ierr )
1520 CALL dormbr(
'P',
'R',
'T', n, n, m, a, lda,
1521 $ work( itaup ), vt, ldvt, work( nwork ),
1522 $ lwork - nwork + 1, ierr )
1531 IF( iscl.EQ.1 )
THEN
1532 IF( anrm.GT.bignum )
1533 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
1535 IF( anrm.LT.smlnum )
1536 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
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 dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, IWORK, INFO)
DGESDD
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF