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