222 SUBROUTINE cgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
223 $ work, lwork, rwork, iwork, info )
232 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
236 REAL RWORK( * ), S( * )
237 COMPLEX A( lda, * ), U( ldu, * ), VT( ldvt, * ),
245 parameter( lquerv = -1 )
247 parameter( czero = ( 0.0e+0, 0.0e+0 ),
248 $ cone = ( 1.0e+0, 0.0e+0 ) )
250 parameter( zero = 0.0e+0, one = 1.0e+0 )
253 LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
254 INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
255 $ iscl, itau, itaup, itauq, iu, ivt, ldwkvt,
256 $ ldwrkl, ldwrkr, ldwrku, maxwrk, minmn, minwrk,
257 $ mnthr1, mnthr2, nrwork, nwork, wrkbl
258 REAL ANRM, BIGNUM, EPS, SMLNUM
273 EXTERNAL clange, slamch, ilaenv, lsame
276 INTRINSIC int, max, min, sqrt
284 mnthr1 = int( minmn*17.0 / 9.0 )
285 mnthr2 = int( minmn*5.0 / 3.0 )
286 wntqa = lsame( jobz,
'A' )
287 wntqs = lsame( jobz,
'S' )
288 wntqas = wntqa .OR. wntqs
289 wntqo = lsame( jobz,
'O' )
290 wntqn = lsame( jobz,
'N' )
294 IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) )
THEN
296 ELSE IF( m.LT.0 )
THEN
298 ELSE IF( n.LT.0 )
THEN
300 ELSE IF( lda.LT.max( 1, m ) )
THEN
302 ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
303 $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) )
THEN
305 ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
306 $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
307 $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) )
THEN
319 IF( info.EQ.0 .AND. m.GT.0 .AND. n.GT.0 )
THEN
329 IF( m.GE.mnthr1 )
THEN
334 maxwrk = n + n*ilaenv( 1,
'CGEQRF',
' ', m, n, -1,
336 maxwrk = max( maxwrk, 2*n+2*n*
337 $ ilaenv( 1,
'CGEBRD',
' ', n, n, -1, -1 ) )
339 ELSE IF( wntqo )
THEN
343 wrkbl = n + n*ilaenv( 1,
'CGEQRF',
' ', m, n, -1, -1 )
344 wrkbl = max( wrkbl, n+n*ilaenv( 1,
'CUNGQR',
' ', m,
346 wrkbl = max( wrkbl, 2*n+2*n*
347 $ ilaenv( 1,
'CGEBRD',
' ', n, n, -1, -1 ) )
348 wrkbl = max( wrkbl, 2*n+n*
349 $ ilaenv( 1,
'CUNMBR',
'QLN', n, n, n, -1 ) )
350 wrkbl = max( wrkbl, 2*n+n*
351 $ ilaenv( 1,
'CUNMBR',
'PRC', n, n, n, -1 ) )
352 maxwrk = m*n + n*n + wrkbl
354 ELSE IF( wntqs )
THEN
358 wrkbl = n + n*ilaenv( 1,
'CGEQRF',
' ', m, n, -1, -1 )
359 wrkbl = max( wrkbl, n+n*ilaenv( 1,
'CUNGQR',
' ', m,
361 wrkbl = max( wrkbl, 2*n+2*n*
362 $ ilaenv( 1,
'CGEBRD',
' ', n, n, -1, -1 ) )
363 wrkbl = max( wrkbl, 2*n+n*
364 $ ilaenv( 1,
'CUNMBR',
'QLN', n, n, n, -1 ) )
365 wrkbl = max( wrkbl, 2*n+n*
366 $ ilaenv( 1,
'CUNMBR',
'PRC', n, n, n, -1 ) )
369 ELSE IF( wntqa )
THEN
373 wrkbl = n + n*ilaenv( 1,
'CGEQRF',
' ', m, n, -1, -1 )
374 wrkbl = max( wrkbl, n+m*ilaenv( 1,
'CUNGQR',
' ', m,
376 wrkbl = max( wrkbl, 2*n+2*n*
377 $ ilaenv( 1,
'CGEBRD',
' ', n, n, -1, -1 ) )
378 wrkbl = max( wrkbl, 2*n+n*
379 $ ilaenv( 1,
'CUNMBR',
'QLN', n, n, n, -1 ) )
380 wrkbl = max( wrkbl, 2*n+n*
381 $ ilaenv( 1,
'CUNMBR',
'PRC', n, n, n, -1 ) )
383 minwrk = n*n + 2*n + m
385 ELSE IF( m.GE.mnthr2 )
THEN
389 maxwrk = 2*n + ( m+n )*ilaenv( 1,
'CGEBRD',
' ', m, n,
393 maxwrk = max( maxwrk, 2*n+n*
394 $ ilaenv( 1,
'CUNGBR',
'P', n, n, n, -1 ) )
395 maxwrk = max( maxwrk, 2*n+n*
396 $ ilaenv( 1,
'CUNGBR',
'Q', m, n, n, -1 ) )
397 maxwrk = maxwrk + m*n
398 minwrk = minwrk + n*n
399 ELSE IF( wntqs )
THEN
400 maxwrk = max( maxwrk, 2*n+n*
401 $ ilaenv( 1,
'CUNGBR',
'P', n, n, n, -1 ) )
402 maxwrk = max( maxwrk, 2*n+n*
403 $ ilaenv( 1,
'CUNGBR',
'Q', m, n, n, -1 ) )
404 ELSE IF( wntqa )
THEN
405 maxwrk = max( maxwrk, 2*n+n*
406 $ ilaenv( 1,
'CUNGBR',
'P', n, n, n, -1 ) )
407 maxwrk = max( maxwrk, 2*n+m*
408 $ ilaenv( 1,
'CUNGBR',
'Q', m, m, n, -1 ) )
414 maxwrk = 2*n + ( m+n )*ilaenv( 1,
'CGEBRD',
' ', m, n,
418 maxwrk = max( maxwrk, 2*n+n*
419 $ ilaenv( 1,
'CUNMBR',
'PRC', n, n, n, -1 ) )
420 maxwrk = max( maxwrk, 2*n+n*
421 $ ilaenv( 1,
'CUNMBR',
'QLN', m, n, n, -1 ) )
422 maxwrk = maxwrk + m*n
423 minwrk = minwrk + n*n
424 ELSE IF( wntqs )
THEN
425 maxwrk = max( maxwrk, 2*n+n*
426 $ ilaenv( 1,
'CUNMBR',
'PRC', n, n, n, -1 ) )
427 maxwrk = max( maxwrk, 2*n+n*
428 $ ilaenv( 1,
'CUNMBR',
'QLN', m, n, n, -1 ) )
429 ELSE IF( wntqa )
THEN
430 maxwrk = max( maxwrk, 2*n+n*
431 $ ilaenv( 1,
'CUNGBR',
'PRC', n, n, n, -1 ) )
432 maxwrk = max( maxwrk, 2*n+m*
433 $ ilaenv( 1,
'CUNGBR',
'QLN', m, m, n, -1 ) )
445 IF( n.GE.mnthr1 )
THEN
450 maxwrk = m + m*ilaenv( 1,
'CGELQF',
' ', m, n, -1,
452 maxwrk = max( maxwrk, 2*m+2*m*
453 $ ilaenv( 1,
'CGEBRD',
' ', m, m, -1, -1 ) )
455 ELSE IF( wntqo )
THEN
459 wrkbl = m + m*ilaenv( 1,
'CGELQF',
' ', m, n, -1, -1 )
460 wrkbl = max( wrkbl, m+m*ilaenv( 1,
'CUNGLQ',
' ', m,
462 wrkbl = max( wrkbl, 2*m+2*m*
463 $ ilaenv( 1,
'CGEBRD',
' ', m, m, -1, -1 ) )
464 wrkbl = max( wrkbl, 2*m+m*
465 $ ilaenv( 1,
'CUNMBR',
'PRC', m, m, m, -1 ) )
466 wrkbl = max( wrkbl, 2*m+m*
467 $ ilaenv( 1,
'CUNMBR',
'QLN', m, m, m, -1 ) )
468 maxwrk = m*n + m*m + wrkbl
470 ELSE IF( wntqs )
THEN
474 wrkbl = m + m*ilaenv( 1,
'CGELQF',
' ', m, n, -1, -1 )
475 wrkbl = max( wrkbl, m+m*ilaenv( 1,
'CUNGLQ',
' ', m,
477 wrkbl = max( wrkbl, 2*m+2*m*
478 $ ilaenv( 1,
'CGEBRD',
' ', m, m, -1, -1 ) )
479 wrkbl = max( wrkbl, 2*m+m*
480 $ ilaenv( 1,
'CUNMBR',
'PRC', m, m, m, -1 ) )
481 wrkbl = max( wrkbl, 2*m+m*
482 $ ilaenv( 1,
'CUNMBR',
'QLN', m, m, m, -1 ) )
485 ELSE IF( wntqa )
THEN
489 wrkbl = m + m*ilaenv( 1,
'CGELQF',
' ', m, n, -1, -1 )
490 wrkbl = max( wrkbl, m+n*ilaenv( 1,
'CUNGLQ',
' ', n,
492 wrkbl = max( wrkbl, 2*m+2*m*
493 $ ilaenv( 1,
'CGEBRD',
' ', m, m, -1, -1 ) )
494 wrkbl = max( wrkbl, 2*m+m*
495 $ ilaenv( 1,
'CUNMBR',
'PRC', m, m, m, -1 ) )
496 wrkbl = max( wrkbl, 2*m+m*
497 $ ilaenv( 1,
'CUNMBR',
'QLN', m, m, m, -1 ) )
499 minwrk = m*m + 2*m + n
501 ELSE IF( n.GE.mnthr2 )
THEN
505 maxwrk = 2*m + ( m+n )*ilaenv( 1,
'CGEBRD',
' ', m, n,
509 maxwrk = max( maxwrk, 2*m+m*
510 $ ilaenv( 1,
'CUNGBR',
'P', m, n, m, -1 ) )
511 maxwrk = max( maxwrk, 2*m+m*
512 $ ilaenv( 1,
'CUNGBR',
'Q', m, m, n, -1 ) )
513 maxwrk = maxwrk + m*n
514 minwrk = minwrk + m*m
515 ELSE IF( wntqs )
THEN
516 maxwrk = max( maxwrk, 2*m+m*
517 $ ilaenv( 1,
'CUNGBR',
'P', m, n, m, -1 ) )
518 maxwrk = max( maxwrk, 2*m+m*
519 $ ilaenv( 1,
'CUNGBR',
'Q', m, m, n, -1 ) )
520 ELSE IF( wntqa )
THEN
521 maxwrk = max( maxwrk, 2*m+n*
522 $ ilaenv( 1,
'CUNGBR',
'P', n, n, m, -1 ) )
523 maxwrk = max( maxwrk, 2*m+m*
524 $ ilaenv( 1,
'CUNGBR',
'Q', m, m, n, -1 ) )
530 maxwrk = 2*m + ( m+n )*ilaenv( 1,
'CGEBRD',
' ', m, n,
534 maxwrk = max( maxwrk, 2*m+m*
535 $ ilaenv( 1,
'CUNMBR',
'PRC', m, n, m, -1 ) )
536 maxwrk = max( maxwrk, 2*m+m*
537 $ ilaenv( 1,
'CUNMBR',
'QLN', m, m, n, -1 ) )
538 maxwrk = maxwrk + m*n
539 minwrk = minwrk + m*m
540 ELSE IF( wntqs )
THEN
541 maxwrk = max( maxwrk, 2*m+m*
542 $ ilaenv( 1,
'CUNGBR',
'PRC', m, n, m, -1 ) )
543 maxwrk = max( maxwrk, 2*m+m*
544 $ ilaenv( 1,
'CUNGBR',
'QLN', m, m, n, -1 ) )
545 ELSE IF( wntqa )
THEN
546 maxwrk = max( maxwrk, 2*m+n*
547 $ ilaenv( 1,
'CUNGBR',
'PRC', n, n, m, -1 ) )
548 maxwrk = max( maxwrk, 2*m+m*
549 $ ilaenv( 1,
'CUNGBR',
'QLN', m, m, n, -1 ) )
553 maxwrk = max( maxwrk, minwrk )
557 IF( lwork.LT.minwrk .AND. lwork.NE.lquerv )
564 CALL xerbla(
'CGESDD', -info )
567 IF( lwork.EQ.lquerv )
569 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
576 smlnum = sqrt( slamch(
'S' ) ) / eps
577 bignum = one / smlnum
581 anrm = clange(
'M', m, n, a, lda, dum )
583 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
585 CALL clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
586 ELSE IF( anrm.GT.bignum )
THEN
588 CALL clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
597 IF( m.GE.mnthr1 )
THEN
611 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
612 $ lwork-nwork+1, ierr )
616 CALL claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
627 CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
628 $ work( itaup ), work( nwork ), lwork-nwork+1,
636 CALL sbdsdc(
'U',
'N', n, s, rwork( ie ), dum, 1, dum, 1,
637 $ dum, idum, rwork( nrwork ), iwork, info )
639 ELSE IF( wntqo )
THEN
651 IF( lwork.GE.m*n+n*n+3*n )
THEN
657 ldwrkr = ( lwork-n*n-3*n ) / n
666 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
667 $ lwork-nwork+1, ierr )
671 CALL clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
672 CALL claset(
'L', n-1, n-1, czero, czero, work( ir+1 ),
679 CALL cungqr( m, n, n, a, lda, work( itau ),
680 $ work( nwork ), lwork-nwork+1, ierr )
690 CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
691 $ work( itauq ), work( itaup ), work( nwork ),
692 $ lwork-nwork+1, ierr )
703 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
704 $ n, rwork( irvt ), n, dum, idum,
705 $ rwork( nrwork ), iwork, info )
712 CALL clacp2(
'F', n, n, rwork( iru ), n, work( iu ),
714 CALL cunmbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
715 $ work( itauq ), work( iu ), ldwrku,
716 $ work( nwork ), lwork-nwork+1, ierr )
723 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
724 CALL cunmbr(
'P',
'R',
'C', n, n, n, work( ir ), ldwrkr,
725 $ work( itaup ), vt, ldvt, work( nwork ),
726 $ lwork-nwork+1, ierr )
733 DO 10 i = 1, m, ldwrkr
734 chunk = min( m-i+1, ldwrkr )
735 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
736 $ lda, work( iu ), ldwrku, czero,
737 $ work( ir ), ldwrkr )
738 CALL clacpy(
'F', chunk, n, work( ir ), ldwrkr,
742 ELSE IF( wntqs )
THEN
760 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
761 $ lwork-nwork+1, ierr )
765 CALL clacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
766 CALL claset(
'L', n-1, n-1, czero, czero, work( ir+1 ),
773 CALL cungqr( m, n, n, a, lda, work( itau ),
774 $ work( nwork ), lwork-nwork+1, ierr )
784 CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
785 $ work( itauq ), work( itaup ), work( nwork ),
786 $ lwork-nwork+1, ierr )
797 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
798 $ n, rwork( irvt ), n, dum, idum,
799 $ rwork( nrwork ), iwork, info )
806 CALL clacp2(
'F', n, n, rwork( iru ), n, u, ldu )
807 CALL cunmbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
808 $ work( itauq ), u, ldu, work( nwork ),
809 $ lwork-nwork+1, ierr )
816 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
817 CALL cunmbr(
'P',
'R',
'C', n, n, n, work( ir ), ldwrkr,
818 $ work( itaup ), vt, ldvt, work( nwork ),
819 $ lwork-nwork+1, ierr )
826 CALL clacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
827 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda, work( ir ),
828 $ ldwrkr, czero, u, ldu )
830 ELSE IF( wntqa )
THEN
848 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
849 $ lwork-nwork+1, ierr )
850 CALL clacpy(
'L', m, n, a, lda, u, ldu )
856 CALL cungqr( m, m, n, u, ldu, work( itau ),
857 $ work( nwork ), lwork-nwork+1, ierr )
861 CALL claset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
872 CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
873 $ work( itaup ), work( nwork ), lwork-nwork+1,
885 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
886 $ n, rwork( irvt ), n, dum, idum,
887 $ rwork( nrwork ), iwork, info )
894 CALL clacp2(
'F', n, n, rwork( iru ), n, work( iu ),
896 CALL cunmbr(
'Q',
'L',
'N', n, n, n, a, lda,
897 $ work( itauq ), work( iu ), ldwrku,
898 $ work( nwork ), lwork-nwork+1, ierr )
905 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
906 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
907 $ work( itaup ), vt, ldvt, work( nwork ),
908 $ lwork-nwork+1, ierr )
915 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu, work( iu ),
916 $ ldwrku, czero, a, lda )
920 CALL clacpy(
'F', m, n, a, lda, u, ldu )
924 ELSE IF( m.GE.mnthr2 )
THEN
942 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
943 $ work( itaup ), work( nwork ), lwork-nwork+1,
951 CALL sbdsdc(
'U',
'N', n, s, rwork( ie ), dum, 1, dum, 1,
952 $ dum, idum, rwork( nrwork ), iwork, info )
953 ELSE IF( wntqo )
THEN
963 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
964 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
965 $ work( nwork ), lwork-nwork+1, ierr )
971 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
972 $ work( nwork ), lwork-nwork+1, ierr )
974 IF( lwork.GE.m*n+3*n )
THEN
983 ldwrku = ( lwork-3*n ) / n
985 nwork = iu + ldwrku*n
993 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
994 $ n, rwork( irvt ), n, dum, idum,
995 $ rwork( nrwork ), iwork, info )
1002 CALL clarcm( n, n, rwork( irvt ), n, vt, ldvt,
1003 $ work( iu ), ldwrku, rwork( nrwork ) )
1004 CALL clacpy(
'F', n, n, work( iu ), ldwrku, vt, ldvt )
1012 DO 20 i = 1, m, ldwrku
1013 chunk = min( m-i+1, ldwrku )
1014 CALL clacrm( chunk, n, a( i, 1 ), lda, rwork( iru ),
1015 $ n, work( iu ), ldwrku, rwork( nrwork ) )
1016 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
1020 ELSE IF( wntqs )
THEN
1026 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1027 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1028 $ work( nwork ), lwork-nwork+1, ierr )
1034 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1035 CALL cungbr(
'Q', m, n, n, u, ldu, work( itauq ),
1036 $ work( nwork ), lwork-nwork+1, ierr )
1047 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1048 $ n, rwork( irvt ), n, dum, idum,
1049 $ rwork( nrwork ), iwork, info )
1056 CALL clarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1058 CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
1066 CALL clacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1068 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1075 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1076 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1077 $ work( nwork ), lwork-nwork+1, ierr )
1083 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1084 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1085 $ work( nwork ), lwork-nwork+1, ierr )
1096 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1097 $ n, rwork( irvt ), n, dum, idum,
1098 $ rwork( nrwork ), iwork, info )
1105 CALL clarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1107 CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
1115 CALL clacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1117 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1138 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1139 $ work( itaup ), work( nwork ), lwork-nwork+1,
1147 CALL sbdsdc(
'U',
'N', n, s, rwork( ie ), dum, 1, dum, 1,
1148 $ dum, idum, rwork( nrwork ), iwork, info )
1149 ELSE IF( wntqo )
THEN
1154 IF( lwork.GE.m*n+3*n )
THEN
1163 ldwrku = ( lwork-3*n ) / n
1165 nwork = iu + ldwrku*n
1173 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1174 $ n, rwork( irvt ), n, dum, idum,
1175 $ rwork( nrwork ), iwork, info )
1182 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1183 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1184 $ work( itaup ), vt, ldvt, work( nwork ),
1185 $ lwork-nwork+1, ierr )
1187 IF( lwork.GE.m*n+3*n )
THEN
1195 CALL claset(
'F', m, n, czero, czero, work( iu ),
1197 CALL clacp2(
'F', n, n, rwork( iru ), n, work( iu ),
1199 CALL cunmbr(
'Q',
'L',
'N', m, n, n, a, lda,
1200 $ work( itauq ), work( iu ), ldwrku,
1201 $ work( nwork ), lwork-nwork+1, ierr )
1202 CALL clacpy(
'F', m, n, work( iu ), ldwrku, a, lda )
1209 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
1210 $ work( nwork ), lwork-nwork+1, ierr )
1218 DO 30 i = 1, m, ldwrku
1219 chunk = min( m-i+1, ldwrku )
1220 CALL clacrm( chunk, n, a( i, 1 ), lda,
1221 $ rwork( iru ), n, work( iu ), ldwrku,
1223 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
1228 ELSE IF( wntqs )
THEN
1239 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1240 $ n, rwork( irvt ), n, dum, idum,
1241 $ rwork( nrwork ), iwork, info )
1248 CALL claset(
'F', m, n, czero, czero, u, ldu )
1249 CALL clacp2(
'F', n, n, rwork( iru ), n, u, ldu )
1250 CALL cunmbr(
'Q',
'L',
'N', m, n, n, a, lda,
1251 $ work( itauq ), u, ldu, work( nwork ),
1252 $ lwork-nwork+1, ierr )
1259 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1260 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1261 $ work( itaup ), vt, ldvt, work( nwork ),
1262 $ lwork-nwork+1, ierr )
1274 CALL sbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1275 $ n, rwork( irvt ), n, dum, idum,
1276 $ rwork( nrwork ), iwork, info )
1280 CALL claset(
'F', m, m, czero, czero, u, ldu )
1282 CALL claset(
'F', m-n, m-n, czero, cone,
1283 $ u( n+1, n+1 ), ldu )
1291 CALL clacp2(
'F', n, n, rwork( iru ), n, u, ldu )
1292 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
1293 $ work( itauq ), u, ldu, work( nwork ),
1294 $ lwork-nwork+1, ierr )
1301 CALL clacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1302 CALL cunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1303 $ work( itaup ), vt, ldvt, work( nwork ),
1304 $ lwork-nwork+1, ierr )
1315 IF( n.GE.mnthr1 )
THEN
1329 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1330 $ lwork-nwork+1, ierr )
1334 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
1345 CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1346 $ work( itaup ), work( nwork ), lwork-nwork+1,
1354 CALL sbdsdc(
'U',
'N', m, s, rwork( ie ), dum, 1, dum, 1,
1355 $ dum, idum, rwork( nrwork ), iwork, info )
1357 ELSE IF( wntqo )
THEN
1369 IF( lwork.GE.m*n+m*m+3*m )
THEN
1380 chunk = ( lwork-m*m-3*m ) / m
1382 itau = il + ldwrkl*chunk
1389 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1390 $ lwork-nwork+1, ierr )
1394 CALL clacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1395 CALL claset(
'U', m-1, m-1, czero, czero,
1396 $ work( il+ldwrkl ), ldwrkl )
1402 CALL cunglq( m, n, m, a, lda, work( itau ),
1403 $ work( nwork ), lwork-nwork+1, ierr )
1413 CALL cgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1414 $ work( itauq ), work( itaup ), work( nwork ),
1415 $ lwork-nwork+1, ierr )
1426 CALL sbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1427 $ m, rwork( irvt ), m, dum, idum,
1428 $ rwork( nrwork ), iwork, info )
1435 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1436 CALL cunmbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1437 $ work( itauq ), u, ldu, work( nwork ),
1438 $ lwork-nwork+1, ierr )
1445 CALL clacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1447 CALL cunmbr(
'P',
'R',
'C', m, m, m, work( il ), ldwrkl,
1448 $ work( itaup ), work( ivt ), ldwkvt,
1449 $ work( nwork ), lwork-nwork+1, ierr )
1456 DO 40 i = 1, n, chunk
1457 blk = min( n-i+1, chunk )
1458 CALL cgemm(
'N',
'N', m, blk, m, cone, work( ivt ), m,
1459 $ a( 1, i ), lda, czero, work( il ),
1461 CALL clacpy(
'F', m, blk, work( il ), ldwrkl,
1465 ELSE IF( wntqs )
THEN
1476 itau = il + ldwrkl*m
1483 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1484 $ lwork-nwork+1, ierr )
1488 CALL clacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1489 CALL claset(
'U', m-1, m-1, czero, czero,
1490 $ work( il+ldwrkl ), ldwrkl )
1496 CALL cunglq( m, n, m, a, lda, work( itau ),
1497 $ work( nwork ), lwork-nwork+1, ierr )
1507 CALL cgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1508 $ work( itauq ), work( itaup ), work( nwork ),
1509 $ lwork-nwork+1, ierr )
1520 CALL sbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1521 $ m, rwork( irvt ), m, dum, idum,
1522 $ rwork( nrwork ), iwork, info )
1529 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1530 CALL cunmbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1531 $ work( itauq ), u, ldu, work( nwork ),
1532 $ lwork-nwork+1, ierr )
1539 CALL clacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
1540 CALL cunmbr(
'P',
'R',
'C', m, m, m, work( il ), ldwrkl,
1541 $ work( itaup ), vt, ldvt, work( nwork ),
1542 $ lwork-nwork+1, ierr )
1549 CALL clacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1550 CALL cgemm(
'N',
'N', m, n, m, cone, work( il ), ldwrkl,
1551 $ a, lda, czero, vt, ldvt )
1553 ELSE IF( wntqa )
THEN
1564 itau = ivt + ldwkvt*m
1571 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
1572 $ lwork-nwork+1, ierr )
1573 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
1579 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
1580 $ work( nwork ), lwork-nwork+1, ierr )
1584 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
1595 CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1596 $ work( itaup ), work( nwork ), lwork-nwork+1,
1608 CALL sbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1609 $ m, rwork( irvt ), m, dum, idum,
1610 $ rwork( nrwork ), iwork, info )
1617 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1618 CALL cunmbr(
'Q',
'L',
'N', m, m, m, a, lda,
1619 $ work( itauq ), u, ldu, work( nwork ),
1620 $ lwork-nwork+1, ierr )
1627 CALL clacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1629 CALL cunmbr(
'P',
'R',
'C', m, m, m, a, lda,
1630 $ work( itaup ), work( ivt ), ldwkvt,
1631 $ work( nwork ), lwork-nwork+1, ierr )
1638 CALL cgemm(
'N',
'N', m, n, m, cone, work( ivt ),
1639 $ ldwkvt, vt, ldvt, czero, a, lda )
1643 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
1647 ELSE IF( n.GE.mnthr2 )
THEN
1666 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1667 $ work( itaup ), work( nwork ), lwork-nwork+1,
1676 CALL sbdsdc(
'L',
'N', m, s, rwork( ie ), dum, 1, dum, 1,
1677 $ dum, idum, rwork( nrwork ), iwork, info )
1678 ELSE IF( wntqo )
THEN
1688 CALL clacpy(
'L', m, m, a, lda, u, ldu )
1689 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1690 $ work( nwork ), lwork-nwork+1, ierr )
1696 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
1697 $ work( nwork ), lwork-nwork+1, ierr )
1700 IF( lwork.GE.m*n+3*m )
THEN
1704 nwork = ivt + ldwkvt*n
1710 chunk = ( lwork-3*m ) / m
1711 nwork = ivt + ldwkvt*chunk
1720 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1721 $ m, rwork( irvt ), m, dum, idum,
1722 $ rwork( nrwork ), iwork, info )
1729 CALL clacrm( m, m, u, ldu, rwork( iru ), m, work( ivt ),
1730 $ ldwkvt, rwork( nrwork ) )
1731 CALL clacpy(
'F', m, m, work( ivt ), ldwkvt, u, ldu )
1739 DO 50 i = 1, n, chunk
1740 blk = min( n-i+1, chunk )
1741 CALL clarcm( m, blk, rwork( irvt ), m, a( 1, i ), lda,
1742 $ work( ivt ), ldwkvt, rwork( nrwork ) )
1743 CALL clacpy(
'F', m, blk, work( ivt ), ldwkvt,
1746 ELSE IF( wntqs )
THEN
1752 CALL clacpy(
'L', m, m, a, lda, u, ldu )
1753 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1754 $ work( nwork ), lwork-nwork+1, ierr )
1760 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
1761 CALL cungbr(
'P', m, n, m, vt, ldvt, work( itaup ),
1762 $ work( nwork ), lwork-nwork+1, ierr )
1773 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1774 $ m, rwork( irvt ), m, dum, idum,
1775 $ rwork( nrwork ), iwork, info )
1782 CALL clacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1784 CALL clacpy(
'F', m, m, a, lda, u, ldu )
1792 CALL clarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1794 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
1801 CALL clacpy(
'L', m, m, a, lda, u, ldu )
1802 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1803 $ work( nwork ), lwork-nwork+1, ierr )
1809 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
1810 CALL cungbr(
'P', n, n, m, vt, ldvt, work( itaup ),
1811 $ work( nwork ), lwork-nwork+1, ierr )
1822 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1823 $ m, rwork( irvt ), m, dum, idum,
1824 $ rwork( nrwork ), iwork, info )
1831 CALL clacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1833 CALL clacpy(
'F', m, m, a, lda, u, ldu )
1840 CALL clarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1842 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
1863 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1864 $ work( itaup ), work( nwork ), lwork-nwork+1,
1872 CALL sbdsdc(
'L',
'N', m, s, rwork( ie ), dum, 1, dum, 1,
1873 $ dum, idum, rwork( nrwork ), iwork, info )
1874 ELSE IF( wntqo )
THEN
1877 IF( lwork.GE.m*n+3*m )
THEN
1881 CALL claset(
'F', m, n, czero, czero, work( ivt ),
1883 nwork = ivt + ldwkvt*n
1888 chunk = ( lwork-3*m ) / m
1889 nwork = ivt + ldwkvt*chunk
1901 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1902 $ m, rwork( irvt ), m, dum, idum,
1903 $ rwork( nrwork ), iwork, info )
1910 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1911 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
1912 $ work( itauq ), u, ldu, work( nwork ),
1913 $ lwork-nwork+1, ierr )
1915 IF( lwork.GE.m*n+3*m )
THEN
1923 CALL clacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1925 CALL cunmbr(
'P',
'R',
'C', m, n, m, a, lda,
1926 $ work( itaup ), work( ivt ), ldwkvt,
1927 $ work( nwork ), lwork-nwork+1, ierr )
1928 CALL clacpy(
'F', m, n, work( ivt ), ldwkvt, a, lda )
1935 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
1936 $ work( nwork ), lwork-nwork+1, ierr )
1944 DO 60 i = 1, n, chunk
1945 blk = min( n-i+1, chunk )
1946 CALL clarcm( m, blk, rwork( irvt ), m, a( 1, i ),
1947 $ lda, work( ivt ), ldwkvt,
1949 CALL clacpy(
'F', m, blk, work( ivt ), ldwkvt,
1953 ELSE IF( wntqs )
THEN
1964 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1965 $ m, rwork( irvt ), m, dum, idum,
1966 $ rwork( nrwork ), iwork, info )
1973 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1974 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
1975 $ work( itauq ), u, ldu, work( nwork ),
1976 $ lwork-nwork+1, ierr )
1983 CALL claset(
'F', m, n, czero, czero, vt, ldvt )
1984 CALL clacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
1985 CALL cunmbr(
'P',
'R',
'C', m, n, m, a, lda,
1986 $ work( itaup ), vt, ldvt, work( nwork ),
1987 $ lwork-nwork+1, ierr )
2000 CALL sbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
2001 $ m, rwork( irvt ), m, dum, idum,
2002 $ rwork( nrwork ), iwork, info )
2009 CALL clacp2(
'F', m, m, rwork( iru ), m, u, ldu )
2010 CALL cunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
2011 $ work( itauq ), u, ldu, work( nwork ),
2012 $ lwork-nwork+1, ierr )
2016 CALL claset(
'F', n, n, czero, cone, vt, ldvt )
2023 CALL clacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
2024 CALL cunmbr(
'P',
'R',
'C', n, n, m, a, lda,
2025 $ work( itaup ), vt, ldvt, work( nwork ),
2026 $ lwork-nwork+1, ierr )
2035 IF( iscl.EQ.1 )
THEN
2036 IF( anrm.GT.bignum )
2037 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
2039 IF( info.NE.0 .AND. anrm.GT.bignum )
2040 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
2041 $ rwork( ie ), minmn, ierr )
2042 IF( anrm.LT.smlnum )
2043 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
2045 IF( info.NE.0 .AND. anrm.LT.smlnum )
2046 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
2047 $ rwork( ie ), minmn, ierr )
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
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 cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
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 clacp2(UPLO, M, N, A, LDA, B, LDB)
CLACP2 copies all or part of a real two-dimensional array to a complex array.
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
subroutine xerbla(SRNAME, INFO)
XERBLA
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 cgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
CGESDD
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine clacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
CLACRM multiplies a complex matrix by a square real matrix.
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
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 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
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR