LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zgsvj0()

subroutine zgsvj0 ( character*1 jobv,
integer m,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( n ) d,
double precision, dimension( n ) sva,
integer mv,
complex*16, dimension( ldv, * ) v,
integer ldv,
double precision eps,
double precision sfmin,
double precision tol,
integer nsweep,
complex*16, dimension( lwork ) work,
integer lwork,
integer info )

ZGSVJ0 pre-processor for the routine zgesvj.

Download ZGSVJ0 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZGSVJ0 is called from ZGESVJ as a pre-processor and that is its main
!> purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but
!> it does not check convergence (stopping criterion). Few tuning
!> parameters (marked by [TP]) are available for the implementer.
!>
Parameters
[in]JOBV
!>          JOBV is CHARACTER*1
!>          Specifies whether the output from this procedure is used
!>          to compute the matrix V:
!>          = 'V': the product of the Jacobi rotations is accumulated
!>                 by postmultiplying the N-by-N array V.
!>                (See the description of V.)
!>          = 'A': the product of the Jacobi rotations is accumulated
!>                 by postmultiplying the MV-by-N array V.
!>                (See the descriptions of MV and V.)
!>          = 'N': the Jacobi rotations are not accumulated.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the input matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the input matrix A.
!>          M >= N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, M-by-N matrix A, such that A*diag(D) represents
!>          the input matrix.
!>          On exit,
!>          A_onexit * diag(D_onexit) represents the input matrix A*diag(D)
!>          post-multiplied by a sequence of Jacobi rotations, where the
!>          rotation threshold and the total number of sweeps are given in
!>          TOL and NSWEEP, respectively.
!>          (See the descriptions of D, TOL and NSWEEP.)
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]D
!>          D is COMPLEX*16 array, dimension (N)
!>          The array D accumulates the scaling factors from the complex scaled
!>          Jacobi rotations.
!>          On entry, A*diag(D) represents the input matrix.
!>          On exit, A_onexit*diag(D_onexit) represents the input matrix
!>          post-multiplied by a sequence of Jacobi rotations, where the
!>          rotation threshold and the total number of sweeps are given in
!>          TOL and NSWEEP, respectively.
!>          (See the descriptions of A, TOL and NSWEEP.)
!>
[in,out]SVA
!>          SVA is DOUBLE PRECISION array, dimension (N)
!>          On entry, SVA contains the Euclidean norms of the columns of
!>          the matrix A*diag(D).
!>          On exit, SVA contains the Euclidean norms of the columns of
!>          the matrix A_onexit*diag(D_onexit).
!>
[in]MV
!>          MV is INTEGER
!>          If JOBV = 'A', then MV rows of V are post-multiplied by a
!>                           sequence of Jacobi rotations.
!>          If JOBV = 'N',   then MV is not referenced.
!>
[in,out]V
!>          V is COMPLEX*16 array, dimension (LDV,N)
!>          If JOBV = 'V' then N rows of V are post-multiplied by a
!>                           sequence of Jacobi rotations.
!>          If JOBV = 'A' then MV rows of V are post-multiplied by a
!>                           sequence of Jacobi rotations.
!>          If JOBV = 'N',   then V is not referenced.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V,  LDV >= 1.
!>          If JOBV = 'V', LDV >= N.
!>          If JOBV = 'A', LDV >= MV.
!>
[in]EPS
!>          EPS is DOUBLE PRECISION
!>          EPS = DLAMCH('Epsilon')
!>
[in]SFMIN
!>          SFMIN is DOUBLE PRECISION
!>          SFMIN = DLAMCH('Safe Minimum')
!>
[in]TOL
!>          TOL is DOUBLE PRECISION
!>          TOL is the threshold for Jacobi rotations. For a pair
!>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
!>          applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.
!>
[in]NSWEEP
!>          NSWEEP is INTEGER
!>          NSWEEP is the number of sweeps of Jacobi rotations to be
!>          performed.
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          LWORK is the dimension of WORK. LWORK >= M.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, then the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
ZGSVJ0 is used just to enable ZGESVJ to call a simplified version of itself to work on a submatrix of the original matrix.

Contributor: Zlatko Drmac (Zagreb, Croatia)

Bugs, Examples and Comments:
Please report all bugs and send interesting test examples and comments to drmac.nosp@m.@mat.nosp@m.h.hr. Thank you.

Definition at line 214 of file zgsvj0.f.

216*
217* -- LAPACK computational routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221 IMPLICIT NONE
222* .. Scalar Arguments ..
223 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
224 DOUBLE PRECISION EPS, SFMIN, TOL
225 CHARACTER*1 JOBV
226* ..
227* .. Array Arguments ..
228 COMPLEX*16 A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
229 DOUBLE PRECISION SVA( N )
230* ..
231*
232* =====================================================================
233*
234* .. Local Parameters ..
235 DOUBLE PRECISION ZERO, HALF, ONE
236 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
237 COMPLEX*16 CZERO, CONE
238 parameter( czero = (0.0d0, 0.0d0), cone = (1.0d0, 0.0d0) )
239* ..
240* .. Local Scalars ..
241 COMPLEX*16 AAPQ, OMPQ
242 DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
243 $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
244 $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
245 $ THSIGN
246 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
247 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
248 $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
249 LOGICAL APPLV, ROTOK, RSVEC
250* ..
251* ..
252* .. Intrinsic Functions ..
253 INTRINSIC abs, max, conjg, dble, min, sign, sqrt
254* ..
255* .. External Functions ..
256 DOUBLE PRECISION DZNRM2
257 COMPLEX*16 ZDOTC
258 INTEGER IDAMAX
259 LOGICAL LSAME
260 EXTERNAL idamax, lsame, zdotc, dznrm2
261* ..
262* ..
263* .. External Subroutines ..
264* ..
265* from BLAS
266 EXTERNAL zcopy, zrot, zswap, zaxpy
267* from LAPACK
268 EXTERNAL zlascl, zlassq, xerbla
269* ..
270* .. Executable Statements ..
271*
272* Test the input parameters.
273*
274 applv = lsame( jobv, 'A' )
275 rsvec = lsame( jobv, 'V' )
276 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
277 info = -1
278 ELSE IF( m.LT.0 ) THEN
279 info = -2
280 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
281 info = -3
282 ELSE IF( lda.LT.m ) THEN
283 info = -5
284 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
285 info = -8
286 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
287 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
288 info = -10
289 ELSE IF( tol.LE.eps ) THEN
290 info = -13
291 ELSE IF( nsweep.LT.0 ) THEN
292 info = -14
293 ELSE IF( lwork.LT.m ) THEN
294 info = -16
295 ELSE
296 info = 0
297 END IF
298*
299* #:(
300 IF( info.NE.0 ) THEN
301 CALL xerbla( 'ZGSVJ0', -info )
302 RETURN
303 END IF
304*
305 IF( rsvec ) THEN
306 mvl = n
307 ELSE IF( applv ) THEN
308 mvl = mv
309 END IF
310 rsvec = rsvec .OR. applv
311
312 rooteps = sqrt( eps )
313 rootsfmin = sqrt( sfmin )
314 small = sfmin / eps
315 big = one / sfmin
316 rootbig = one / rootsfmin
317 bigtheta = one / rooteps
318 roottol = sqrt( tol )
319*
320* .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
321*
322 emptsw = ( n*( n-1 ) ) / 2
323 notrot = 0
324*
325* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
326*
327
328 swband = 0
329*[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
330* if ZGESVJ is used as a computational routine in the preconditioned
331* Jacobi SVD algorithm ZGEJSV. For sweeps i=1:SWBAND the procedure
332* works on pivots inside a band-like region around the diagonal.
333* The boundaries are determined dynamically, based on the number of
334* pivots above a threshold.
335*
336 kbl = min( 8, n )
337*[TP] KBL is a tuning parameter that defines the tile size in the
338* tiling of the p-q loops of pivot pairs. In general, an optimal
339* value of KBL depends on the matrix dimensions and on the
340* parameters of the computer's memory.
341*
342 nbl = n / kbl
343 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
344*
345 blskip = kbl**2
346*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
347*
348 rowskip = min( 5, kbl )
349*[TP] ROWSKIP is a tuning parameter.
350*
351 lkahead = 1
352*[TP] LKAHEAD is a tuning parameter.
353*
354* Quasi block transformations, using the lower (upper) triangular
355* structure of the input matrix. The quasi-block-cycling usually
356* invokes cubic convergence. Big part of this cycle is done inside
357* canonical subspaces of dimensions less than M.
358*
359*
360* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
361*
362 DO 1993 i = 1, nsweep
363*
364* .. go go go ...
365*
366 mxaapq = zero
367 mxsinj = zero
368 iswrot = 0
369*
370 notrot = 0
371 pskipped = 0
372*
373* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
374* 1 <= p < q <= N. This is the first step toward a blocked implementation
375* of the rotations. New implementation, based on block transformations,
376* is under development.
377*
378 DO 2000 ibr = 1, nbl
379*
380 igl = ( ibr-1 )*kbl + 1
381*
382 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
383*
384 igl = igl + ir1*kbl
385*
386 DO 2001 p = igl, min( igl+kbl-1, n-1 )
387*
388* .. de Rijk's pivoting
389*
390 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
391 IF( p.NE.q ) THEN
392 CALL zswap( m, a( 1, p ), 1, a( 1, q ), 1 )
393 IF( rsvec )CALL zswap( mvl, v( 1, p ), 1,
394 $ v( 1, q ), 1 )
395 temp1 = sva( p )
396 sva( p ) = sva( q )
397 sva( q ) = temp1
398 aapq = d(p)
399 d(p) = d(q)
400 d(q) = aapq
401 END IF
402*
403 IF( ir1.EQ.0 ) THEN
404*
405* Column norms are periodically updated by explicit
406* norm computation.
407* Caveat:
408* Unfortunately, some BLAS implementations compute SNCRM2(M,A(1,p),1)
409* as SQRT(S=ZDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
410* overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
411* underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
412* Hence, DZNRM2 cannot be trusted, not even in the case when
413* the true norm is far from the under(over)flow boundaries.
414* If properly implemented DZNRM2 is available, the IF-THEN-ELSE-END IF
415* below should be replaced with "AAPP = DZNRM2( M, A(1,p), 1 )".
416*
417 IF( ( sva( p ).LT.rootbig ) .AND.
418 $ ( sva( p ).GT.rootsfmin ) ) THEN
419 sva( p ) = dznrm2( m, a( 1, p ), 1 )
420 ELSE
421 temp1 = zero
422 aapp = one
423 CALL zlassq( m, a( 1, p ), 1, temp1, aapp )
424 sva( p ) = temp1*sqrt( aapp )
425 END IF
426 aapp = sva( p )
427 ELSE
428 aapp = sva( p )
429 END IF
430*
431 IF( aapp.GT.zero ) THEN
432*
433 pskipped = 0
434*
435 DO 2002 q = p + 1, min( igl+kbl-1, n )
436*
437 aaqq = sva( q )
438*
439 IF( aaqq.GT.zero ) THEN
440*
441 aapp0 = aapp
442 IF( aaqq.GE.one ) THEN
443 rotok = ( small*aapp ).LE.aaqq
444 IF( aapp.LT.( big / aaqq ) ) THEN
445 aapq = ( zdotc( m, a( 1, p ), 1,
446 $ a( 1, q ), 1 ) / aaqq ) / aapp
447 ELSE
448 CALL zcopy( m, a( 1, p ), 1,
449 $ work, 1 )
450 CALL zlascl( 'G', 0, 0, aapp, one,
451 $ m, 1, work, lda, ierr )
452 aapq = zdotc( m, work, 1,
453 $ a( 1, q ), 1 ) / aaqq
454 END IF
455 ELSE
456 rotok = aapp.LE.( aaqq / small )
457 IF( aapp.GT.( small / aaqq ) ) THEN
458 aapq = ( zdotc( m, a( 1, p ), 1,
459 $ a( 1, q ), 1 ) / aapp ) / aaqq
460 ELSE
461 CALL zcopy( m, a( 1, q ), 1,
462 $ work, 1 )
463 CALL zlascl( 'G', 0, 0, aaqq,
464 $ one, m, 1,
465 $ work, lda, ierr )
466 aapq = zdotc( m, a( 1, p ), 1,
467 $ work, 1 ) / aapp
468 END IF
469 END IF
470*
471* AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
472 aapq1 = -abs(aapq)
473 mxaapq = max( mxaapq, -aapq1 )
474*
475* TO rotate or NOT to rotate, THAT is the question ...
476*
477 IF( abs( aapq1 ).GT.tol ) THEN
478 ompq = aapq / abs(aapq)
479*
480* .. rotate
481*[RTD] ROTATED = ROTATED + ONE
482*
483 IF( ir1.EQ.0 ) THEN
484 notrot = 0
485 pskipped = 0
486 iswrot = iswrot + 1
487 END IF
488*
489 IF( rotok ) THEN
490*
491 aqoap = aaqq / aapp
492 apoaq = aapp / aaqq
493 theta = -half*abs( aqoap-apoaq )/aapq1
494*
495 IF( abs( theta ).GT.bigtheta ) THEN
496*
497 t = half / theta
498 cs = one
499
500 CALL zrot( m, a(1,p), 1, a(1,q),
501 $ 1,
502 $ cs, conjg(ompq)*t )
503 IF ( rsvec ) THEN
504 CALL zrot( mvl, v(1,p), 1,
505 $ v(1,q), 1, cs, conjg(ompq)*t )
506 END IF
507
508 sva( q ) = aaqq*sqrt( max( zero,
509 $ one+t*apoaq*aapq1 ) )
510 aapp = aapp*sqrt( max( zero,
511 $ one-t*aqoap*aapq1 ) )
512 mxsinj = max( mxsinj, abs( t ) )
513*
514 ELSE
515*
516* .. choose correct signum for THETA and rotate
517*
518 thsign = -sign( one, aapq1 )
519 t = one / ( theta+thsign*
520 $ sqrt( one+theta*theta ) )
521 cs = sqrt( one / ( one+t*t ) )
522 sn = t*cs
523*
524 mxsinj = max( mxsinj, abs( sn ) )
525 sva( q ) = aaqq*sqrt( max( zero,
526 $ one+t*apoaq*aapq1 ) )
527 aapp = aapp*sqrt( max( zero,
528 $ one-t*aqoap*aapq1 ) )
529*
530 CALL zrot( m, a(1,p), 1, a(1,q),
531 $ 1,
532 $ cs, conjg(ompq)*sn )
533 IF ( rsvec ) THEN
534 CALL zrot( mvl, v(1,p), 1,
535 $ v(1,q), 1, cs, conjg(ompq)*sn )
536 END IF
537 END IF
538 d(p) = -d(q) * ompq
539*
540 ELSE
541* .. have to use modified Gram-Schmidt like transformation
542 CALL zcopy( m, a( 1, p ), 1,
543 $ work, 1 )
544 CALL zlascl( 'G', 0, 0, aapp, one,
545 $ m,
546 $ 1, work, lda,
547 $ ierr )
548 CALL zlascl( 'G', 0, 0, aaqq, one,
549 $ m,
550 $ 1, a( 1, q ), lda, ierr )
551 CALL zaxpy( m, -aapq, work, 1,
552 $ a( 1, q ), 1 )
553 CALL zlascl( 'G', 0, 0, one, aaqq,
554 $ m,
555 $ 1, a( 1, q ), lda, ierr )
556 sva( q ) = aaqq*sqrt( max( zero,
557 $ one-aapq1*aapq1 ) )
558 mxsinj = max( mxsinj, sfmin )
559 END IF
560* END IF ROTOK THEN ... ELSE
561*
562* In the case of cancellation in updating SVA(q), SVA(p)
563* recompute SVA(q), SVA(p).
564*
565 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
566 $ THEN
567 IF( ( aaqq.LT.rootbig ) .AND.
568 $ ( aaqq.GT.rootsfmin ) ) THEN
569 sva( q ) = dznrm2( m, a( 1, q ),
570 $ 1 )
571 ELSE
572 t = zero
573 aaqq = one
574 CALL zlassq( m, a( 1, q ), 1, t,
575 $ aaqq )
576 sva( q ) = t*sqrt( aaqq )
577 END IF
578 END IF
579 IF( ( aapp / aapp0 ).LE.rooteps ) THEN
580 IF( ( aapp.LT.rootbig ) .AND.
581 $ ( aapp.GT.rootsfmin ) ) THEN
582 aapp = dznrm2( m, a( 1, p ), 1 )
583 ELSE
584 t = zero
585 aapp = one
586 CALL zlassq( m, a( 1, p ), 1, t,
587 $ aapp )
588 aapp = t*sqrt( aapp )
589 END IF
590 sva( p ) = aapp
591 END IF
592*
593 ELSE
594* A(:,p) and A(:,q) already numerically orthogonal
595 IF( ir1.EQ.0 )notrot = notrot + 1
596*[RTD] SKIPPED = SKIPPED + 1
597 pskipped = pskipped + 1
598 END IF
599 ELSE
600* A(:,q) is zero column
601 IF( ir1.EQ.0 )notrot = notrot + 1
602 pskipped = pskipped + 1
603 END IF
604*
605 IF( ( i.LE.swband ) .AND.
606 $ ( pskipped.GT.rowskip ) ) THEN
607 IF( ir1.EQ.0 )aapp = -aapp
608 notrot = 0
609 GO TO 2103
610 END IF
611*
612 2002 CONTINUE
613* END q-LOOP
614*
615 2103 CONTINUE
616* bailed out of q-loop
617*
618 sva( p ) = aapp
619*
620 ELSE
621 sva( p ) = aapp
622 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
623 $ notrot = notrot + min( igl+kbl-1, n ) - p
624 END IF
625*
626 2001 CONTINUE
627* end of the p-loop
628* end of doing the block ( ibr, ibr )
629 1002 CONTINUE
630* end of ir1-loop
631*
632* ... go to the off diagonal blocks
633*
634 igl = ( ibr-1 )*kbl + 1
635*
636 DO 2010 jbc = ibr + 1, nbl
637*
638 jgl = ( jbc-1 )*kbl + 1
639*
640* doing the block at ( ibr, jbc )
641*
642 ijblsk = 0
643 DO 2100 p = igl, min( igl+kbl-1, n )
644*
645 aapp = sva( p )
646 IF( aapp.GT.zero ) THEN
647*
648 pskipped = 0
649*
650 DO 2200 q = jgl, min( jgl+kbl-1, n )
651*
652 aaqq = sva( q )
653 IF( aaqq.GT.zero ) THEN
654 aapp0 = aapp
655*
656* .. M x 2 Jacobi SVD ..
657*
658* Safe Gram matrix computation
659*
660 IF( aaqq.GE.one ) THEN
661 IF( aapp.GE.aaqq ) THEN
662 rotok = ( small*aapp ).LE.aaqq
663 ELSE
664 rotok = ( small*aaqq ).LE.aapp
665 END IF
666 IF( aapp.LT.( big / aaqq ) ) THEN
667 aapq = ( zdotc( m, a( 1, p ), 1,
668 $ a( 1, q ), 1 ) / aaqq ) / aapp
669 ELSE
670 CALL zcopy( m, a( 1, p ), 1,
671 $ work, 1 )
672 CALL zlascl( 'G', 0, 0, aapp,
673 $ one, m, 1,
674 $ work, lda, ierr )
675 aapq = zdotc( m, work, 1,
676 $ a( 1, q ), 1 ) / aaqq
677 END IF
678 ELSE
679 IF( aapp.GE.aaqq ) THEN
680 rotok = aapp.LE.( aaqq / small )
681 ELSE
682 rotok = aaqq.LE.( aapp / small )
683 END IF
684 IF( aapp.GT.( small / aaqq ) ) THEN
685 aapq = ( zdotc( m, a( 1, p ), 1,
686 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
687 $ / min(aaqq,aapp)
688 ELSE
689 CALL zcopy( m, a( 1, q ), 1,
690 $ work, 1 )
691 CALL zlascl( 'G', 0, 0, aaqq,
692 $ one, m, 1,
693 $ work, lda, ierr )
694 aapq = zdotc( m, a( 1, p ), 1,
695 $ work, 1 ) / aapp
696 END IF
697 END IF
698*
699* AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
700 aapq1 = -abs(aapq)
701 mxaapq = max( mxaapq, -aapq1 )
702*
703* TO rotate or NOT to rotate, THAT is the question ...
704*
705 IF( abs( aapq1 ).GT.tol ) THEN
706 ompq = aapq / abs(aapq)
707 notrot = 0
708*[RTD] ROTATED = ROTATED + 1
709 pskipped = 0
710 iswrot = iswrot + 1
711*
712 IF( rotok ) THEN
713*
714 aqoap = aaqq / aapp
715 apoaq = aapp / aaqq
716 theta = -half*abs( aqoap-apoaq )/ aapq1
717 IF( aaqq.GT.aapp0 )theta = -theta
718*
719 IF( abs( theta ).GT.bigtheta ) THEN
720 t = half / theta
721 cs = one
722 CALL zrot( m, a(1,p), 1, a(1,q),
723 $ 1,
724 $ cs, conjg(ompq)*t )
725 IF( rsvec ) THEN
726 CALL zrot( mvl, v(1,p), 1,
727 $ v(1,q), 1, cs, conjg(ompq)*t )
728 END IF
729 sva( q ) = aaqq*sqrt( max( zero,
730 $ one+t*apoaq*aapq1 ) )
731 aapp = aapp*sqrt( max( zero,
732 $ one-t*aqoap*aapq1 ) )
733 mxsinj = max( mxsinj, abs( t ) )
734 ELSE
735*
736* .. choose correct signum for THETA and rotate
737*
738 thsign = -sign( one, aapq1 )
739 IF( aaqq.GT.aapp0 )thsign = -thsign
740 t = one / ( theta+thsign*
741 $ sqrt( one+theta*theta ) )
742 cs = sqrt( one / ( one+t*t ) )
743 sn = t*cs
744 mxsinj = max( mxsinj, abs( sn ) )
745 sva( q ) = aaqq*sqrt( max( zero,
746 $ one+t*apoaq*aapq1 ) )
747 aapp = aapp*sqrt( max( zero,
748 $ one-t*aqoap*aapq1 ) )
749*
750 CALL zrot( m, a(1,p), 1, a(1,q),
751 $ 1,
752 $ cs, conjg(ompq)*sn )
753 IF( rsvec ) THEN
754 CALL zrot( mvl, v(1,p), 1,
755 $ v(1,q), 1, cs, conjg(ompq)*sn )
756 END IF
757 END IF
758 d(p) = -d(q) * ompq
759*
760 ELSE
761* .. have to use modified Gram-Schmidt like transformation
762 IF( aapp.GT.aaqq ) THEN
763 CALL zcopy( m, a( 1, p ), 1,
764 $ work, 1 )
765 CALL zlascl( 'G', 0, 0, aapp,
766 $ one,
767 $ m, 1, work,lda,
768 $ ierr )
769 CALL zlascl( 'G', 0, 0, aaqq,
770 $ one,
771 $ m, 1, a( 1, q ), lda,
772 $ ierr )
773 CALL zaxpy( m, -aapq, work,
774 $ 1, a( 1, q ), 1 )
775 CALL zlascl( 'G', 0, 0, one,
776 $ aaqq,
777 $ m, 1, a( 1, q ), lda,
778 $ ierr )
779 sva( q ) = aaqq*sqrt( max( zero,
780 $ one-aapq1*aapq1 ) )
781 mxsinj = max( mxsinj, sfmin )
782 ELSE
783 CALL zcopy( m, a( 1, q ), 1,
784 $ work, 1 )
785 CALL zlascl( 'G', 0, 0, aaqq,
786 $ one,
787 $ m, 1, work,lda,
788 $ ierr )
789 CALL zlascl( 'G', 0, 0, aapp,
790 $ one,
791 $ m, 1, a( 1, p ), lda,
792 $ ierr )
793 CALL zaxpy( m, -conjg(aapq),
794 $ work, 1, a( 1, p ), 1 )
795 CALL zlascl( 'G', 0, 0, one,
796 $ aapp,
797 $ m, 1, a( 1, p ), lda,
798 $ ierr )
799 sva( p ) = aapp*sqrt( max( zero,
800 $ one-aapq1*aapq1 ) )
801 mxsinj = max( mxsinj, sfmin )
802 END IF
803 END IF
804* END IF ROTOK THEN ... ELSE
805*
806* In the case of cancellation in updating SVA(q), SVA(p)
807* .. recompute SVA(q), SVA(p)
808 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
809 $ THEN
810 IF( ( aaqq.LT.rootbig ) .AND.
811 $ ( aaqq.GT.rootsfmin ) ) THEN
812 sva( q ) = dznrm2( m, a( 1, q ),
813 $ 1)
814 ELSE
815 t = zero
816 aaqq = one
817 CALL zlassq( m, a( 1, q ), 1, t,
818 $ aaqq )
819 sva( q ) = t*sqrt( aaqq )
820 END IF
821 END IF
822 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
823 IF( ( aapp.LT.rootbig ) .AND.
824 $ ( aapp.GT.rootsfmin ) ) THEN
825 aapp = dznrm2( m, a( 1, p ), 1 )
826 ELSE
827 t = zero
828 aapp = one
829 CALL zlassq( m, a( 1, p ), 1, t,
830 $ aapp )
831 aapp = t*sqrt( aapp )
832 END IF
833 sva( p ) = aapp
834 END IF
835* end of OK rotation
836 ELSE
837 notrot = notrot + 1
838*[RTD] SKIPPED = SKIPPED + 1
839 pskipped = pskipped + 1
840 ijblsk = ijblsk + 1
841 END IF
842 ELSE
843 notrot = notrot + 1
844 pskipped = pskipped + 1
845 ijblsk = ijblsk + 1
846 END IF
847*
848 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
849 $ THEN
850 sva( p ) = aapp
851 notrot = 0
852 GO TO 2011
853 END IF
854 IF( ( i.LE.swband ) .AND.
855 $ ( pskipped.GT.rowskip ) ) THEN
856 aapp = -aapp
857 notrot = 0
858 GO TO 2203
859 END IF
860*
861 2200 CONTINUE
862* end of the q-loop
863 2203 CONTINUE
864*
865 sva( p ) = aapp
866*
867 ELSE
868*
869 IF( aapp.EQ.zero )notrot = notrot +
870 $ min( jgl+kbl-1, n ) - jgl + 1
871 IF( aapp.LT.zero )notrot = 0
872*
873 END IF
874*
875 2100 CONTINUE
876* end of the p-loop
877 2010 CONTINUE
878* end of the jbc-loop
879 2011 CONTINUE
880*2011 bailed out of the jbc-loop
881 DO 2012 p = igl, min( igl+kbl-1, n )
882 sva( p ) = abs( sva( p ) )
883 2012 CONTINUE
884***
885 2000 CONTINUE
886*2000 :: end of the ibr-loop
887*
888* .. update SVA(N)
889 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
890 $ THEN
891 sva( n ) = dznrm2( m, a( 1, n ), 1 )
892 ELSE
893 t = zero
894 aapp = one
895 CALL zlassq( m, a( 1, n ), 1, t, aapp )
896 sva( n ) = t*sqrt( aapp )
897 END IF
898*
899* Additional steering devices
900*
901 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
902 $ ( iswrot.LE.n ) ) )swband = i
903*
904 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( dble( n ) )*
905 $ tol ) .AND. ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
906 GO TO 1994
907 END IF
908*
909 IF( notrot.GE.emptsw )GO TO 1994
910*
911 1993 CONTINUE
912* end i=1:NSWEEP loop
913*
914* #:( Reaching this point means that the procedure has not converged.
915 info = nsweep - 1
916 GO TO 1995
917*
918 1994 CONTINUE
919* #:) Reaching this point means numerical convergence after the i-th
920* sweep.
921*
922 info = 0
923* #:) INFO = 0 confirms successful iterations.
924 1995 CONTINUE
925*
926* Sort the vector SVA() of column norms.
927 DO 5991 p = 1, n - 1
928 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
929 IF( p.NE.q ) THEN
930 temp1 = sva( p )
931 sva( p ) = sva( q )
932 sva( q ) = temp1
933 aapq = d( p )
934 d( p ) = d( q )
935 d( q ) = aapq
936 CALL zswap( m, a( 1, p ), 1, a( 1, q ), 1 )
937 IF( rsvec )CALL zswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
938 END IF
939 5991 CONTINUE
940*
941 RETURN
942* ..
943* .. END OF ZGSVJ0
944* ..
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zdotc
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
Definition zdotc.f:83
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
zlascl
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.
Definition zlascl.f:142
zlassq
subroutine zlassq(n, x, incx, scale, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.f90:122
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dznrm2
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
zrot
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.
Definition zrot.f:101
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: