236 SUBROUTINE cgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
237 $ eps, sfmin, tol, nsweep, work, lwork, info )
246 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
250 COMPLEX A( lda, * ), D( n ), V( ldv, * ), WORK( lwork )
258 parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
262 REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
263 $ bigtheta, cs, mxaapq, mxsinj, rootbig,
264 $ rooteps, rootsfmin, roottol, small, sn, t,
265 $ temp1, theta, thsign
266 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
267 $ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
268 $ p, pskipped, q, rowskip, swband
269 LOGICAL APPLV, ROTOK, RSVEC
273 INTRINSIC abs, max, conjg,
REAL, MIN, SIGN, SQRT
280 EXTERNAL isamax, lsame, cdotc, scnrm2
292 applv = lsame( jobv,
'A' )
293 rsvec = lsame( jobv,
'V' )
294 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv,
'N' ) ) )
THEN
296 ELSE IF( m.LT.0 )
THEN
298 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
300 ELSE IF( n1.LT.0 )
THEN
302 ELSE IF( lda.LT.m )
THEN
304 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
306 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
307 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
309 ELSE IF( tol.LE.eps )
THEN
311 ELSE IF( nsweep.LT.0 )
THEN
313 ELSE IF( lwork.LT.m )
THEN
321 CALL xerbla(
'CGSVJ1', -info )
327 ELSE IF( applv )
THEN
330 rsvec = rsvec .OR. applv
332 rooteps = sqrt( eps )
333 rootsfmin = sqrt( sfmin )
336 rootbig = one / rootsfmin
338 bigtheta = one / rooteps
339 roottol = sqrt( tol )
352 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
356 nblc = ( n-n1 ) / kbl
357 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
358 blskip = ( kbl**2 ) + 1
361 rowskip = min( 5, kbl )
377 DO 1993 i = 1, nsweep
393 DO 2000 ibr = 1, nblr
395 igl = ( ibr-1 )*kbl + 1
401 igl = ( ibr-1 )*kbl + 1
404 DO 2010 jbc = 1, nblc
406 jgl = ( jbc-1 )*kbl + n1 + 1
411 DO 2100 p = igl, min( igl+kbl-1, n1 )
414 IF( aapp.GT.zero )
THEN
418 DO 2200 q = jgl, min( jgl+kbl-1, n )
421 IF( aaqq.GT.zero )
THEN
428 IF( aaqq.GE.one )
THEN
429 IF( aapp.GE.aaqq )
THEN
430 rotok = ( small*aapp ).LE.aaqq
432 rotok = ( small*aaqq ).LE.aapp
434 IF( aapp.LT.( big / aaqq ) )
THEN
435 aapq = ( cdotc( m, a( 1, p ), 1,
436 $ a( 1, q ), 1 ) / aaqq ) / aapp
438 CALL ccopy( m, a( 1, p ), 1,
440 CALL clascl(
'G', 0, 0, aapp,
443 aapq = cdotc( m, work, 1,
444 $ a( 1, q ), 1 ) / aaqq
447 IF( aapp.GE.aaqq )
THEN
448 rotok = aapp.LE.( aaqq / small )
450 rotok = aaqq.LE.( aapp / small )
452 IF( aapp.GT.( small / aaqq ) )
THEN
453 aapq = ( cdotc( m, a( 1, p ), 1,
454 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
457 CALL ccopy( m, a( 1, q ), 1,
459 CALL clascl(
'G', 0, 0, aaqq,
462 aapq = cdotc( m, a( 1, p ), 1,
469 mxaapq = max( mxaapq, -aapq1 )
473 IF( abs( aapq1 ).GT.tol )
THEN
474 ompq = aapq / abs(aapq)
484 theta = -half*abs( aqoap-apoaq )/ aapq1
485 IF( aaqq.GT.aapp0 )theta = -theta
487 IF( abs( theta ).GT.bigtheta )
THEN
490 CALL crot( m, a(1,p), 1, a(1,q), 1,
491 $ cs, conjg(ompq)*t )
493 CALL crot( mvl, v(1,p), 1,
494 $ v(1,q), 1, cs, conjg(ompq)*t )
496 sva( q ) = aaqq*sqrt( max( zero,
497 $ one+t*apoaq*aapq1 ) )
498 aapp = aapp*sqrt( max( zero,
499 $ one-t*aqoap*aapq1 ) )
500 mxsinj = max( mxsinj, abs( t ) )
505 thsign = -sign( one, aapq1 )
506 IF( aaqq.GT.aapp0 )thsign = -thsign
507 t = one / ( theta+thsign*
508 $ sqrt( one+theta*theta ) )
509 cs = sqrt( one / ( one+t*t ) )
511 mxsinj = max( mxsinj, abs( sn ) )
512 sva( q ) = aaqq*sqrt( max( zero,
513 $ one+t*apoaq*aapq1 ) )
514 aapp = aapp*sqrt( max( zero,
515 $ one-t*aqoap*aapq1 ) )
517 CALL crot( m, a(1,p), 1, a(1,q), 1,
518 $ cs, conjg(ompq)*sn )
520 CALL crot( mvl, v(1,p), 1,
521 $ v(1,q), 1, cs, conjg(ompq)*sn )
528 IF( aapp.GT.aaqq )
THEN
529 CALL ccopy( m, a( 1, p ), 1,
531 CALL clascl(
'G', 0, 0, aapp, one,
534 CALL clascl(
'G', 0, 0, aaqq, one,
535 $ m, 1, a( 1, q ), lda,
537 CALL caxpy( m, -aapq, work,
539 CALL clascl(
'G', 0, 0, one, aaqq,
540 $ m, 1, a( 1, q ), lda,
542 sva( q ) = aaqq*sqrt( max( zero,
543 $ one-aapq1*aapq1 ) )
544 mxsinj = max( mxsinj, sfmin )
546 CALL ccopy( m, a( 1, q ), 1,
548 CALL clascl(
'G', 0, 0, aaqq, one,
551 CALL clascl(
'G', 0, 0, aapp, one,
552 $ m, 1, a( 1, p ), lda,
554 CALL caxpy( m, -conjg(aapq),
555 $ work, 1, a( 1, p ), 1 )
556 CALL clascl(
'G', 0, 0, one, aapp,
557 $ m, 1, a( 1, p ), lda,
559 sva( p ) = aapp*sqrt( max( zero,
560 $ one-aapq1*aapq1 ) )
561 mxsinj = max( mxsinj, sfmin )
568 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
570 IF( ( aaqq.LT.rootbig ) .AND.
571 $ ( aaqq.GT.rootsfmin ) )
THEN
572 sva( q ) = scnrm2( m, a( 1, q ), 1)
576 CALL classq( m, a( 1, q ), 1, t,
578 sva( q ) = t*sqrt( aaqq )
581 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
582 IF( ( aapp.LT.rootbig ) .AND.
583 $ ( aapp.GT.rootsfmin ) )
THEN
584 aapp = scnrm2( m, a( 1, p ), 1 )
588 CALL classq( m, a( 1, p ), 1, t,
590 aapp = t*sqrt( aapp )
598 pskipped = pskipped + 1
603 pskipped = pskipped + 1
607 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
613 IF( ( i.LE.swband ) .AND.
614 $ ( pskipped.GT.rowskip ) )
THEN
628 IF( aapp.EQ.zero )notrot = notrot +
629 $ min( jgl+kbl-1, n ) - jgl + 1
630 IF( aapp.LT.zero )notrot = 0
640 DO 2012 p = igl, min( igl+kbl-1, n )
641 sva( p ) = abs( sva( p ) )
648 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
650 sva( n ) = scnrm2( m, a( 1, n ), 1 )
654 CALL classq( m, a( 1, n ), 1, t, aapp )
655 sva( n ) = t*sqrt( aapp )
660 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
661 $ ( iswrot.LE.n ) ) )swband = i
663 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt(
REAL( N ) )*
664 $ tol ) .AND. (
REAL( n )*MXAAPQ*MXSINJ.LT.TOL ) ) then
668 IF( notrot.GE.emptsw )
GO TO 1994
687 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
695 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
696 IF( rsvec )
CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
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
subroutine cgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
CGSVJ1 pre-processor for the routine cgesvj, applies Jacobi rotations targeting only particular pivot...