226 SUBROUTINE cgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
227 $ work, lwork, rwork, iwork, info )
237 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
241 REAL RWORK( * ), S( * )
242 COMPLEX A( lda, * ), U( ldu, * ), VT( ldvt, * ),
250 parameter ( czero = ( 0.0e+0, 0.0e+0 ),
251 $ cone = ( 1.0e+0, 0.0e+0 ) )
253 parameter ( zero = 0.0e+0, one = 1.0e+0 )
256 LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
257 INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
258 $ iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
259 $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
260 $ mnthr1, mnthr2, nrwork, nwork, wrkbl
261 INTEGER LWORK_CGEBRD_MN, LWORK_CGEBRD_MM,
262 $ lwork_cgebrd_nn, lwork_cgelqf_mn,
264 $ lwork_cungbr_p_mn, lwork_cungbr_p_nn,
265 $ lwork_cungbr_q_mn, lwork_cungbr_q_mm,
266 $ lwork_cunglq_mn, lwork_cunglq_nn,
267 $ lwork_cungqr_mm, lwork_cungqr_mn,
268 $ lwork_cunmbr_prc_mm, lwork_cunmbr_qln_mm,
269 $ lwork_cunmbr_prc_mn, lwork_cunmbr_qln_mn,
270 $ lwork_cunmbr_prc_nn, lwork_cunmbr_qln_nn
271 REAL ANRM, BIGNUM, EPS, SMLNUM
286 EXTERNAL lsame, slamch, clange
289 INTRINSIC int, max, min, sqrt
297 mnthr1 = int( minmn*17.0e0 / 9.0e0 )
298 mnthr2 = int( minmn*5.0e0 / 3.0e0 )
299 wntqa = lsame( jobz,
'A' )
300 wntqs = lsame( jobz,
'S' )
301 wntqas = wntqa .OR. wntqs
302 wntqo = lsame( jobz,
'O' )
303 wntqn = lsame( jobz,
'N' )
304 lquery = ( lwork.EQ.-1 )
308 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) )
THEN
310 ELSE IF( m.LT.0 )
THEN
312 ELSE IF( n.LT.0 )
THEN
314 ELSE IF( lda.LT.max( 1, m ) )
THEN
316 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
317 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) )
THEN
319 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
320 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
321 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) )
THEN
336 IF( m.GE.n .AND. minmn.GT.0 )
THEN
345 CALL cgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
346 $ cdum(1), cdum(1), -1, ierr )
347 lwork_cgebrd_mn = int( cdum(1) )
349 CALL cgebrd( n, n, cdum(1), n, dum(1), dum(1), cdum(1),
350 $ cdum(1), cdum(1), -1, ierr )
351 lwork_cgebrd_nn = int( cdum(1) )
353 CALL cgeqrf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
354 lwork_cgeqrf_mn = int( cdum(1) )
356 CALL cungbr(
'P', n, n, n, cdum(1), n, cdum(1), cdum(1),
358 lwork_cungbr_p_nn = int( cdum(1) )
360 CALL cungbr(
'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
362 lwork_cungbr_q_mm = int( cdum(1) )
364 CALL cungbr(
'Q', m, n, n, cdum(1), m, cdum(1), cdum(1),
366 lwork_cungbr_q_mn = int( cdum(1) )
368 CALL cungqr( m, m, n, cdum(1), m, cdum(1), cdum(1),
370 lwork_cungqr_mm = int( cdum(1) )
372 CALL cungqr( m, n, n, cdum(1), m, cdum(1), cdum(1),
374 lwork_cungqr_mn = int( cdum(1) )
376 CALL cunmbr(
'P',
'R',
'C', n, n, n, cdum(1), n, cdum(1),
377 $ cdum(1), n, cdum(1), -1, ierr )
378 lwork_cunmbr_prc_nn = int( cdum(1) )
380 CALL cunmbr(
'Q',
'L',
'N', m, m, n, cdum(1), m, cdum(1),
381 $ cdum(1), m, cdum(1), -1, ierr )
382 lwork_cunmbr_qln_mm = int( cdum(1) )
384 CALL cunmbr(
'Q',
'L',
'N', m, n, n, cdum(1), m, cdum(1),
385 $ cdum(1), m, cdum(1), -1, ierr )
386 lwork_cunmbr_qln_mn = int( cdum(1) )
388 CALL cunmbr(
'Q',
'L',
'N', n, n, n, cdum(1), n, cdum(1),
389 $ cdum(1), n, cdum(1), -1, ierr )
390 lwork_cunmbr_qln_nn = int( cdum(1) )
392 IF( m.GE.mnthr1 )
THEN
397 maxwrk = n + lwork_cgeqrf_mn
398 maxwrk = max( maxwrk, 2*n + lwork_cgebrd_nn )
400 ELSE IF( wntqo )
THEN
404 wrkbl = n + lwork_cgeqrf_mn
405 wrkbl = max( wrkbl, n + lwork_cungqr_mn )
406 wrkbl = max( wrkbl, 2*n + lwork_cgebrd_nn )
407 wrkbl = max( wrkbl, 2*n + lwork_cunmbr_qln_nn )
408 wrkbl = max( wrkbl, 2*n + lwork_cunmbr_prc_nn )
409 maxwrk = m*n + n*n + wrkbl
411 ELSE IF( wntqs )
THEN
415 wrkbl = n + lwork_cgeqrf_mn
416 wrkbl = max( wrkbl, n + lwork_cungqr_mn )
417 wrkbl = max( wrkbl, 2*n + lwork_cgebrd_nn )
418 wrkbl = max( wrkbl, 2*n + lwork_cunmbr_qln_nn )
419 wrkbl = max( wrkbl, 2*n + lwork_cunmbr_prc_nn )
422 ELSE IF( wntqa )
THEN
426 wrkbl = n + lwork_cgeqrf_mn
427 wrkbl = max( wrkbl, n + lwork_cungqr_mm )
428 wrkbl = max( wrkbl, 2*n + lwork_cgebrd_nn )
429 wrkbl = max( wrkbl, 2*n + lwork_cunmbr_qln_nn )
430 wrkbl = max( wrkbl, 2*n + lwork_cunmbr_prc_nn )
432 minwrk = n*n + max( 3*n, n + m )
434 ELSE IF( m.GE.mnthr2 )
THEN
438 maxwrk = 2*n + lwork_cgebrd_mn
442 maxwrk = max( maxwrk, 2*n + lwork_cungbr_p_nn )
443 maxwrk = max( maxwrk, 2*n + lwork_cungbr_q_mn )
444 maxwrk = maxwrk + m*n
445 minwrk = minwrk + n*n
446 ELSE IF( wntqs )
THEN
448 maxwrk = max( maxwrk, 2*n + lwork_cungbr_p_nn )
449 maxwrk = max( maxwrk, 2*n + lwork_cungbr_q_mn )
450 ELSE IF( wntqa )
THEN
452 maxwrk = max( maxwrk, 2*n + lwork_cungbr_p_nn )
453 maxwrk = max( maxwrk, 2*n + lwork_cungbr_q_mm )
459 maxwrk = 2*n + lwork_cgebrd_mn
463 maxwrk = max( maxwrk, 2*n + lwork_cunmbr_prc_nn )
464 maxwrk = max( maxwrk, 2*n + lwork_cunmbr_qln_mn )
465 maxwrk = maxwrk + m*n
466 minwrk = minwrk + n*n
467 ELSE IF( wntqs )
THEN
469 maxwrk = max( maxwrk, 2*n + lwork_cunmbr_qln_mn )
470 maxwrk = max( maxwrk, 2*n + lwork_cunmbr_prc_nn )
471 ELSE IF( wntqa )
THEN
473 maxwrk = max( maxwrk, 2*n + lwork_cunmbr_qln_mm )
474 maxwrk = max( maxwrk, 2*n + lwork_cunmbr_prc_nn )
477 ELSE IF( minmn.GT.0 )
THEN
486 CALL cgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
487 $ cdum(1), cdum(1), -1, ierr )
488 lwork_cgebrd_mn = int( cdum(1) )
490 CALL cgebrd( m, m, cdum(1), m, dum(1), dum(1), cdum(1),
491 $ cdum(1), cdum(1), -1, ierr )
492 lwork_cgebrd_mm = int( cdum(1) )
494 CALL cgelqf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
495 lwork_cgelqf_mn = int( cdum(1) )
497 CALL cungbr(
'P', m, n, m, cdum(1), m, cdum(1), cdum(1),
499 lwork_cungbr_p_mn = int( cdum(1) )
501 CALL cungbr(
'P', n, n, m, cdum(1), n, cdum(1), cdum(1),
503 lwork_cungbr_p_nn = int( cdum(1) )
505 CALL cungbr(
'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
507 lwork_cungbr_q_mm = int( cdum(1) )
509 CALL cunglq( m, n, m, cdum(1), m, cdum(1), cdum(1),
511 lwork_cunglq_mn = int( cdum(1) )
513 CALL cunglq( n, n, m, cdum(1), n, cdum(1), cdum(1),
515 lwork_cunglq_nn = int( cdum(1) )
517 CALL cunmbr(
'P',
'R',
'C', m, m, m, cdum(1), m, cdum(1),
518 $ cdum(1), m, cdum(1), -1, ierr )
519 lwork_cunmbr_prc_mm = int( cdum(1) )
521 CALL cunmbr(
'P',
'R',
'C', m, n, m, cdum(1), m, cdum(1),
522 $ cdum(1), m, cdum(1), -1, ierr )
523 lwork_cunmbr_prc_mn = int( cdum(1) )
525 CALL cunmbr(
'P',
'R',
'C', n, n, m, cdum(1), n, cdum(1),
526 $ cdum(1), n, cdum(1), -1, ierr )
527 lwork_cunmbr_prc_nn = int( cdum(1) )
529 CALL cunmbr(
'Q',
'L',
'N', m, m, m, cdum(1), m, cdum(1),
530 $ cdum(1), m, cdum(1), -1, ierr )
531 lwork_cunmbr_qln_mm = int( cdum(1) )
533 IF( n.GE.mnthr1 )
THEN
538 maxwrk = m + lwork_cgelqf_mn
539 maxwrk = max( maxwrk, 2*m + lwork_cgebrd_mm )
541 ELSE IF( wntqo )
THEN
545 wrkbl = m + lwork_cgelqf_mn
546 wrkbl = max( wrkbl, m + lwork_cunglq_mn )
547 wrkbl = max( wrkbl, 2*m + lwork_cgebrd_mm )
548 wrkbl = max( wrkbl, 2*m + lwork_cunmbr_qln_mm )
549 wrkbl = max( wrkbl, 2*m + lwork_cunmbr_prc_mm )
550 maxwrk = m*n + m*m + wrkbl
552 ELSE IF( wntqs )
THEN
556 wrkbl = m + lwork_cgelqf_mn
557 wrkbl = max( wrkbl, m + lwork_cunglq_mn )
558 wrkbl = max( wrkbl, 2*m + lwork_cgebrd_mm )
559 wrkbl = max( wrkbl, 2*m + lwork_cunmbr_qln_mm )
560 wrkbl = max( wrkbl, 2*m + lwork_cunmbr_prc_mm )
563 ELSE IF( wntqa )
THEN
567 wrkbl = m + lwork_cgelqf_mn
568 wrkbl = max( wrkbl, m + lwork_cunglq_nn )
569 wrkbl = max( wrkbl, 2*m + lwork_cgebrd_mm )
570 wrkbl = max( wrkbl, 2*m + lwork_cunmbr_qln_mm )
571 wrkbl = max( wrkbl, 2*m + lwork_cunmbr_prc_mm )
573 minwrk = m*m + max( 3*m, m + n )
575 ELSE IF( n.GE.mnthr2 )
THEN
579 maxwrk = 2*m + lwork_cgebrd_mn
583 maxwrk = max( maxwrk, 2*m + lwork_cungbr_q_mm )
584 maxwrk = max( maxwrk, 2*m + lwork_cungbr_p_mn )
585 maxwrk = maxwrk + m*n
586 minwrk = minwrk + m*m
587 ELSE IF( wntqs )
THEN
589 maxwrk = max( maxwrk, 2*m + lwork_cungbr_q_mm )
590 maxwrk = max( maxwrk, 2*m + lwork_cungbr_p_mn )
591 ELSE IF( wntqa )
THEN
593 maxwrk = max( maxwrk, 2*m + lwork_cungbr_q_mm )
594 maxwrk = max( maxwrk, 2*m + lwork_cungbr_p_nn )
600 maxwrk = 2*m + lwork_cgebrd_mn
604 maxwrk = max( maxwrk, 2*m + lwork_cunmbr_qln_mm )
605 maxwrk = max( maxwrk, 2*m + lwork_cunmbr_prc_mn )
606 maxwrk = maxwrk + m*n
607 minwrk = minwrk + m*m
608 ELSE IF( wntqs )
THEN
610 maxwrk = max( maxwrk, 2*m + lwork_cunmbr_qln_mm )
611 maxwrk = max( maxwrk, 2*m + lwork_cunmbr_prc_mn )
612 ELSE IF( wntqa )
THEN
614 maxwrk = max( maxwrk, 2*m + lwork_cunmbr_qln_mm )
615 maxwrk = max( maxwrk, 2*m + lwork_cunmbr_prc_nn )
619 maxwrk = max( maxwrk, minwrk )
623 IF( lwork.LT.minwrk .AND. .NOT. lquery )
THEN
629 CALL xerbla(
'CGESDD', -info )
631 ELSE IF( lquery )
THEN
637 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
644 smlnum = sqrt( slamch(
'S' ) ) / eps
645 bignum = one / smlnum
649 anrm = clange(
'M', m, n, a, lda, dum )
651 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
653 CALL clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
654 ELSE IF( anrm.GT.bignum )
THEN
656 CALL clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
665 IF( m.GE.mnthr1 )
THEN
680 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
681 $ lwork-nwork+1, ierr )
685 CALL claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
697 CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
698 $ work( itaup ), work( nwork ), lwork-nwork+1,
706 CALL sbdsdc(
'U',
'N', n, s, rwork( ie ), dum,1,dum,1,
707 $ dum, idum, rwork( nrwork ), iwork, info )
709 ELSE IF( wntqo )
THEN
721 IF( lwork .GE. m*n + n*n + 3*n )
THEN
727 ldwrkr = ( lwork - n*n - 3*n ) / n
737 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
738 $ lwork-nwork+1, ierr )
742 CALL clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
743 CALL claset(
'L', n-1, n-1, czero, czero, work( ir+1 ),
751 CALL cungqr( m, n, n, a, lda, work( itau ),
752 $ work( nwork ), lwork-nwork+1, ierr )
763 CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
764 $ work( itauq ), work( itaup ), work( nwork ),
765 $ lwork-nwork+1, ierr )
776 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
777 $ n, rwork( irvt ), n, dum, idum,
778 $ rwork( nrwork ), iwork, info )
786 CALL clacp2(
'F', n, n, rwork( iru ), n, work( iu ),
788 CALL cunmbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
789 $ work( itauq ), work( iu ), ldwrku,
790 $ work( nwork ), lwork-nwork+1, ierr )
798 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
799 CALL cunmbr(
'P',
'R',
'C', n, n, n, work( ir ), ldwrkr,
800 $ work( itaup ), vt, ldvt, work( nwork ),
801 $ lwork-nwork+1, ierr )
809 DO 10 i = 1, m, ldwrkr
810 chunk = min( m-i+1, ldwrkr )
811 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
812 $ lda, work( iu ), ldwrku, czero,
813 $ work( ir ), ldwrkr )
814 CALL clacpy(
'F', chunk, n, work( ir ), ldwrkr,
818 ELSE IF( wntqs )
THEN
837 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
838 $ lwork-nwork+1, ierr )
842 CALL clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
843 CALL claset(
'L', n-1, n-1, czero, czero, work( ir+1 ),
851 CALL cungqr( m, n, n, a, lda, work( itau ),
852 $ work( nwork ), lwork-nwork+1, ierr )
863 CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
864 $ work( itauq ), work( itaup ), work( nwork ),
865 $ lwork-nwork+1, ierr )
876 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
877 $ n, rwork( irvt ), n, dum, idum,
878 $ rwork( nrwork ), iwork, info )
886 CALL clacp2(
'F', n, n, rwork( iru ), n, u, ldu )
887 CALL cunmbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
888 $ work( itauq ), u, ldu, work( nwork ),
889 $ lwork-nwork+1, ierr )
897 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
898 CALL cunmbr(
'P',
'R',
'C', n, n, n, work( ir ), ldwrkr,
899 $ work( itaup ), vt, ldvt, work( nwork ),
900 $ lwork-nwork+1, ierr )
907 CALL clacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
908 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda, work( ir ),
909 $ ldwrkr, czero, u, ldu )
911 ELSE IF( wntqa )
THEN
930 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
931 $ lwork-nwork+1, ierr )
932 CALL clacpy(
'L', m, n, a, lda, u, ldu )
939 CALL cungqr( m, m, n, u, ldu, work( itau ),
940 $ work( nwork ), lwork-nwork+1, ierr )
944 CALL claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
956 CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
957 $ work( itaup ), work( nwork ), lwork-nwork+1,
969 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
970 $ n, rwork( irvt ), n, dum, idum,
971 $ rwork( nrwork ), iwork, info )
979 CALL clacp2(
'F', n, n, rwork( iru ), n, work( iu ),
981 CALL cunmbr(
'Q',
'L',
'N', n, n, n, a, lda,
982 $ work( itauq ), work( iu ), ldwrku,
983 $ work( nwork ), lwork-nwork+1, ierr )
991 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
992 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
993 $ work( itaup ), vt, ldvt, work( nwork ),
994 $ lwork-nwork+1, ierr )
1001 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu, work( iu ),
1002 $ ldwrku, czero, a, lda )
1006 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1010 ELSE IF( m.GE.mnthr2 )
THEN
1029 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1030 $ work( itaup ), work( nwork ), lwork-nwork+1,
1039 CALL sbdsdc(
'U',
'N', n, s, rwork( ie ), dum, 1,dum,1,
1040 $ dum, idum, rwork( nrwork ), iwork, info )
1041 ELSE IF( wntqo )
THEN
1053 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1054 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1055 $ work( nwork ), lwork-nwork+1, ierr )
1062 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
1063 $ work( nwork ), lwork-nwork+1, ierr )
1065 IF( lwork .GE. m*n + 3*n )
THEN
1074 ldwrku = ( lwork - 3*n ) / n
1076 nwork = iu + ldwrku*n
1084 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1085 $ n, rwork( irvt ), n, dum, idum,
1086 $ rwork( nrwork ), iwork, info )
1093 CALL clarcm( n, n, rwork( irvt ), n, vt, ldvt,
1094 $ work( iu ), ldwrku, rwork( nrwork ) )
1095 CALL clacpy(
'F', n, n, work( iu ), ldwrku, vt, ldvt )
1105 DO 20 i = 1, m, ldwrku
1106 chunk = min( m-i+1, ldwrku )
1107 CALL clacrm( chunk, n, a( i, 1 ), lda, rwork( iru ),
1108 $ n, work( iu ), ldwrku, rwork( nrwork ) )
1109 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
1113 ELSE IF( wntqs )
THEN
1121 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1122 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1123 $ work( nwork ), lwork-nwork+1, ierr )
1130 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1131 CALL cungbr(
'Q', m, n, n, u, ldu, work( itauq ),
1132 $ work( nwork ), lwork-nwork+1, ierr )
1143 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1144 $ n, rwork( irvt ), n, dum, idum,
1145 $ rwork( nrwork ), iwork, info )
1152 CALL clarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1154 CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
1162 CALL clacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1164 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1173 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1174 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1175 $ work( nwork ), lwork-nwork+1, ierr )
1182 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1183 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1184 $ work( nwork ), lwork-nwork+1, ierr )
1195 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1196 $ n, rwork( irvt ), n, dum, idum,
1197 $ rwork( nrwork ), iwork, info )
1204 CALL clarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1206 CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
1214 CALL clacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1216 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1238 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1239 $ work( itaup ), work( nwork ), lwork-nwork+1,
1248 CALL sbdsdc(
'U',
'N', n, s, rwork( ie ), dum,1,dum,1,
1249 $ dum, idum, rwork( nrwork ), iwork, info )
1250 ELSE IF( wntqo )
THEN
1255 IF( lwork .GE. m*n + 3*n )
THEN
1264 ldwrku = ( lwork - 3*n ) / n
1266 nwork = iu + ldwrku*n
1275 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1276 $ n, rwork( irvt ), n, dum, idum,
1277 $ rwork( nrwork ), iwork, info )
1285 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1286 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1287 $ work( itaup ), vt, ldvt, work( nwork ),
1288 $ lwork-nwork+1, ierr )
1290 IF( lwork .GE. m*n + 3*n )
THEN
1300 CALL claset(
'F', m, n, czero, czero, work( iu ),
1302 CALL clacp2(
'F', n, n, rwork( iru ), n, work( iu ),
1304 CALL cunmbr(
'Q',
'L',
'N', m, n, n, a, lda,
1305 $ work( itauq ), work( iu ), ldwrku,
1306 $ work( nwork ), lwork-nwork+1, ierr )
1307 CALL clacpy(
'F', m, n, work( iu ), ldwrku, a, lda )
1316 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
1317 $ work( nwork ), lwork-nwork+1, ierr )
1327 DO 30 i = 1, m, ldwrku
1328 chunk = min( m-i+1, ldwrku )
1329 CALL clacrm( chunk, n, a( i, 1 ), lda,
1330 $ rwork( iru ), n, work( iu ), ldwrku,
1332 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
1337 ELSE IF( wntqs )
THEN
1349 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1350 $ n, rwork( irvt ), n, dum, idum,
1351 $ rwork( nrwork ), iwork, info )
1359 CALL claset(
'F', m, n, czero, czero, u, ldu )
1360 CALL clacp2(
'F', n, n, rwork( iru ), n, u, ldu )
1361 CALL cunmbr(
'Q',
'L',
'N', m, n, n, a, lda,
1362 $ work( itauq ), u, ldu, work( nwork ),
1363 $ lwork-nwork+1, ierr )
1371 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1372 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1373 $ work( itaup ), vt, ldvt, work( nwork ),
1374 $ lwork-nwork+1, ierr )
1387 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1388 $ n, rwork( irvt ), n, dum, idum,
1389 $ rwork( nrwork ), iwork, info )
1393 CALL claset(
'F', m, m, czero, czero, u, ldu )
1395 CALL claset(
'F', m-n, m-n, czero, cone,
1396 $ u( n+1, n+1 ), ldu )
1405 CALL clacp2(
'F', n, n, rwork( iru ), n, u, ldu )
1406 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
1407 $ work( itauq ), u, ldu, work( nwork ),
1408 $ lwork-nwork+1, ierr )
1416 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1417 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1418 $ work( itaup ), vt, ldvt, work( nwork ),
1419 $ lwork-nwork+1, ierr )
1430 IF( n.GE.mnthr1 )
THEN
1445 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1446 $ lwork-nwork+1, ierr )
1450 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
1462 CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1463 $ work( itaup ), work( nwork ), lwork-nwork+1,
1471 CALL sbdsdc(
'U',
'N', m, s, rwork( ie ), dum,1,dum,1,
1472 $ dum, idum, rwork( nrwork ), iwork, info )
1474 ELSE IF( wntqo )
THEN
1486 IF( lwork .GE. m*n + m*m + 3*m )
THEN
1497 chunk = ( lwork - m*m - 3*m ) / m
1499 itau = il + ldwrkl*chunk
1507 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1508 $ lwork-nwork+1, ierr )
1512 CALL clacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1513 CALL claset(
'U', m-1, m-1, czero, czero,
1514 $ work( il+ldwrkl ), ldwrkl )
1521 CALL cunglq( m, n, m, a, lda, work( itau ),
1522 $ work( nwork ), lwork-nwork+1, ierr )
1533 CALL cgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1534 $ work( itauq ), work( itaup ), work( nwork ),
1535 $ lwork-nwork+1, ierr )
1546 CALL sbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1547 $ m, rwork( irvt ), m, dum, idum,
1548 $ rwork( nrwork ), iwork, info )
1556 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1557 CALL cunmbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1558 $ work( itauq ), u, ldu, work( nwork ),
1559 $ lwork-nwork+1, ierr )
1567 CALL clacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1569 CALL cunmbr(
'P',
'R',
'C', m, m, m, work( il ), ldwrkl,
1570 $ work( itaup ), work( ivt ), ldwkvt,
1571 $ work( nwork ), lwork-nwork+1, ierr )
1579 DO 40 i = 1, n, chunk
1580 blk = min( n-i+1, chunk )
1581 CALL cgemm(
'N',
'N', m, blk, m, cone, work( ivt ), m,
1582 $ a( 1, i ), lda, czero, work( il ),
1584 CALL clacpy(
'F', m, blk, work( il ), ldwrkl,
1588 ELSE IF( wntqs )
THEN
1599 itau = il + ldwrkl*m
1607 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1608 $ lwork-nwork+1, ierr )
1612 CALL clacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1613 CALL claset(
'U', m-1, m-1, czero, czero,
1614 $ work( il+ldwrkl ), ldwrkl )
1621 CALL cunglq( m, n, m, a, lda, work( itau ),
1622 $ work( nwork ), lwork-nwork+1, ierr )
1633 CALL cgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1634 $ work( itauq ), work( itaup ), work( nwork ),
1635 $ lwork-nwork+1, ierr )
1646 CALL sbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1647 $ m, rwork( irvt ), m, dum, idum,
1648 $ rwork( nrwork ), iwork, info )
1656 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1657 CALL cunmbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1658 $ work( itauq ), u, ldu, work( nwork ),
1659 $ lwork-nwork+1, ierr )
1667 CALL clacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
1668 CALL cunmbr(
'P',
'R',
'C', m, m, m, work( il ), ldwrkl,
1669 $ work( itaup ), vt, ldvt, work( nwork ),
1670 $ lwork-nwork+1, ierr )
1677 CALL clacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1678 CALL cgemm(
'N',
'N', m, n, m, cone, work( il ), ldwrkl,
1679 $ a, lda, czero, vt, ldvt )
1681 ELSE IF( wntqa )
THEN
1692 itau = ivt + ldwkvt*m
1700 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1701 $ lwork-nwork+1, ierr )
1702 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
1709 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
1710 $ work( nwork ), lwork-nwork+1, ierr )
1714 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
1726 CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1727 $ work( itaup ), work( nwork ), lwork-nwork+1,
1739 CALL sbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1740 $ m, rwork( irvt ), m, dum, idum,
1741 $ rwork( nrwork ), iwork, info )
1749 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1750 CALL cunmbr(
'Q',
'L',
'N', m, m, m, a, lda,
1751 $ work( itauq ), u, ldu, work( nwork ),
1752 $ lwork-nwork+1, ierr )
1760 CALL clacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1762 CALL cunmbr(
'P',
'R',
'C', m, m, m, a, lda,
1763 $ work( itaup ), work( ivt ), ldwkvt,
1764 $ work( nwork ), lwork-nwork+1, ierr )
1771 CALL cgemm(
'N',
'N', m, n, m, cone, work( ivt ), ldwkvt,
1772 $ vt, ldvt, czero, a, lda )
1776 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
1780 ELSE IF( n.GE.mnthr2 )
THEN
1799 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1800 $ work( itaup ), work( nwork ), lwork-nwork+1,
1810 CALL sbdsdc(
'L',
'N', m, s, rwork( ie ), dum,1,dum,1,
1811 $ dum, idum, rwork( nrwork ), iwork, info )
1812 ELSE IF( wntqo )
THEN
1824 CALL clacpy(
'L', m, m, a, lda, u, ldu )
1825 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1826 $ work( nwork ), lwork-nwork+1, ierr )
1833 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
1834 $ work( nwork ), lwork-nwork+1, ierr )
1837 IF( lwork .GE. m*n + 3*m )
THEN
1841 nwork = ivt + ldwkvt*n
1847 chunk = ( lwork - 3*m ) / m
1848 nwork = ivt + ldwkvt*chunk
1857 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1858 $ m, rwork( irvt ), m, dum, idum,
1859 $ rwork( nrwork ), iwork, info )
1866 CALL clacrm( m, m, u, ldu, rwork( iru ), m, work( ivt ),
1867 $ ldwkvt, rwork( nrwork ) )
1868 CALL clacpy(
'F', m, m, work( ivt ), ldwkvt, u, ldu )
1878 DO 50 i = 1, n, chunk
1879 blk = min( n-i+1, chunk )
1880 CALL clarcm( m, blk, rwork( irvt ), m, a( 1, i ), lda,
1881 $ work( ivt ), ldwkvt, rwork( nrwork ) )
1882 CALL clacpy(
'F', m, blk, work( ivt ), ldwkvt,
1885 ELSE IF( wntqs )
THEN
1893 CALL clacpy(
'L', m, m, a, lda, u, ldu )
1894 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1895 $ work( nwork ), lwork-nwork+1, ierr )
1902 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
1903 CALL cungbr(
'P', m, n, m, vt, ldvt, work( itaup ),
1904 $ work( nwork ), lwork-nwork+1, ierr )
1915 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1916 $ m, rwork( irvt ), m, dum, idum,
1917 $ rwork( nrwork ), iwork, info )
1924 CALL clacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1926 CALL clacpy(
'F', m, m, a, lda, u, ldu )
1934 CALL clarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1936 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
1945 CALL clacpy(
'L', m, m, a, lda, u, ldu )
1946 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1947 $ work( nwork ), lwork-nwork+1, ierr )
1954 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
1955 CALL cungbr(
'P', n, n, m, vt, ldvt, work( itaup ),
1956 $ work( nwork ), lwork-nwork+1, ierr )
1967 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1968 $ m, rwork( irvt ), m, dum, idum,
1969 $ rwork( nrwork ), iwork, info )
1976 CALL clacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1978 CALL clacpy(
'F', m, m, a, lda, u, ldu )
1986 CALL clarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1988 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
2010 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2011 $ work( itaup ), work( nwork ), lwork-nwork+1,
2020 CALL sbdsdc(
'L',
'N', m, s, rwork( ie ), dum,1,dum,1,
2021 $ dum, idum, rwork( nrwork ), iwork, info )
2022 ELSE IF( wntqo )
THEN
2026 IF( lwork .GE. m*n + 3*m )
THEN
2030 CALL claset(
'F', m, n, czero, czero, work( ivt ),
2032 nwork = ivt + ldwkvt*n
2037 chunk = ( lwork - 3*m ) / m
2038 nwork = ivt + ldwkvt*chunk
2050 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
2051 $ m, rwork( irvt ), m, dum, idum,
2052 $ rwork( nrwork ), iwork, info )
2060 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
2061 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
2062 $ work( itauq ), u, ldu, work( nwork ),
2063 $ lwork-nwork+1, ierr )
2065 IF( lwork .GE. m*n + 3*m )
THEN
2075 CALL clacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
2077 CALL cunmbr(
'P',
'R',
'C', m, n, m, a, lda,
2078 $ work( itaup ), work( ivt ), ldwkvt,
2079 $ work( nwork ), lwork-nwork+1, ierr )
2080 CALL clacpy(
'F', m, n, work( ivt ), ldwkvt, a, lda )
2089 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
2090 $ work( nwork ), lwork-nwork+1, ierr )
2100 DO 60 i = 1, n, chunk
2101 blk = min( n-i+1, chunk )
2102 CALL clarcm( m, blk, rwork( irvt ), m, a( 1, i ),
2103 $ lda, work( ivt ), ldwkvt,
2105 CALL clacpy(
'F', m, blk, work( ivt ), ldwkvt,
2109 ELSE IF( wntqs )
THEN
2121 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
2122 $ m, rwork( irvt ), m, dum, idum,
2123 $ rwork( nrwork ), iwork, info )
2131 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
2132 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
2133 $ work( itauq ), u, ldu, work( nwork ),
2134 $ lwork-nwork+1, ierr )
2142 CALL claset(
'F', m, n, czero, czero, vt, ldvt )
2143 CALL clacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
2144 CALL cunmbr(
'P',
'R',
'C', m, n, m, a, lda,
2145 $ work( itaup ), vt, ldvt, work( nwork ),
2146 $ lwork-nwork+1, ierr )
2160 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
2161 $ m, rwork( irvt ), m, dum, idum,
2162 $ rwork( nrwork ), iwork, info )
2170 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
2171 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
2172 $ work( itauq ), u, ldu, work( nwork ),
2173 $ lwork-nwork+1, ierr )
2177 CALL claset(
'F', n, n, czero, cone, vt, ldvt )
2185 CALL clacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
2186 CALL cunmbr(
'P',
'R',
'C', n, n, m, a, lda,
2187 $ work( itaup ), vt, ldvt, work( nwork ),
2188 $ lwork-nwork+1, ierr )
2197 IF( iscl.EQ.1 )
THEN
2198 IF( anrm.GT.bignum )
2199 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
2201 IF( info.NE.0 .AND. anrm.GT.bignum )
2202 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
2203 $ rwork( ie ), minmn, ierr )
2204 IF( anrm.LT.smlnum )
2205 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
2207 IF( info.NE.0 .AND. anrm.LT.smlnum )
2208 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
2209 $ rwork( ie ), minmn, ierr )
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
subroutine clacp2(UPLO, M, N, A, LDA, B, LDB)
CLACP2 copies all or part of a real two-dimensional array to a complex array.
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
subroutine clarcm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLARCM copies all or part of a real two-dimensional array to a complex array.
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
CGESDD
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM