222 SUBROUTINE zgesdd( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
223 $ lwork, rwork, iwork, info )
232 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
236 DOUBLE PRECISION RWORK( * ), S( * )
237 COMPLEX*16 A( lda, * ), U( ldu, * ), VT( ldvt, * ),
245 parameter( lquerv = -1 )
246 COMPLEX*16 CZERO, CONE
247 parameter( czero = ( 0.0d+0, 0.0d+0 ),
248 $ cone = ( 1.0d+0, 0.0d+0 ) )
249 DOUBLE PRECISION ZERO, ONE
250 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
262 DOUBLE PRECISION DUM( 1 )
272 DOUBLE PRECISION DLAMCH, ZLANGE
273 EXTERNAL lsame, ilaenv, dlamch, zlange
276 INTRINSIC int, max, min, sqrt
284 mnthr1 = int( minmn*17.0d0 / 9.0d0 )
285 mnthr2 = int( minmn*5.0d0 / 3.0d0 )
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,
'ZGEQRF',
' ', m, n, -1,
336 maxwrk = max( maxwrk, 2*n+2*n*
337 $ ilaenv( 1,
'ZGEBRD',
' ', n, n, -1, -1 ) )
339 ELSE IF( wntqo )
THEN
343 wrkbl = n + n*ilaenv( 1,
'ZGEQRF',
' ', m, n, -1, -1 )
344 wrkbl = max( wrkbl, n+n*ilaenv( 1,
'ZUNGQR',
' ', m,
346 wrkbl = max( wrkbl, 2*n+2*n*
347 $ ilaenv( 1,
'ZGEBRD',
' ', n, n, -1, -1 ) )
348 wrkbl = max( wrkbl, 2*n+n*
349 $ ilaenv( 1,
'ZUNMBR',
'QLN', n, n, n, -1 ) )
350 wrkbl = max( wrkbl, 2*n+n*
351 $ ilaenv( 1,
'ZUNMBR',
'PRC', n, n, n, -1 ) )
352 maxwrk = m*n + n*n + wrkbl
354 ELSE IF( wntqs )
THEN
358 wrkbl = n + n*ilaenv( 1,
'ZGEQRF',
' ', m, n, -1, -1 )
359 wrkbl = max( wrkbl, n+n*ilaenv( 1,
'ZUNGQR',
' ', m,
361 wrkbl = max( wrkbl, 2*n+2*n*
362 $ ilaenv( 1,
'ZGEBRD',
' ', n, n, -1, -1 ) )
363 wrkbl = max( wrkbl, 2*n+n*
364 $ ilaenv( 1,
'ZUNMBR',
'QLN', n, n, n, -1 ) )
365 wrkbl = max( wrkbl, 2*n+n*
366 $ ilaenv( 1,
'ZUNMBR',
'PRC', n, n, n, -1 ) )
369 ELSE IF( wntqa )
THEN
373 wrkbl = n + n*ilaenv( 1,
'ZGEQRF',
' ', m, n, -1, -1 )
374 wrkbl = max( wrkbl, n+m*ilaenv( 1,
'ZUNGQR',
' ', m,
376 wrkbl = max( wrkbl, 2*n+2*n*
377 $ ilaenv( 1,
'ZGEBRD',
' ', n, n, -1, -1 ) )
378 wrkbl = max( wrkbl, 2*n+n*
379 $ ilaenv( 1,
'ZUNMBR',
'QLN', n, n, n, -1 ) )
380 wrkbl = max( wrkbl, 2*n+n*
381 $ ilaenv( 1,
'ZUNMBR',
'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,
'ZGEBRD',
' ', m, n,
393 maxwrk = max( maxwrk, 2*n+n*
394 $ ilaenv( 1,
'ZUNGBR',
'P', n, n, n, -1 ) )
395 maxwrk = max( maxwrk, 2*n+n*
396 $ ilaenv( 1,
'ZUNGBR',
'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,
'ZUNGBR',
'P', n, n, n, -1 ) )
402 maxwrk = max( maxwrk, 2*n+n*
403 $ ilaenv( 1,
'ZUNGBR',
'Q', m, n, n, -1 ) )
404 ELSE IF( wntqa )
THEN
405 maxwrk = max( maxwrk, 2*n+n*
406 $ ilaenv( 1,
'ZUNGBR',
'P', n, n, n, -1 ) )
407 maxwrk = max( maxwrk, 2*n+m*
408 $ ilaenv( 1,
'ZUNGBR',
'Q', m, m, n, -1 ) )
414 maxwrk = 2*n + ( m+n )*ilaenv( 1,
'ZGEBRD',
' ', m, n,
418 maxwrk = max( maxwrk, 2*n+n*
419 $ ilaenv( 1,
'ZUNMBR',
'PRC', n, n, n, -1 ) )
420 maxwrk = max( maxwrk, 2*n+n*
421 $ ilaenv( 1,
'ZUNMBR',
'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,
'ZUNMBR',
'PRC', n, n, n, -1 ) )
427 maxwrk = max( maxwrk, 2*n+n*
428 $ ilaenv( 1,
'ZUNMBR',
'QLN', m, n, n, -1 ) )
429 ELSE IF( wntqa )
THEN
430 maxwrk = max( maxwrk, 2*n+n*
431 $ ilaenv( 1,
'ZUNGBR',
'PRC', n, n, n, -1 ) )
432 maxwrk = max( maxwrk, 2*n+m*
433 $ ilaenv( 1,
'ZUNGBR',
'QLN', m, m, n, -1 ) )
445 IF( n.GE.mnthr1 )
THEN
450 maxwrk = m + m*ilaenv( 1,
'ZGELQF',
' ', m, n, -1,
452 maxwrk = max( maxwrk, 2*m+2*m*
453 $ ilaenv( 1,
'ZGEBRD',
' ', m, m, -1, -1 ) )
455 ELSE IF( wntqo )
THEN
459 wrkbl = m + m*ilaenv( 1,
'ZGELQF',
' ', m, n, -1, -1 )
460 wrkbl = max( wrkbl, m+m*ilaenv( 1,
'ZUNGLQ',
' ', m,
462 wrkbl = max( wrkbl, 2*m+2*m*
463 $ ilaenv( 1,
'ZGEBRD',
' ', m, m, -1, -1 ) )
464 wrkbl = max( wrkbl, 2*m+m*
465 $ ilaenv( 1,
'ZUNMBR',
'PRC', m, m, m, -1 ) )
466 wrkbl = max( wrkbl, 2*m+m*
467 $ ilaenv( 1,
'ZUNMBR',
'QLN', m, m, m, -1 ) )
468 maxwrk = m*n + m*m + wrkbl
470 ELSE IF( wntqs )
THEN
474 wrkbl = m + m*ilaenv( 1,
'ZGELQF',
' ', m, n, -1, -1 )
475 wrkbl = max( wrkbl, m+m*ilaenv( 1,
'ZUNGLQ',
' ', m,
477 wrkbl = max( wrkbl, 2*m+2*m*
478 $ ilaenv( 1,
'ZGEBRD',
' ', m, m, -1, -1 ) )
479 wrkbl = max( wrkbl, 2*m+m*
480 $ ilaenv( 1,
'ZUNMBR',
'PRC', m, m, m, -1 ) )
481 wrkbl = max( wrkbl, 2*m+m*
482 $ ilaenv( 1,
'ZUNMBR',
'QLN', m, m, m, -1 ) )
485 ELSE IF( wntqa )
THEN
489 wrkbl = m + m*ilaenv( 1,
'ZGELQF',
' ', m, n, -1, -1 )
490 wrkbl = max( wrkbl, m+n*ilaenv( 1,
'ZUNGLQ',
' ', n,
492 wrkbl = max( wrkbl, 2*m+2*m*
493 $ ilaenv( 1,
'ZGEBRD',
' ', m, m, -1, -1 ) )
494 wrkbl = max( wrkbl, 2*m+m*
495 $ ilaenv( 1,
'ZUNMBR',
'PRC', m, m, m, -1 ) )
496 wrkbl = max( wrkbl, 2*m+m*
497 $ ilaenv( 1,
'ZUNMBR',
'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,
'ZGEBRD',
' ', m, n,
509 maxwrk = max( maxwrk, 2*m+m*
510 $ ilaenv( 1,
'ZUNGBR',
'P', m, n, m, -1 ) )
511 maxwrk = max( maxwrk, 2*m+m*
512 $ ilaenv( 1,
'ZUNGBR',
'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,
'ZUNGBR',
'P', m, n, m, -1 ) )
518 maxwrk = max( maxwrk, 2*m+m*
519 $ ilaenv( 1,
'ZUNGBR',
'Q', m, m, n, -1 ) )
520 ELSE IF( wntqa )
THEN
521 maxwrk = max( maxwrk, 2*m+n*
522 $ ilaenv( 1,
'ZUNGBR',
'P', n, n, m, -1 ) )
523 maxwrk = max( maxwrk, 2*m+m*
524 $ ilaenv( 1,
'ZUNGBR',
'Q', m, m, n, -1 ) )
530 maxwrk = 2*m + ( m+n )*ilaenv( 1,
'ZGEBRD',
' ', m, n,
534 maxwrk = max( maxwrk, 2*m+m*
535 $ ilaenv( 1,
'ZUNMBR',
'PRC', m, n, m, -1 ) )
536 maxwrk = max( maxwrk, 2*m+m*
537 $ ilaenv( 1,
'ZUNMBR',
'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,
'ZUNGBR',
'PRC', m, n, m, -1 ) )
543 maxwrk = max( maxwrk, 2*m+m*
544 $ ilaenv( 1,
'ZUNGBR',
'QLN', m, m, n, -1 ) )
545 ELSE IF( wntqa )
THEN
546 maxwrk = max( maxwrk, 2*m+n*
547 $ ilaenv( 1,
'ZUNGBR',
'PRC', n, n, m, -1 ) )
548 maxwrk = max( maxwrk, 2*m+m*
549 $ ilaenv( 1,
'ZUNGBR',
'QLN', m, m, n, -1 ) )
553 maxwrk = max( maxwrk, minwrk )
557 IF( lwork.LT.minwrk .AND. lwork.NE.lquerv )
564 CALL xerbla(
'ZGESDD', -info )
567 IF( lwork.EQ.lquerv )
569 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
576 smlnum = sqrt( dlamch(
'S' ) ) / eps
577 bignum = one / smlnum
581 anrm = zlange(
'M', m, n, a, lda, dum )
583 IF( anrm.GT.zero .AND. anrm.LT.smlnum )
THEN
585 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
586 ELSE IF( anrm.GT.bignum )
THEN
588 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
597 IF( m.GE.mnthr1 )
THEN
611 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
612 $ lwork-nwork+1, ierr )
616 CALL zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
627 CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
628 $ work( itaup ), work( nwork ), lwork-nwork+1,
636 CALL dbdsdc(
'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 zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
667 $ lwork-nwork+1, ierr )
671 CALL zlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
672 CALL zlaset(
'L', n-1, n-1, czero, czero, work( ir+1 ),
679 CALL zungqr( m, n, n, a, lda, work( itau ),
680 $ work( nwork ), lwork-nwork+1, ierr )
690 CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
691 $ work( itauq ), work( itaup ), work( nwork ),
692 $ lwork-nwork+1, ierr )
703 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
704 $ n, rwork( irvt ), n, dum, idum,
705 $ rwork( nrwork ), iwork, info )
712 CALL zlacp2(
'F', n, n, rwork( iru ), n, work( iu ),
714 CALL zunmbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
715 $ work( itauq ), work( iu ), ldwrku,
716 $ work( nwork ), lwork-nwork+1, ierr )
723 CALL zlacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
724 CALL zunmbr(
'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 zgemm(
'N',
'N', chunk, n, n, cone, a( i, 1 ),
736 $ lda, work( iu ), ldwrku, czero,
737 $ work( ir ), ldwrkr )
738 CALL zlacpy(
'F', chunk, n, work( ir ), ldwrkr,
742 ELSE IF( wntqs )
THEN
760 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
761 $ lwork-nwork+1, ierr )
765 CALL zlacpy(
'U', n, n, a, lda, work( ir ), ldwrkr )
766 CALL zlaset(
'L', n-1, n-1, czero, czero, work( ir+1 ),
773 CALL zungqr( m, n, n, a, lda, work( itau ),
774 $ work( nwork ), lwork-nwork+1, ierr )
784 CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
785 $ work( itauq ), work( itaup ), work( nwork ),
786 $ lwork-nwork+1, ierr )
797 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
798 $ n, rwork( irvt ), n, dum, idum,
799 $ rwork( nrwork ), iwork, info )
806 CALL zlacp2(
'F', n, n, rwork( iru ), n, u, ldu )
807 CALL zunmbr(
'Q',
'L',
'N', n, n, n, work( ir ), ldwrkr,
808 $ work( itauq ), u, ldu, work( nwork ),
809 $ lwork-nwork+1, ierr )
816 CALL zlacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
817 CALL zunmbr(
'P',
'R',
'C', n, n, n, work( ir ), ldwrkr,
818 $ work( itaup ), vt, ldvt, work( nwork ),
819 $ lwork-nwork+1, ierr )
826 CALL zlacpy(
'F', n, n, u, ldu, work( ir ), ldwrkr )
827 CALL zgemm(
'N',
'N', m, n, n, cone, a, lda, work( ir ),
828 $ ldwrkr, czero, u, ldu )
830 ELSE IF( wntqa )
THEN
848 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
849 $ lwork-nwork+1, ierr )
850 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
856 CALL zungqr( m, m, n, u, ldu, work( itau ),
857 $ work( nwork ), lwork-nwork+1, ierr )
861 CALL zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
872 CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
873 $ work( itaup ), work( nwork ), lwork-nwork+1,
885 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
886 $ n, rwork( irvt ), n, dum, idum,
887 $ rwork( nrwork ), iwork, info )
894 CALL zlacp2(
'F', n, n, rwork( iru ), n, work( iu ),
896 CALL zunmbr(
'Q',
'L',
'N', n, n, n, a, lda,
897 $ work( itauq ), work( iu ), ldwrku,
898 $ work( nwork ), lwork-nwork+1, ierr )
905 CALL zlacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
906 CALL zunmbr(
'P',
'R',
'C', n, n, n, a, lda,
907 $ work( itaup ), vt, ldvt, work( nwork ),
908 $ lwork-nwork+1, ierr )
915 CALL zgemm(
'N',
'N', m, n, n, cone, u, ldu, work( iu ),
916 $ ldwrku, czero, a, lda )
920 CALL zlacpy(
'F', m, n, a, lda, u, ldu )
924 ELSE IF( m.GE.mnthr2 )
THEN
942 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
943 $ work( itaup ), work( nwork ), lwork-nwork+1,
951 CALL dbdsdc(
'U',
'N', n, s, rwork( ie ), dum, 1, dum, 1,
952 $ dum, idum, rwork( nrwork ), iwork, info )
953 ELSE IF( wntqo )
THEN
963 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
964 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
965 $ work( nwork ), lwork-nwork+1, ierr )
971 CALL zungbr(
'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 dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
994 $ n, rwork( irvt ), n, dum, idum,
995 $ rwork( nrwork ), iwork, info )
1002 CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt,
1003 $ work( iu ), ldwrku, rwork( nrwork ) )
1004 CALL zlacpy(
'F', n, n, work( iu ), ldwrku, vt, ldvt )
1012 DO 20 i = 1, m, ldwrku
1013 chunk = min( m-i+1, ldwrku )
1014 CALL zlacrm( chunk, n, a( i, 1 ), lda, rwork( iru ),
1015 $ n, work( iu ), ldwrku, rwork( nrwork ) )
1016 CALL zlacpy(
'F', chunk, n, work( iu ), ldwrku,
1020 ELSE IF( wntqs )
THEN
1026 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
1027 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1028 $ work( nwork ), lwork-nwork+1, ierr )
1034 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1035 CALL zungbr(
'Q', m, n, n, u, ldu, work( itauq ),
1036 $ work( nwork ), lwork-nwork+1, ierr )
1047 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1048 $ n, rwork( irvt ), n, dum, idum,
1049 $ rwork( nrwork ), iwork, info )
1056 CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1058 CALL zlacpy(
'F', n, n, a, lda, vt, ldvt )
1066 CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1068 CALL zlacpy(
'F', m, n, a, lda, u, ldu )
1075 CALL zlacpy(
'U', n, n, a, lda, vt, ldvt )
1076 CALL zungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1077 $ work( nwork ), lwork-nwork+1, ierr )
1083 CALL zlacpy(
'L', m, n, a, lda, u, ldu )
1084 CALL zungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1085 $ work( nwork ), lwork-nwork+1, ierr )
1096 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1097 $ n, rwork( irvt ), n, dum, idum,
1098 $ rwork( nrwork ), iwork, info )
1105 CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1107 CALL zlacpy(
'F', n, n, a, lda, vt, ldvt )
1115 CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1117 CALL zlacpy(
'F', m, n, a, lda, u, ldu )
1138 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1139 $ work( itaup ), work( nwork ), lwork-nwork+1,
1147 CALL dbdsdc(
'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 dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1174 $ n, rwork( irvt ), n, dum, idum,
1175 $ rwork( nrwork ), iwork, info )
1182 CALL zlacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1183 CALL zunmbr(
'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 zlaset(
'F', m, n, czero, czero, work( iu ),
1197 CALL zlacp2(
'F', n, n, rwork( iru ), n, work( iu ),
1199 CALL zunmbr(
'Q',
'L',
'N', m, n, n, a, lda,
1200 $ work( itauq ), work( iu ), ldwrku,
1201 $ work( nwork ), lwork-nwork+1, ierr )
1202 CALL zlacpy(
'F', m, n, work( iu ), ldwrku, a, lda )
1209 CALL zungbr(
'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 zlacrm( chunk, n, a( i, 1 ), lda,
1221 $ rwork( iru ), n, work( iu ), ldwrku,
1223 CALL zlacpy(
'F', chunk, n, work( iu ), ldwrku,
1228 ELSE IF( wntqs )
THEN
1239 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1240 $ n, rwork( irvt ), n, dum, idum,
1241 $ rwork( nrwork ), iwork, info )
1248 CALL zlaset(
'F', m, n, czero, czero, u, ldu )
1249 CALL zlacp2(
'F', n, n, rwork( iru ), n, u, ldu )
1250 CALL zunmbr(
'Q',
'L',
'N', m, n, n, a, lda,
1251 $ work( itauq ), u, ldu, work( nwork ),
1252 $ lwork-nwork+1, ierr )
1259 CALL zlacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1260 CALL zunmbr(
'P',
'R',
'C', n, n, n, a, lda,
1261 $ work( itaup ), vt, ldvt, work( nwork ),
1262 $ lwork-nwork+1, ierr )
1274 CALL dbdsdc(
'U',
'I', n, s, rwork( ie ), rwork( iru ),
1275 $ n, rwork( irvt ), n, dum, idum,
1276 $ rwork( nrwork ), iwork, info )
1280 CALL zlaset(
'F', m, m, czero, czero, u, ldu )
1282 CALL zlaset(
'F', m-n, m-n, czero, cone,
1283 $ u( n+1, n+1 ), ldu )
1291 CALL zlacp2(
'F', n, n, rwork( iru ), n, u, ldu )
1292 CALL zunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
1293 $ work( itauq ), u, ldu, work( nwork ),
1294 $ lwork-nwork+1, ierr )
1301 CALL zlacp2(
'F', n, n, rwork( irvt ), n, vt, ldvt )
1302 CALL zunmbr(
'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 zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1330 $ lwork-nwork+1, ierr )
1334 CALL zlaset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
1345 CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1346 $ work( itaup ), work( nwork ), lwork-nwork+1,
1354 CALL dbdsdc(
'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 zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1390 $ lwork-nwork+1, ierr )
1394 CALL zlacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1395 CALL zlaset(
'U', m-1, m-1, czero, czero,
1396 $ work( il+ldwrkl ), ldwrkl )
1402 CALL zunglq( m, n, m, a, lda, work( itau ),
1403 $ work( nwork ), lwork-nwork+1, ierr )
1413 CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1414 $ work( itauq ), work( itaup ), work( nwork ),
1415 $ lwork-nwork+1, ierr )
1426 CALL dbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1427 $ m, rwork( irvt ), m, dum, idum,
1428 $ rwork( nrwork ), iwork, info )
1435 CALL zlacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1436 CALL zunmbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1437 $ work( itauq ), u, ldu, work( nwork ),
1438 $ lwork-nwork+1, ierr )
1445 CALL zlacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1447 CALL zunmbr(
'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 zgemm(
'N',
'N', m, blk, m, cone, work( ivt ), m,
1459 $ a( 1, i ), lda, czero, work( il ),
1461 CALL zlacpy(
'F', m, blk, work( il ), ldwrkl,
1465 ELSE IF( wntqs )
THEN
1476 itau = il + ldwrkl*m
1483 CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1484 $ lwork-nwork+1, ierr )
1488 CALL zlacpy(
'L', m, m, a, lda, work( il ), ldwrkl )
1489 CALL zlaset(
'U', m-1, m-1, czero, czero,
1490 $ work( il+ldwrkl ), ldwrkl )
1496 CALL zunglq( m, n, m, a, lda, work( itau ),
1497 $ work( nwork ), lwork-nwork+1, ierr )
1507 CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1508 $ work( itauq ), work( itaup ), work( nwork ),
1509 $ lwork-nwork+1, ierr )
1520 CALL dbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1521 $ m, rwork( irvt ), m, dum, idum,
1522 $ rwork( nrwork ), iwork, info )
1529 CALL zlacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1530 CALL zunmbr(
'Q',
'L',
'N', m, m, m, work( il ), ldwrkl,
1531 $ work( itauq ), u, ldu, work( nwork ),
1532 $ lwork-nwork+1, ierr )
1539 CALL zlacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
1540 CALL zunmbr(
'P',
'R',
'C', m, m, m, work( il ), ldwrkl,
1541 $ work( itaup ), vt, ldvt, work( nwork ),
1542 $ lwork-nwork+1, ierr )
1549 CALL zlacpy(
'F', m, m, vt, ldvt, work( il ), ldwrkl )
1550 CALL zgemm(
'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 zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1572 $ lwork-nwork+1, ierr )
1573 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
1579 CALL zunglq( n, n, m, vt, ldvt, work( itau ),
1580 $ work( nwork ), lwork-nwork+1, ierr )
1584 CALL zlaset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
1595 CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1596 $ work( itaup ), work( nwork ), lwork-nwork+1,
1608 CALL dbdsdc(
'U',
'I', m, s, rwork( ie ), rwork( iru ),
1609 $ m, rwork( irvt ), m, dum, idum,
1610 $ rwork( nrwork ), iwork, info )
1617 CALL zlacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1618 CALL zunmbr(
'Q',
'L',
'N', m, m, m, a, lda,
1619 $ work( itauq ), u, ldu, work( nwork ),
1620 $ lwork-nwork+1, ierr )
1627 CALL zlacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1629 CALL zunmbr(
'P',
'R',
'C', m, m, m, a, lda,
1630 $ work( itaup ), work( ivt ), ldwkvt,
1631 $ work( nwork ), lwork-nwork+1, ierr )
1638 CALL zgemm(
'N',
'N', m, n, m, cone, work( ivt ), ldwkvt,
1639 $ vt, ldvt, czero, a, lda )
1643 CALL zlacpy(
'F', m, n, a, lda, vt, ldvt )
1647 ELSE IF( n.GE.mnthr2 )
THEN
1666 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1667 $ work( itaup ), work( nwork ), lwork-nwork+1,
1676 CALL dbdsdc(
'L',
'N', m, s, rwork( ie ), dum, 1, dum, 1,
1677 $ dum, idum, rwork( nrwork ), iwork, info )
1678 ELSE IF( wntqo )
THEN
1688 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
1689 CALL zungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1690 $ work( nwork ), lwork-nwork+1, ierr )
1696 CALL zungbr(
'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 dbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1721 $ m, rwork( irvt ), m, dum, idum,
1722 $ rwork( nrwork ), iwork, info )
1729 CALL zlacrm( m, m, u, ldu, rwork( iru ), m, work( ivt ),
1730 $ ldwkvt, rwork( nrwork ) )
1731 CALL zlacpy(
'F', m, m, work( ivt ), ldwkvt, u, ldu )
1739 DO 50 i = 1, n, chunk
1740 blk = min( n-i+1, chunk )
1741 CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ), lda,
1742 $ work( ivt ), ldwkvt, rwork( nrwork ) )
1743 CALL zlacpy(
'F', m, blk, work( ivt ), ldwkvt,
1746 ELSE IF( wntqs )
THEN
1752 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
1753 CALL zungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1754 $ work( nwork ), lwork-nwork+1, ierr )
1760 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
1761 CALL zungbr(
'P', m, n, m, vt, ldvt, work( itaup ),
1762 $ work( nwork ), lwork-nwork+1, ierr )
1773 CALL dbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1774 $ m, rwork( irvt ), m, dum, idum,
1775 $ rwork( nrwork ), iwork, info )
1782 CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1784 CALL zlacpy(
'F', m, m, a, lda, u, ldu )
1792 CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1794 CALL zlacpy(
'F', m, n, a, lda, vt, ldvt )
1801 CALL zlacpy(
'L', m, m, a, lda, u, ldu )
1802 CALL zungbr(
'Q', m, m, n, u, ldu, work( itauq ),
1803 $ work( nwork ), lwork-nwork+1, ierr )
1809 CALL zlacpy(
'U', m, n, a, lda, vt, ldvt )
1810 CALL zungbr(
'P', n, n, m, vt, ldvt, work( itaup ),
1811 $ work( nwork ), lwork-nwork+1, ierr )
1822 CALL dbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1823 $ m, rwork( irvt ), m, dum, idum,
1824 $ rwork( nrwork ), iwork, info )
1831 CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1833 CALL zlacpy(
'F', m, m, a, lda, u, ldu )
1840 CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1842 CALL zlacpy(
'F', m, n, a, lda, vt, ldvt )
1863 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1864 $ work( itaup ), work( nwork ), lwork-nwork+1,
1872 CALL dbdsdc(
'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 zlaset(
'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 dbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1902 $ m, rwork( irvt ), m, dum, idum,
1903 $ rwork( nrwork ), iwork, info )
1910 CALL zlacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1911 CALL zunmbr(
'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 zlacp2(
'F', m, m, rwork( irvt ), m, work( ivt ),
1925 CALL zunmbr(
'P',
'R',
'C', m, n, m, a, lda,
1926 $ work( itaup ), work( ivt ), ldwkvt,
1927 $ work( nwork ), lwork-nwork+1, ierr )
1928 CALL zlacpy(
'F', m, n, work( ivt ), ldwkvt, a, lda )
1935 CALL zungbr(
'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 zlarcm( m, blk, rwork( irvt ), m, a( 1, i ),
1947 $ lda, work( ivt ), ldwkvt,
1949 CALL zlacpy(
'F', m, blk, work( ivt ), ldwkvt,
1953 ELSE IF( wntqs )
THEN
1964 CALL dbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
1965 $ m, rwork( irvt ), m, dum, idum,
1966 $ rwork( nrwork ), iwork, info )
1973 CALL zlacp2(
'F', m, m, rwork( iru ), m, u, ldu )
1974 CALL zunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
1975 $ work( itauq ), u, ldu, work( nwork ),
1976 $ lwork-nwork+1, ierr )
1983 CALL zlaset(
'F', m, n, czero, czero, vt, ldvt )
1984 CALL zlacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
1985 CALL zunmbr(
'P',
'R',
'C', m, n, m, a, lda,
1986 $ work( itaup ), vt, ldvt, work( nwork ),
1987 $ lwork-nwork+1, ierr )
2000 CALL dbdsdc(
'L',
'I', m, s, rwork( ie ), rwork( iru ),
2001 $ m, rwork( irvt ), m, dum, idum,
2002 $ rwork( nrwork ), iwork, info )
2009 CALL zlacp2(
'F', m, m, rwork( iru ), m, u, ldu )
2010 CALL zunmbr(
'Q',
'L',
'N', m, m, n, a, lda,
2011 $ work( itauq ), u, ldu, work( nwork ),
2012 $ lwork-nwork+1, ierr )
2016 CALL zlaset(
'F', n, n, czero, cone, vt, ldvt )
2023 CALL zlacp2(
'F', m, m, rwork( irvt ), m, vt, ldvt )
2024 CALL zunmbr(
'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 dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
2039 IF( info.NE.0 .AND. anrm.GT.bignum )
2040 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
2041 $ rwork( ie ), minmn, ierr )
2042 IF( anrm.LT.smlnum )
2043 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
2045 IF( info.NE.0 .AND. anrm.LT.smlnum )
2046 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
2047 $ rwork( ie ), minmn, ierr )
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlacp2(UPLO, M, N, A, LDA, B, LDB)
ZLACP2 copies all or part of a real two-dimensional array to a complex array.
subroutine zgesdd(JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
ZGESDD
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLACRM multiplies a complex matrix by a square real matrix.
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
subroutine zlarcm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLARCM copies all or part of a real two-dimensional array to a complex array.
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD