241
242
243
244
245
246
247 CHARACTER UPLO
248 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
249
250
251 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
252 $ VT( LDVT, * ), WORK( * )
253
254
255
256
257
258 DOUBLE PRECISION ZERO
259 parameter( zero = 0.0d0 )
260 DOUBLE PRECISION ONE
261 parameter( one = 1.0d0 )
262 DOUBLE PRECISION NEGONE
263 parameter( negone = -1.0d0 )
264 DOUBLE PRECISION HNDRTH
265 parameter( hndrth = 0.01d0 )
266 DOUBLE PRECISION TEN
267 parameter( ten = 10.0d0 )
268 DOUBLE PRECISION HNDRD
269 parameter( hndrd = 100.0d0 )
270 DOUBLE PRECISION MEIGTH
271 parameter( meigth = -0.125d0 )
272 INTEGER MAXITR
273 parameter( maxitr = 6 )
274
275
276 LOGICAL LOWER, ROTATE
277 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
278 $ MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
279 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
280 $ OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
281 $ SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
282 $ SN, THRESH, TOL, TOLMUL, UNFL
283
284
285 LOGICAL LSAME
286 DOUBLE PRECISION DLAMCH
288
289
292
293
294 INTRINSIC abs, dble, max, min, sign, sqrt
295
296
297
298
299
300 info = 0
301 lower =
lsame( uplo,
'L' )
302 IF( .NOT.
lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
303 info = -1
304 ELSE IF( n.LT.0 ) THEN
305 info = -2
306 ELSE IF( ncvt.LT.0 ) THEN
307 info = -3
308 ELSE IF( nru.LT.0 ) THEN
309 info = -4
310 ELSE IF( ncc.LT.0 ) THEN
311 info = -5
312 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
313 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
314 info = -9
315 ELSE IF( ldu.LT.max( 1, nru ) ) THEN
316 info = -11
317 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
318 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
319 info = -13
320 END IF
321 IF( info.NE.0 ) THEN
322 CALL xerbla(
'DBDSQR', -info )
323 RETURN
324 END IF
325 IF( n.EQ.0 )
326 $ RETURN
327 IF( n.EQ.1 )
328 $ GO TO 160
329
330
331
332 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
333
334
335
336 IF( .NOT.rotate ) THEN
337 CALL dlasq1( n, d, e, work, info )
338
339
340
341 IF( info .NE. 2 ) RETURN
342 info = 0
343 END IF
344
345 nm1 = n - 1
346 nm12 = nm1 + nm1
347 nm13 = nm12 + nm1
348 idir = 0
349
350
351
353 unfl =
dlamch(
'Safe minimum' )
354
355
356
357
358 IF( lower ) THEN
359 DO 10 i = 1, n - 1
360 CALL dlartg( d( i ), e( i ), cs, sn, r )
361 d( i ) = r
362 e( i ) = sn*d( i+1 )
363 d( i+1 ) = cs*d( i+1 )
364 work( i ) = cs
365 work( nm1+i ) = sn
366 10 CONTINUE
367
368
369
370 IF( nru.GT.0 )
371 $
CALL dlasr(
'R',
'V',
'F', nru, n, work( 1 ), work( n ), u,
372 $ ldu )
373 IF( ncc.GT.0 )
374 $
CALL dlasr(
'L',
'V',
'F', n, ncc, work( 1 ), work( n ), c,
375 $ ldc )
376 END IF
377
378
379
380
381
382 tolmul = max( ten, min( hndrd, eps**meigth ) )
383 tol = tolmul*eps
384
385
386
387 smax = zero
388 DO 20 i = 1, n
389 smax = max( smax, abs( d( i ) ) )
390 20 CONTINUE
391 DO 30 i = 1, n - 1
392 smax = max( smax, abs( e( i ) ) )
393 30 CONTINUE
394 sminl = zero
395 IF( tol.GE.zero ) THEN
396
397
398
399 sminoa = abs( d( 1 ) )
400 IF( sminoa.EQ.zero )
401 $ GO TO 50
402 mu = sminoa
403 DO 40 i = 2, n
404 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
405 sminoa = min( sminoa, mu )
406 IF( sminoa.EQ.zero )
407 $ GO TO 50
408 40 CONTINUE
409 50 CONTINUE
410 sminoa = sminoa / sqrt( dble( n ) )
411 thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
412 ELSE
413
414
415
416 thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
417 END IF
418
419
420
421
422
423 maxitdivn = maxitr*n
424 iterdivn = 0
425 iter = -1
426 oldll = -1
427 oldm = -1
428
429
430
431 m = n
432
433
434
435 60 CONTINUE
436
437
438
439 IF( m.LE.1 )
440 $ GO TO 160
441
442 IF( iter.GE.n ) THEN
443 iter = iter - n
444 iterdivn = iterdivn + 1
445 IF( iterdivn.GE.maxitdivn )
446 $ GO TO 200
447 END IF
448
449
450
451 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
452 $ d( m ) = zero
453 smax = abs( d( m ) )
454 smin = smax
455 DO 70 lll = 1, m - 1
456 ll = m - lll
457 abss = abs( d( ll ) )
458 abse = abs( e( ll ) )
459 IF( tol.LT.zero .AND. abss.LE.thresh )
460 $ d( ll ) = zero
461 IF( abse.LE.thresh )
462 $ GO TO 80
463 smin = min( smin, abss )
464 smax = max( smax, abss, abse )
465 70 CONTINUE
466 ll = 0
467 GO TO 90
468 80 CONTINUE
469 e( ll ) = zero
470
471
472
473 IF( ll.EQ.m-1 ) THEN
474
475
476
477 m = m - 1
478 GO TO 60
479 END IF
480 90 CONTINUE
481 ll = ll + 1
482
483
484
485 IF( ll.EQ.m-1 ) THEN
486
487
488
489 CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
490 $ cosr, sinl, cosl )
491 d( m-1 ) = sigmx
492 e( m-1 ) = zero
493 d( m ) = sigmn
494
495
496
497 IF( ncvt.GT.0 )
498 $
CALL drot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
499 $ sinr )
500 IF( nru.GT.0 )
501 $
CALL drot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
502 IF( ncc.GT.0 )
503 $
CALL drot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
504 $ sinl )
505 m = m - 2
506 GO TO 60
507 END IF
508
509
510
511
512 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
513 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
514
515
516
517 idir = 1
518 ELSE
519
520
521
522 idir = 2
523 END IF
524 END IF
525
526
527
528 IF( idir.EQ.1 ) THEN
529
530
531
532
533 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
534 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
535 e( m-1 ) = zero
536 GO TO 60
537 END IF
538
539 IF( tol.GE.zero ) THEN
540
541
542
543
544 mu = abs( d( ll ) )
545 sminl = mu
546 DO 100 lll = ll, m - 1
547 IF( abs( e( lll ) ).LE.tol*mu ) THEN
548 e( lll ) = zero
549 GO TO 60
550 END IF
551 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
552 sminl = min( sminl, mu )
553 100 CONTINUE
554 END IF
555
556 ELSE
557
558
559
560
561 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
562 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
563 e( ll ) = zero
564 GO TO 60
565 END IF
566
567 IF( tol.GE.zero ) THEN
568
569
570
571
572 mu = abs( d( m ) )
573 sminl = mu
574 DO 110 lll = m - 1, ll, -1
575 IF( abs( e( lll ) ).LE.tol*mu ) THEN
576 e( lll ) = zero
577 GO TO 60
578 END IF
579 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
580 sminl = min( sminl, mu )
581 110 CONTINUE
582 END IF
583 END IF
584 oldll = ll
585 oldm = m
586
587
588
589
590 IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
591 $ max( eps, hndrth*tol ) ) THEN
592
593
594
595 shift = zero
596 ELSE
597
598
599
600 IF( idir.EQ.1 ) THEN
601 sll = abs( d( ll ) )
602 CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
603 ELSE
604 sll = abs( d( m ) )
605 CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
606 END IF
607
608
609
610 IF( sll.GT.zero ) THEN
611 IF( ( shift / sll )**2.LT.eps )
612 $ shift = zero
613 END IF
614 END IF
615
616
617
618 iter = iter + m - ll
619
620
621
622 IF( shift.EQ.zero ) THEN
623 IF( idir.EQ.1 ) THEN
624
625
626
627
628 cs = one
629 oldcs = one
630 DO 120 i = ll, m - 1
631 CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
632 IF( i.GT.ll )
633 $ e( i-1 ) = oldsn*r
634 CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
635 work( i-ll+1 ) = cs
636 work( i-ll+1+nm1 ) = sn
637 work( i-ll+1+nm12 ) = oldcs
638 work( i-ll+1+nm13 ) = oldsn
639 120 CONTINUE
640 h = d( m )*cs
641 d( m ) = h*oldcs
642 e( m-1 ) = h*oldsn
643
644
645
646 IF( ncvt.GT.0 )
647 $
CALL dlasr(
'L',
'V',
'F', m-ll+1, ncvt, work( 1 ),
648 $ work( n ), vt( ll, 1 ), ldvt )
649 IF( nru.GT.0 )
650 $
CALL dlasr(
'R',
'V',
'F', nru, m-ll+1, work( nm12+1 ),
651 $ work( nm13+1 ), u( 1, ll ), ldu )
652 IF( ncc.GT.0 )
653 $
CALL dlasr(
'L',
'V',
'F', m-ll+1, ncc, work( nm12+1 ),
654 $ work( nm13+1 ), c( ll, 1 ), ldc )
655
656
657
658 IF( abs( e( m-1 ) ).LE.thresh )
659 $ e( m-1 ) = zero
660
661 ELSE
662
663
664
665
666 cs = one
667 oldcs = one
668 DO 130 i = m, ll + 1, -1
669 CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
670 IF( i.LT.m )
671 $ e( i ) = oldsn*r
672 CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
673 work( i-ll ) = cs
674 work( i-ll+nm1 ) = -sn
675 work( i-ll+nm12 ) = oldcs
676 work( i-ll+nm13 ) = -oldsn
677 130 CONTINUE
678 h = d( ll )*cs
679 d( ll ) = h*oldcs
680 e( ll ) = h*oldsn
681
682
683
684 IF( ncvt.GT.0 )
685 $
CALL dlasr(
'L',
'V',
'B', m-ll+1, ncvt, work( nm12+1 ),
686 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
687 IF( nru.GT.0 )
688 $
CALL dlasr(
'R',
'V',
'B', nru, m-ll+1, work( 1 ),
689 $ work( n ), u( 1, ll ), ldu )
690 IF( ncc.GT.0 )
691 $
CALL dlasr(
'L',
'V',
'B', m-ll+1, ncc, work( 1 ),
692 $ work( n ), c( ll, 1 ), ldc )
693
694
695
696 IF( abs( e( ll ) ).LE.thresh )
697 $ e( ll ) = zero
698 END IF
699 ELSE
700
701
702
703 IF( idir.EQ.1 ) THEN
704
705
706
707
708 f = ( abs( d( ll ) )-shift )*
709 $ ( sign( one, d( ll ) )+shift / d( ll ) )
710 g = e( ll )
711 DO 140 i = ll, m - 1
712 CALL dlartg( f, g, cosr, sinr, r )
713 IF( i.GT.ll )
714 $ e( i-1 ) = r
715 f = cosr*d( i ) + sinr*e( i )
716 e( i ) = cosr*e( i ) - sinr*d( i )
717 g = sinr*d( i+1 )
718 d( i+1 ) = cosr*d( i+1 )
719 CALL dlartg( f, g, cosl, sinl, r )
720 d( i ) = r
721 f = cosl*e( i ) + sinl*d( i+1 )
722 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
723 IF( i.LT.m-1 ) THEN
724 g = sinl*e( i+1 )
725 e( i+1 ) = cosl*e( i+1 )
726 END IF
727 work( i-ll+1 ) = cosr
728 work( i-ll+1+nm1 ) = sinr
729 work( i-ll+1+nm12 ) = cosl
730 work( i-ll+1+nm13 ) = sinl
731 140 CONTINUE
732 e( m-1 ) = f
733
734
735
736 IF( ncvt.GT.0 )
737 $
CALL dlasr(
'L',
'V',
'F', m-ll+1, ncvt, work( 1 ),
738 $ work( n ), vt( ll, 1 ), ldvt )
739 IF( nru.GT.0 )
740 $
CALL dlasr(
'R',
'V',
'F', nru, m-ll+1, work( nm12+1 ),
741 $ work( nm13+1 ), u( 1, ll ), ldu )
742 IF( ncc.GT.0 )
743 $
CALL dlasr(
'L',
'V',
'F', m-ll+1, ncc, work( nm12+1 ),
744 $ work( nm13+1 ), c( ll, 1 ), ldc )
745
746
747
748 IF( abs( e( m-1 ) ).LE.thresh )
749 $ e( m-1 ) = zero
750
751 ELSE
752
753
754
755
756 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
757 $ d( m ) )
758 g = e( m-1 )
759 DO 150 i = m, ll + 1, -1
760 CALL dlartg( f, g, cosr, sinr, r )
761 IF( i.LT.m )
762 $ e( i ) = r
763 f = cosr*d( i ) + sinr*e( i-1 )
764 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
765 g = sinr*d( i-1 )
766 d( i-1 ) = cosr*d( i-1 )
767 CALL dlartg( f, g, cosl, sinl, r )
768 d( i ) = r
769 f = cosl*e( i-1 ) + sinl*d( i-1 )
770 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
771 IF( i.GT.ll+1 ) THEN
772 g = sinl*e( i-2 )
773 e( i-2 ) = cosl*e( i-2 )
774 END IF
775 work( i-ll ) = cosr
776 work( i-ll+nm1 ) = -sinr
777 work( i-ll+nm12 ) = cosl
778 work( i-ll+nm13 ) = -sinl
779 150 CONTINUE
780 e( ll ) = f
781
782
783
784 IF( abs( e( ll ) ).LE.thresh )
785 $ e( ll ) = zero
786
787
788
789 IF( ncvt.GT.0 )
790 $
CALL dlasr(
'L',
'V',
'B', m-ll+1, ncvt, work( nm12+1 ),
791 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
792 IF( nru.GT.0 )
793 $
CALL dlasr(
'R',
'V',
'B', nru, m-ll+1, work( 1 ),
794 $ work( n ), u( 1, ll ), ldu )
795 IF( ncc.GT.0 )
796 $
CALL dlasr(
'L',
'V',
'B', m-ll+1, ncc, work( 1 ),
797 $ work( n ), c( ll, 1 ), ldc )
798 END IF
799 END IF
800
801
802
803 GO TO 60
804
805
806
807 160 CONTINUE
808 DO 170 i = 1, n
809 IF( d( i ).LT.zero ) THEN
810 d( i ) = -d( i )
811
812
813
814 IF( ncvt.GT.0 )
815 $
CALL dscal( ncvt, negone, vt( i, 1 ), ldvt )
816 END IF
817 170 CONTINUE
818
819
820
821
822 DO 190 i = 1, n - 1
823
824
825
826 isub = 1
827 smin = d( 1 )
828 DO 180 j = 2, n + 1 - i
829 IF( d( j ).LE.smin ) THEN
830 isub = j
831 smin = d( j )
832 END IF
833 180 CONTINUE
834 IF( isub.NE.n+1-i ) THEN
835
836
837
838 d( isub ) = d( n+1-i )
839 d( n+1-i ) = smin
840 IF( ncvt.GT.0 )
841 $
CALL dswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
842 $ ldvt )
843 IF( nru.GT.0 )
844 $
CALL dswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
845 IF( ncc.GT.0 )
846 $
CALL dswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
847 END IF
848 190 CONTINUE
849 GO TO 220
850
851
852
853 200 CONTINUE
854 info = 0
855 DO 210 i = 1, n - 1
856 IF( e( i ).NE.zero )
857 $ info = info + 1
858 210 CONTINUE
859 220 CONTINUE
860 RETURN
861
862
863
double precision function dlamch(CMACH)
DLAMCH
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine dlasq1(N, D, E, WORK, INFO)
DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP