218 SUBROUTINE cgsvj0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
219 $ sfmin, tol, nsweep, work, lwork, info )
228 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
233 COMPLEX A( lda, * ), D( n ), V( ldv, * ), WORK( lwork )
241 parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
243 parameter( czero = (0.0e0, 0.0e0), cone = (1.0e0, 0.0e0) )
247 REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
248 $ bigtheta, cs, mxaapq, mxsinj, rootbig, rooteps,
249 $ rootsfmin, roottol, small, sn, t, temp1, theta,
251 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
252 $ iswrot, jbc, jgl, kbl, lkahead, mvl, nbl,
253 $ notrot, p, pskipped, q, rowskip, swband
254 LOGICAL APPLV, ROTOK, RSVEC
258 INTRINSIC abs, max, conjg,
REAL, MIN, SIGN, SQRT
265 EXTERNAL isamax, lsame, cdotc, scnrm2
279 applv = lsame( jobv,
'A' )
280 rsvec = lsame( jobv,
'V' )
281 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv,
'N' ) ) )
THEN
283 ELSE IF( m.LT.0 )
THEN
285 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
287 ELSE IF( lda.LT.m )
THEN
289 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
291 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
292 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
294 ELSE IF( tol.LE.eps )
THEN
296 ELSE IF( nsweep.LT.0 )
THEN
298 ELSE IF( lwork.LT.m )
THEN
306 CALL xerbla(
'CGSVJ0', -info )
312 ELSE IF( applv )
THEN
315 rsvec = rsvec .OR. applv
317 rooteps = sqrt( eps )
318 rootsfmin = sqrt( sfmin )
321 rootbig = one / rootsfmin
322 bigtheta = one / rooteps
323 roottol = sqrt( tol )
327 emptsw = ( n*( n-1 ) ) / 2
348 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
353 rowskip = min( 5, kbl )
367 DO 1993 i = 1, nsweep
385 igl = ( ibr-1 )*kbl + 1
387 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
391 DO 2001 p = igl, min( igl+kbl-1, n-1 )
395 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
397 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
398 IF( rsvec )
CALL cswap( mvl, v( 1, p ), 1,
422 IF( ( sva( p ).LT.rootbig ) .AND.
423 $ ( sva( p ).GT.rootsfmin ) )
THEN
424 sva( p ) = scnrm2( m, a( 1, p ), 1 )
428 CALL classq( m, a( 1, p ), 1, temp1, aapp )
429 sva( p ) = temp1*sqrt( aapp )
436 IF( aapp.GT.zero )
THEN
440 DO 2002 q = p + 1, min( igl+kbl-1, n )
444 IF( aaqq.GT.zero )
THEN
447 IF( aaqq.GE.one )
THEN
448 rotok = ( small*aapp ).LE.aaqq
449 IF( aapp.LT.( big / aaqq ) )
THEN
450 aapq = ( cdotc( m, a( 1, p ), 1,
451 $ a( 1, q ), 1 ) / aaqq ) / aapp
453 CALL ccopy( m, a( 1, p ), 1,
455 CALL clascl(
'G', 0, 0, aapp, one,
456 $ m, 1, work, lda, ierr )
457 aapq = cdotc( m, work, 1,
458 $ a( 1, q ), 1 ) / aaqq
461 rotok = aapp.LE.( aaqq / small )
462 IF( aapp.GT.( small / aaqq ) )
THEN
463 aapq = ( cdotc( m, a( 1, p ), 1,
464 $ a( 1, q ), 1 ) / aapp ) / aaqq
466 CALL ccopy( m, a( 1, q ), 1,
468 CALL clascl(
'G', 0, 0, aaqq,
471 aapq = cdotc( m, a( 1, p ), 1,
478 mxaapq = max( mxaapq, -aapq1 )
482 IF( abs( aapq1 ).GT.tol )
THEN
483 ompq = aapq / abs(aapq)
498 theta = -half*abs( aqoap-apoaq )/aapq1
500 IF( abs( theta ).GT.bigtheta )
THEN
505 CALL crot( m, a(1,p), 1, a(1,q), 1,
506 $ cs, conjg(ompq)*t )
508 CALL crot( mvl, v(1,p), 1,
509 $ v(1,q), 1, cs, conjg(ompq)*t )
512 sva( q ) = aaqq*sqrt( max( zero,
513 $ one+t*apoaq*aapq1 ) )
514 aapp = aapp*sqrt( max( zero,
515 $ one-t*aqoap*aapq1 ) )
516 mxsinj = max( mxsinj, abs( t ) )
522 thsign = -sign( one, aapq1 )
523 t = one / ( theta+thsign*
524 $ sqrt( one+theta*theta ) )
525 cs = sqrt( one / ( one+t*t ) )
528 mxsinj = max( mxsinj, abs( sn ) )
529 sva( q ) = aaqq*sqrt( max( zero,
530 $ one+t*apoaq*aapq1 ) )
531 aapp = aapp*sqrt( max( zero,
532 $ one-t*aqoap*aapq1 ) )
534 CALL crot( m, a(1,p), 1, a(1,q), 1,
535 $ cs, conjg(ompq)*sn )
537 CALL crot( mvl, v(1,p), 1,
538 $ v(1,q), 1, cs, conjg(ompq)*sn )
545 CALL ccopy( m, a( 1, p ), 1,
547 CALL clascl(
'G', 0, 0, aapp, one, m,
550 CALL clascl(
'G', 0, 0, aaqq, one, m,
551 $ 1, a( 1, q ), lda, ierr )
552 CALL caxpy( m, -aapq, work, 1,
554 CALL clascl(
'G', 0, 0, one, aaqq, m,
555 $ 1, a( 1, q ), lda, ierr )
556 sva( q ) = aaqq*sqrt( max( zero,
557 $ one-aapq1*aapq1 ) )
558 mxsinj = max( mxsinj, sfmin )
565 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
567 IF( ( aaqq.LT.rootbig ) .AND.
568 $ ( aaqq.GT.rootsfmin ) )
THEN
569 sva( q ) = scnrm2( m, a( 1, q ), 1 )
573 CALL classq( m, a( 1, q ), 1, t,
575 sva( q ) = t*sqrt( aaqq )
578 IF( ( aapp / aapp0 ).LE.rooteps )
THEN
579 IF( ( aapp.LT.rootbig ) .AND.
580 $ ( aapp.GT.rootsfmin ) )
THEN
581 aapp = scnrm2( m, a( 1, p ), 1 )
585 CALL classq( m, a( 1, p ), 1, t,
587 aapp = t*sqrt( aapp )
594 IF( ir1.EQ.0 )notrot = notrot + 1
596 pskipped = pskipped + 1
600 IF( ir1.EQ.0 )notrot = notrot + 1
601 pskipped = pskipped + 1
604 IF( ( i.LE.swband ) .AND.
605 $ ( pskipped.GT.rowskip ) )
THEN
606 IF( ir1.EQ.0 )aapp = -aapp
621 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
622 $ notrot = notrot + min( igl+kbl-1, n ) - p
633 igl = ( ibr-1 )*kbl + 1
635 DO 2010 jbc = ibr + 1, nbl
637 jgl = ( jbc-1 )*kbl + 1
642 DO 2100 p = igl, min( igl+kbl-1, n )
645 IF( aapp.GT.zero )
THEN
649 DO 2200 q = jgl, min( jgl+kbl-1, n )
652 IF( aaqq.GT.zero )
THEN
659 IF( aaqq.GE.one )
THEN
660 IF( aapp.GE.aaqq )
THEN
661 rotok = ( small*aapp ).LE.aaqq
663 rotok = ( small*aaqq ).LE.aapp
665 IF( aapp.LT.( big / aaqq ) )
THEN
666 aapq = ( cdotc( m, a( 1, p ), 1,
667 $ a( 1, q ), 1 ) / aaqq ) / aapp
669 CALL ccopy( m, a( 1, p ), 1,
671 CALL clascl(
'G', 0, 0, aapp,
674 aapq = cdotc( m, work, 1,
675 $ a( 1, q ), 1 ) / aaqq
678 IF( aapp.GE.aaqq )
THEN
679 rotok = aapp.LE.( aaqq / small )
681 rotok = aaqq.LE.( aapp / small )
683 IF( aapp.GT.( small / aaqq ) )
THEN
684 aapq = ( cdotc( m, a( 1, p ), 1,
685 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
688 CALL ccopy( m, a( 1, q ), 1,
690 CALL clascl(
'G', 0, 0, aaqq,
693 aapq = cdotc( m, a( 1, p ), 1,
700 mxaapq = max( mxaapq, -aapq1 )
704 IF( abs( aapq1 ).GT.tol )
THEN
705 ompq = aapq / abs(aapq)
715 theta = -half*abs( aqoap-apoaq )/ aapq1
716 IF( aaqq.GT.aapp0 )theta = -theta
718 IF( abs( theta ).GT.bigtheta )
THEN
721 CALL crot( m, a(1,p), 1, a(1,q), 1,
722 $ cs, conjg(ompq)*t )
724 CALL crot( mvl, v(1,p), 1,
725 $ v(1,q), 1, cs, conjg(ompq)*t )
727 sva( q ) = aaqq*sqrt( max( zero,
728 $ one+t*apoaq*aapq1 ) )
729 aapp = aapp*sqrt( max( zero,
730 $ one-t*aqoap*aapq1 ) )
731 mxsinj = max( mxsinj, abs( t ) )
736 thsign = -sign( one, aapq1 )
737 IF( aaqq.GT.aapp0 )thsign = -thsign
738 t = one / ( theta+thsign*
739 $ sqrt( one+theta*theta ) )
740 cs = sqrt( one / ( one+t*t ) )
742 mxsinj = max( mxsinj, abs( sn ) )
743 sva( q ) = aaqq*sqrt( max( zero,
744 $ one+t*apoaq*aapq1 ) )
745 aapp = aapp*sqrt( max( zero,
746 $ one-t*aqoap*aapq1 ) )
748 CALL crot( m, a(1,p), 1, a(1,q), 1,
749 $ cs, conjg(ompq)*sn )
751 CALL crot( mvl, v(1,p), 1,
752 $ v(1,q), 1, cs, conjg(ompq)*sn )
759 IF( aapp.GT.aaqq )
THEN
760 CALL ccopy( m, a( 1, p ), 1,
762 CALL clascl(
'G', 0, 0, aapp, one,
765 CALL clascl(
'G', 0, 0, aaqq, one,
766 $ m, 1, a( 1, q ), lda,
768 CALL caxpy( m, -aapq, work,
770 CALL clascl(
'G', 0, 0, one, aaqq,
771 $ m, 1, a( 1, q ), lda,
773 sva( q ) = aaqq*sqrt( max( zero,
774 $ one-aapq1*aapq1 ) )
775 mxsinj = max( mxsinj, sfmin )
777 CALL ccopy( m, a( 1, q ), 1,
779 CALL clascl(
'G', 0, 0, aaqq, one,
782 CALL clascl(
'G', 0, 0, aapp, one,
783 $ m, 1, a( 1, p ), lda,
785 CALL caxpy( m, -conjg(aapq),
786 $ work, 1, a( 1, p ), 1 )
787 CALL clascl(
'G', 0, 0, one, aapp,
788 $ m, 1, a( 1, p ), lda,
790 sva( p ) = aapp*sqrt( max( zero,
791 $ one-aapq1*aapq1 ) )
792 mxsinj = max( mxsinj, sfmin )
799 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
801 IF( ( aaqq.LT.rootbig ) .AND.
802 $ ( aaqq.GT.rootsfmin ) )
THEN
803 sva( q ) = scnrm2( m, a( 1, q ), 1)
807 CALL classq( m, a( 1, q ), 1, t,
809 sva( q ) = t*sqrt( aaqq )
812 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
813 IF( ( aapp.LT.rootbig ) .AND.
814 $ ( aapp.GT.rootsfmin ) )
THEN
815 aapp = scnrm2( m, a( 1, p ), 1 )
819 CALL classq( m, a( 1, p ), 1, t,
821 aapp = t*sqrt( aapp )
829 pskipped = pskipped + 1
834 pskipped = pskipped + 1
838 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
844 IF( ( i.LE.swband ) .AND.
845 $ ( pskipped.GT.rowskip ) )
THEN
859 IF( aapp.EQ.zero )notrot = notrot +
860 $ min( jgl+kbl-1, n ) - jgl + 1
861 IF( aapp.LT.zero )notrot = 0
871 DO 2012 p = igl, min( igl+kbl-1, n )
872 sva( p ) = abs( sva( p ) )
879 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
881 sva( n ) = scnrm2( m, a( 1, n ), 1 )
885 CALL classq( m, a( 1, n ), 1, t, aapp )
886 sva( n ) = t*sqrt( aapp )
891 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
892 $ ( iswrot.LE.n ) ) )swband = i
894 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt(
REAL( N ) )*
895 $ tol ) .AND. (
REAL( n )*MXAAPQ*MXSINJ.LT.TOL ) ) then
899 IF( notrot.GE.emptsw )
GO TO 1994
918 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
926 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
927 IF( rsvec )
CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
subroutine cgsvj0(JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
CGSVJ0 pre-processor for the routine cgesvj.
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
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 xerbla(SRNAME, INFO)
XERBLA
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY