270
271
272
273
274
275
276 CHARACTER JOBU, JOBVT, RANGE
277 INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
278 DOUBLE PRECISION VL, VU
279
280
281 INTEGER IWORK( * )
282 DOUBLE PRECISION S( * ), RWORK( * )
283 COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
284 $ WORK( * )
285
286
287
288
289
290 COMPLEX*16 CZERO, CONE
291 parameter( czero = ( 0.0d0, 0.0d0 ),
292 $ cone = ( 1.0d0, 0.0d0 ) )
293 DOUBLE PRECISION ZERO, ONE
294 parameter( zero = 0.0d0, one = 1.0d0 )
295
296
297 CHARACTER JOBZ, RNGTGK
298 LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
299 INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
300 $ ITAU, ITAUP, ITAUQ, ITEMP, ITEMPR, ITGKZ,
301 $ IUTGK, J, K, MAXWRK, MINMN, MINWRK, MNTHR
302 DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
303
304
305 DOUBLE PRECISION DUM( 1 )
306
307
310
311
312 LOGICAL LSAME
313 INTEGER ILAENV
314 DOUBLE PRECISION DLAMCH, ZLANGE
316
317
318 INTRINSIC max, min, sqrt
319
320
321
322
323
324 ns = 0
325 info = 0
327 lquery = ( lwork.EQ.-1 )
328 minmn = min( m, n )
329
330 wantu =
lsame( jobu,
'V' )
331 wantvt =
lsame( jobvt,
'V' )
332 IF( wantu .OR. wantvt ) THEN
333 jobz = 'V'
334 ELSE
335 jobz = 'N'
336 END IF
337 alls =
lsame( range,
'A' )
338 vals =
lsame( range,
'V' )
339 inds =
lsame( range,
'I' )
340
341 info = 0
342 IF( .NOT.
lsame( jobu,
'V' ) .AND.
343 $ .NOT.
lsame( jobu,
'N' ) )
THEN
344 info = -1
345 ELSE IF( .NOT.
lsame( jobvt,
'V' ) .AND.
346 $ .NOT.
lsame( jobvt,
'N' ) )
THEN
347 info = -2
348 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
349 info = -3
350 ELSE IF( m.LT.0 ) THEN
351 info = -4
352 ELSE IF( n.LT.0 ) THEN
353 info = -5
354 ELSE IF( m.GT.lda ) THEN
355 info = -7
356 ELSE IF( minmn.GT.0 ) THEN
357 IF( vals ) THEN
358 IF( vl.LT.zero ) THEN
359 info = -8
360 ELSE IF( vu.LE.vl ) THEN
361 info = -9
362 END IF
363 ELSE IF( inds ) THEN
364 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
365 info = -10
366 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
367 info = -11
368 END IF
369 END IF
370 IF( info.EQ.0 ) THEN
371 IF( wantu .AND. ldu.LT.m ) THEN
372 info = -15
373 ELSE IF( wantvt ) THEN
374 IF( inds ) THEN
375 IF( ldvt.LT.iu-il+1 ) THEN
376 info = -17
377 END IF
378 ELSE IF( ldvt.LT.minmn ) THEN
379 info = -17
380 END IF
381 END IF
382 END IF
383 END IF
384
385
386
387
388
389
390
391
392 IF( info.EQ.0 ) THEN
393 minwrk = 1
394 maxwrk = 1
395 IF( minmn.GT.0 ) THEN
396 IF( m.GE.n ) THEN
397 mnthr =
ilaenv( 6,
'ZGESVD', jobu // jobvt, m, n, 0, 0 )
398 IF( m.GE.mnthr ) THEN
399
400
401
402 minwrk = n*(n+5)
403 maxwrk = n + n*
ilaenv(1,
'ZGEQRF',
' ',m,n,-1,-1)
404 maxwrk = max(maxwrk,
405 $ n*n+2*n+2*n*
ilaenv(1,
'ZGEBRD',
' ',n,n,-1,-1))
406 IF (wantu .OR. wantvt) THEN
407 maxwrk = max(maxwrk,
408 $ n*n+2*n+n*
ilaenv(1,
'ZUNMQR',
'LN',n,n,n,-1))
409 END IF
410 ELSE
411
412
413
414 minwrk = 3*n + m
415 maxwrk = 2*n + (m+n)*
ilaenv(1,
'ZGEBRD',
' ',m,n,-1,-1)
416 IF (wantu .OR. wantvt) THEN
417 maxwrk = max(maxwrk,
418 $ 2*n+n*
ilaenv(1,
'ZUNMQR',
'LN',n,n,n,-1))
419 END IF
420 END IF
421 ELSE
422 mnthr =
ilaenv( 6,
'ZGESVD', jobu // jobvt, m, n, 0, 0 )
423 IF( n.GE.mnthr ) THEN
424
425
426
427 minwrk = m*(m+5)
428 maxwrk = m + m*
ilaenv(1,
'ZGELQF',
' ',m,n,-1,-1)
429 maxwrk = max(maxwrk,
430 $ m*m+2*m+2*m*
ilaenv(1,
'ZGEBRD',
' ',m,m,-1,-1))
431 IF (wantu .OR. wantvt) THEN
432 maxwrk = max(maxwrk,
433 $ m*m+2*m+m*
ilaenv(1,
'ZUNMQR',
'LN',m,m,m,-1))
434 END IF
435 ELSE
436
437
438
439
440 minwrk = 3*m + n
441 maxwrk = 2*m + (m+n)*
ilaenv(1,
'ZGEBRD',
' ',m,n,-1,-1)
442 IF (wantu .OR. wantvt) THEN
443 maxwrk = max(maxwrk,
444 $ 2*m+m*
ilaenv(1,
'ZUNMQR',
'LN',m,m,m,-1))
445 END IF
446 END IF
447 END IF
448 END IF
449 maxwrk = max( maxwrk, minwrk )
450 work( 1 ) = dcmplx( dble( maxwrk ), zero )
451
452 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
453 info = -19
454 END IF
455 END IF
456
457 IF( info.NE.0 ) THEN
458 CALL xerbla(
'ZGESVDX', -info )
459 RETURN
460 ELSE IF( lquery ) THEN
461 RETURN
462 END IF
463
464
465
466 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
467 RETURN
468 END IF
469
470
471
472 IF( alls ) THEN
473 rngtgk = 'I'
474 iltgk = 1
475 iutgk = min( m, n )
476 ELSE IF( inds ) THEN
477 rngtgk = 'I'
478 iltgk = il
479 iutgk = iu
480 ELSE
481 rngtgk = 'V'
482 iltgk = 0
483 iutgk = 0
484 END IF
485
486
487
489 smlnum = sqrt(
dlamch(
'S' ) ) / eps
490 bignum = one / smlnum
491
492
493
494 anrm =
zlange(
'M', m, n, a, lda, dum )
495 iscl = 0
496 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
497 iscl = 1
498 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
499 ELSE IF( anrm.GT.bignum ) THEN
500 iscl = 1
501 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
502 END IF
503
504 IF( m.GE.n ) THEN
505
506
507
508
509
510 IF( m.GE.mnthr ) THEN
511
512
513
514
515
516
517
518
519
520 itau = 1
521 itemp = itau + n
522 CALL zgeqrf( m, n, a, lda, work( itau ), work( itemp ),
523 $ lwork-itemp+1, info )
524
525
526
527
528 iqrf = itemp
529 itauq = itemp + n*n
530 itaup = itauq + n
531 itemp = itaup + n
532 id = 1
533 ie = id + n
534 itgkz = ie + n
535 CALL zlacpy(
'U', n, n, a, lda, work( iqrf ), n )
536 CALL zlaset(
'L', n-1, n-1, czero, czero,
537 $ work( iqrf+1 ), n )
538 CALL zgebrd( n, n, work( iqrf ), n, rwork( id ),
539 $ rwork( ie ), work( itauq ), work( itaup ),
540 $ work( itemp ), lwork-itemp+1, info )
541 itempr = itgkz + n*(n*2+1)
542
543
544
545
546 CALL dbdsvdx(
'U', jobz, rngtgk, n, rwork( id ),
547 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
548 $ rwork( itgkz ), n*2, rwork( itempr ),
549 $ iwork, info)
550
551
552
553 IF( wantu ) THEN
554 k = itgkz
555 DO i = 1, ns
556 DO j = 1, n
557 u( j, i ) = dcmplx( rwork( k ), zero )
558 k = k + 1
559 END DO
560 k = k + n
561 END DO
562 CALL zlaset(
'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
563
564
565
566
567 CALL zunmbr(
'Q',
'L',
'N', n, ns, n, work( iqrf ), n,
568 $ work( itauq ), u, ldu, work( itemp ),
569 $ lwork-itemp+1, info )
570
571
572
573
574 CALL zunmqr(
'L',
'N', m, ns, n, a, lda,
575 $ work( itau ), u, ldu, work( itemp ),
576 $ lwork-itemp+1, info )
577 END IF
578
579
580
581 IF( wantvt) THEN
582 k = itgkz + n
583 DO i = 1, ns
584 DO j = 1, n
585 vt( i, j ) = dcmplx( rwork( k ), zero )
586 k = k + 1
587 END DO
588 k = k + n
589 END DO
590
591
592
593
594 CALL zunmbr(
'P',
'R',
'C', ns, n, n, work( iqrf ), n,
595 $ work( itaup ), vt, ldvt, work( itemp ),
596 $ lwork-itemp+1, info )
597 END IF
598 ELSE
599
600
601
602
603
604
605
606
607
608 itauq = 1
609 itaup = itauq + n
610 itemp = itaup + n
611 id = 1
612 ie = id + n
613 itgkz = ie + n
614 CALL zgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
615 $ work( itauq ), work( itaup ), work( itemp ),
616 $ lwork-itemp+1, info )
617 itempr = itgkz + n*(n*2+1)
618
619
620
621
622 CALL dbdsvdx(
'U', jobz, rngtgk, n, rwork( id ),
623 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
624 $ rwork( itgkz ), n*2, rwork( itempr ),
625 $ iwork, info)
626
627
628
629 IF( wantu ) THEN
630 k = itgkz
631 DO i = 1, ns
632 DO j = 1, n
633 u( j, i ) = dcmplx( rwork( k ), zero )
634 k = k + 1
635 END DO
636 k = k + n
637 END DO
638 CALL zlaset(
'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
639
640
641
642
643 CALL zunmbr(
'Q',
'L',
'N', m, ns, n, a, lda,
644 $ work( itauq ), u, ldu, work( itemp ),
645 $ lwork-itemp+1, ierr )
646 END IF
647
648
649
650 IF( wantvt) THEN
651 k = itgkz + n
652 DO i = 1, ns
653 DO j = 1, n
654 vt( i, j ) = dcmplx( rwork( k ), zero )
655 k = k + 1
656 END DO
657 k = k + n
658 END DO
659
660
661
662
663 CALL zunmbr(
'P',
'R',
'C', ns, n, n, a, lda,
664 $ work( itaup ), vt, ldvt, work( itemp ),
665 $ lwork-itemp+1, ierr )
666 END IF
667 END IF
668 ELSE
669
670
671
672
673 IF( n.GE.mnthr ) THEN
674
675
676
677
678
679
680
681
682
683 itau = 1
684 itemp = itau + m
685 CALL zgelqf( m, n, a, lda, work( itau ), work( itemp ),
686 $ lwork-itemp+1, info )
687
688
689
690
691 ilqf = itemp
692 itauq = ilqf + m*m
693 itaup = itauq + m
694 itemp = itaup + m
695 id = 1
696 ie = id + m
697 itgkz = ie + m
698 CALL zlacpy(
'L', m, m, a, lda, work( ilqf ), m )
699 CALL zlaset(
'U', m-1, m-1, czero, czero,
700 $ work( ilqf+m ), m )
701 CALL zgebrd( m, m, work( ilqf ), m, rwork( id ),
702 $ rwork( ie ), work( itauq ), work( itaup ),
703 $ work( itemp ), lwork-itemp+1, info )
704 itempr = itgkz + m*(m*2+1)
705
706
707
708
709 CALL dbdsvdx(
'U', jobz, rngtgk, m, rwork( id ),
710 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
711 $ rwork( itgkz ), m*2, rwork( itempr ),
712 $ iwork, info)
713
714
715
716 IF( wantu ) THEN
717 k = itgkz
718 DO i = 1, ns
719 DO j = 1, m
720 u( j, i ) = dcmplx( rwork( k ), zero )
721 k = k + 1
722 END DO
723 k = k + m
724 END DO
725
726
727
728
729 CALL zunmbr(
'Q',
'L',
'N', m, ns, m, work( ilqf ), m,
730 $ work( itauq ), u, ldu, work( itemp ),
731 $ lwork-itemp+1, info )
732 END IF
733
734
735
736 IF( wantvt) THEN
737 k = itgkz + m
738 DO i = 1, ns
739 DO j = 1, m
740 vt( i, j ) = dcmplx( rwork( k ), zero )
741 k = k + 1
742 END DO
743 k = k + m
744 END DO
745 CALL zlaset(
'A', ns, n-m, czero, czero,
746 $ vt( 1,m+1 ), ldvt )
747
748
749
750
751 CALL zunmbr(
'P',
'R',
'C', ns, m, m, work( ilqf ), m,
752 $ work( itaup ), vt, ldvt, work( itemp ),
753 $ lwork-itemp+1, info )
754
755
756
757
758 CALL zunmlq(
'R',
'N', ns, n, m, a, lda,
759 $ work( itau ), vt, ldvt, work( itemp ),
760 $ lwork-itemp+1, info )
761 END IF
762 ELSE
763
764
765
766
767
768
769
770
771
772 itauq = 1
773 itaup = itauq + m
774 itemp = itaup + m
775 id = 1
776 ie = id + m
777 itgkz = ie + m
778 CALL zgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
779 $ work( itauq ), work( itaup ), work( itemp ),
780 $ lwork-itemp+1, info )
781 itempr = itgkz + m*(m*2+1)
782
783
784
785
786 CALL dbdsvdx(
'L', jobz, rngtgk, m, rwork( id ),
787 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
788 $ rwork( itgkz ), m*2, rwork( itempr ),
789 $ iwork, info)
790
791
792
793 IF( wantu ) THEN
794 k = itgkz
795 DO i = 1, ns
796 DO j = 1, m
797 u( j, i ) = dcmplx( rwork( k ), zero )
798 k = k + 1
799 END DO
800 k = k + m
801 END DO
802
803
804
805
806 CALL zunmbr(
'Q',
'L',
'N', m, ns, n, a, lda,
807 $ work( itauq ), u, ldu, work( itemp ),
808 $ lwork-itemp+1, info )
809 END IF
810
811
812
813 IF( wantvt) THEN
814 k = itgkz + m
815 DO i = 1, ns
816 DO j = 1, m
817 vt( i, j ) = dcmplx( rwork( k ), zero )
818 k = k + 1
819 END DO
820 k = k + m
821 END DO
822 CALL zlaset(
'A', ns, n-m, czero, czero,
823 $ vt( 1,m+1 ), ldvt )
824
825
826
827
828 CALL zunmbr(
'P',
'R',
'C', ns, n, m, a, lda,
829 $ work( itaup ), vt, ldvt, work( itemp ),
830 $ lwork-itemp+1, info )
831 END IF
832 END IF
833 END IF
834
835
836
837 IF( iscl.EQ.1 ) THEN
838 IF( anrm.GT.bignum )
839 $
CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1,
840 $ s, minmn, info )
841 IF( anrm.LT.smlnum )
842 $
CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1,
843 $ s, minmn, info )
844 END IF
845
846
847
848 work( 1 ) = dcmplx( dble( maxwrk ), zero )
849
850 RETURN
851
852
853
subroutine xerbla(srname, info)
subroutine dbdsvdx(uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
DBDSVDX
subroutine zgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
ZGEBRD
subroutine zgelqf(m, n, a, lda, tau, work, lwork, info)
ZGELQF
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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 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 zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME
subroutine zunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMBR
subroutine zunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMLQ
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR