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