236 SUBROUTINE zgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
237 $ eps, sfmin, tol, nsweep, work, lwork, info )
246 DOUBLE PRECISION EPS, SFMIN, TOL
247 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
251 COMPLEX*16 A( lda, * ), D( n ), V( ldv, * ), WORK( lwork )
252 DOUBLE PRECISION SVA( n )
258 DOUBLE PRECISION ZERO, HALF, ONE
259 parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
262 COMPLEX*16 AAPQ, OMPQ
263 DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
264 $ bigtheta, cs, mxaapq, mxsinj, rootbig,
265 $ rooteps, rootsfmin, roottol, small, sn, t,
266 $ temp1, theta, thsign
267 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
268 $ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
269 $ p, pskipped, q, rowskip, swband
270 LOGICAL APPLV, ROTOK, RSVEC
274 INTRINSIC abs, conjg, max, dble, min, sign, sqrt
277 DOUBLE PRECISION DZNRM2
281 EXTERNAL idamax, lsame, zdotc, dznrm2
293 applv = lsame( jobv,
'A' )
294 rsvec = lsame( jobv,
'V' )
295 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv,
'N' ) ) )
THEN
297 ELSE IF( m.LT.0 )
THEN
299 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
301 ELSE IF( n1.LT.0 )
THEN
303 ELSE IF( lda.LT.m )
THEN
305 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
307 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
308 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
310 ELSE IF( tol.LE.eps )
THEN
312 ELSE IF( nsweep.LT.0 )
THEN
314 ELSE IF( lwork.LT.m )
THEN
322 CALL xerbla(
'ZGSVJ1', -info )
328 ELSE IF( applv )
THEN
331 rsvec = rsvec .OR. applv
333 rooteps = sqrt( eps )
334 rootsfmin = sqrt( sfmin )
337 rootbig = one / rootsfmin
339 bigtheta = one / rooteps
340 roottol = sqrt( tol )
353 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
357 nblc = ( n-n1 ) / kbl
358 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
359 blskip = ( kbl**2 ) + 1
362 rowskip = min( 5, kbl )
378 DO 1993 i = 1, nsweep
394 DO 2000 ibr = 1, nblr
396 igl = ( ibr-1 )*kbl + 1
402 igl = ( ibr-1 )*kbl + 1
405 DO 2010 jbc = 1, nblc
407 jgl = ( jbc-1 )*kbl + n1 + 1
412 DO 2100 p = igl, min( igl+kbl-1, n1 )
415 IF( aapp.GT.zero )
THEN
419 DO 2200 q = jgl, min( jgl+kbl-1, n )
422 IF( aaqq.GT.zero )
THEN
429 IF( aaqq.GE.one )
THEN
430 IF( aapp.GE.aaqq )
THEN
431 rotok = ( small*aapp ).LE.aaqq
433 rotok = ( small*aaqq ).LE.aapp
435 IF( aapp.LT.( big / aaqq ) )
THEN
436 aapq = ( zdotc( m, a( 1, p ), 1,
437 $ a( 1, q ), 1 ) / aaqq ) / aapp
439 CALL zcopy( m, a( 1, p ), 1,
441 CALL zlascl(
'G', 0, 0, aapp,
444 aapq = zdotc( m, work, 1,
445 $ a( 1, q ), 1 ) / aaqq
448 IF( aapp.GE.aaqq )
THEN
449 rotok = aapp.LE.( aaqq / small )
451 rotok = aaqq.LE.( aapp / small )
453 IF( aapp.GT.( small / aaqq ) )
THEN
454 aapq = ( zdotc( m, a( 1, p ), 1,
455 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
458 CALL zcopy( m, a( 1, q ), 1,
460 CALL zlascl(
'G', 0, 0, aaqq,
463 aapq = zdotc( m, a( 1, p ), 1,
470 mxaapq = max( mxaapq, -aapq1 )
474 IF( abs( aapq1 ).GT.tol )
THEN
475 ompq = aapq / abs(aapq)
485 theta = -half*abs( aqoap-apoaq )/ aapq1
486 IF( aaqq.GT.aapp0 )theta = -theta
488 IF( abs( theta ).GT.bigtheta )
THEN
491 CALL zrot( m, a(1,p), 1, a(1,q), 1,
492 $ cs, conjg(ompq)*t )
494 CALL zrot( mvl, v(1,p), 1,
495 $ v(1,q), 1, cs, conjg(ompq)*t )
497 sva( q ) = aaqq*sqrt( max( zero,
498 $ one+t*apoaq*aapq1 ) )
499 aapp = aapp*sqrt( max( zero,
500 $ one-t*aqoap*aapq1 ) )
501 mxsinj = max( mxsinj, abs( t ) )
506 thsign = -sign( one, aapq1 )
507 IF( aaqq.GT.aapp0 )thsign = -thsign
508 t = one / ( theta+thsign*
509 $ sqrt( one+theta*theta ) )
510 cs = sqrt( one / ( one+t*t ) )
512 mxsinj = max( mxsinj, abs( sn ) )
513 sva( q ) = aaqq*sqrt( max( zero,
514 $ one+t*apoaq*aapq1 ) )
515 aapp = aapp*sqrt( max( zero,
516 $ one-t*aqoap*aapq1 ) )
518 CALL zrot( m, a(1,p), 1, a(1,q), 1,
519 $ cs, conjg(ompq)*sn )
521 CALL zrot( mvl, v(1,p), 1,
522 $ v(1,q), 1, cs, conjg(ompq)*sn )
529 IF( aapp.GT.aaqq )
THEN
530 CALL zcopy( m, a( 1, p ), 1,
532 CALL zlascl(
'G', 0, 0, aapp, one,
535 CALL zlascl(
'G', 0, 0, aaqq, one,
536 $ m, 1, a( 1, q ), lda,
538 CALL zaxpy( m, -aapq, work,
540 CALL zlascl(
'G', 0, 0, one, aaqq,
541 $ m, 1, a( 1, q ), lda,
543 sva( q ) = aaqq*sqrt( max( zero,
544 $ one-aapq1*aapq1 ) )
545 mxsinj = max( mxsinj, sfmin )
547 CALL zcopy( m, a( 1, q ), 1,
549 CALL zlascl(
'G', 0, 0, aaqq, one,
552 CALL zlascl(
'G', 0, 0, aapp, one,
553 $ m, 1, a( 1, p ), lda,
555 CALL zaxpy( m, -conjg(aapq),
556 $ work, 1, a( 1, p ), 1 )
557 CALL zlascl(
'G', 0, 0, one, aapp,
558 $ m, 1, a( 1, p ), lda,
560 sva( p ) = aapp*sqrt( max( zero,
561 $ one-aapq1*aapq1 ) )
562 mxsinj = max( mxsinj, sfmin )
569 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
571 IF( ( aaqq.LT.rootbig ) .AND.
572 $ ( aaqq.GT.rootsfmin ) )
THEN
573 sva( q ) = dznrm2( m, a( 1, q ), 1)
577 CALL zlassq( m, a( 1, q ), 1, t,
579 sva( q ) = t*sqrt( aaqq )
582 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
583 IF( ( aapp.LT.rootbig ) .AND.
584 $ ( aapp.GT.rootsfmin ) )
THEN
585 aapp = dznrm2( m, a( 1, p ), 1 )
589 CALL zlassq( m, a( 1, p ), 1, t,
591 aapp = t*sqrt( aapp )
599 pskipped = pskipped + 1
604 pskipped = pskipped + 1
608 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
614 IF( ( i.LE.swband ) .AND.
615 $ ( pskipped.GT.rowskip ) )
THEN
629 IF( aapp.EQ.zero )notrot = notrot +
630 $ min( jgl+kbl-1, n ) - jgl + 1
631 IF( aapp.LT.zero )notrot = 0
641 DO 2012 p = igl, min( igl+kbl-1, n )
642 sva( p ) = abs( sva( p ) )
649 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
651 sva( n ) = dznrm2( m, a( 1, n ), 1 )
655 CALL zlassq( m, a( 1, n ), 1, t, aapp )
656 sva( n ) = t*sqrt( aapp )
661 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
662 $ ( iswrot.LE.n ) ) )swband = i
664 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( dble( n ) )*
665 $ tol ) .AND. ( dble( n )*mxaapq*mxsinj.LT.tol ) )
THEN
669 IF( notrot.GE.emptsw )
GO TO 1994
688 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
696 CALL zswap( m, a( 1, p ), 1, a( 1, q ), 1 )
697 IF( rsvec )
CALL zswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
ZGSVJ1 pre-processor for the routine zgesvj, applies Jacobi rotations targeting only particular pivot...
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY