LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 postmulyiplying the N-by-N array V.
                (See the description of V.)
          = 'A': the product of the Jacobi rotations is accumulated
                 by postmulyiplying 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).

 \param[in] MV
 \verbatim
          MV is INTEGER
          If JOBV .EQ. 'A', then MV rows of V are post-multipled 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 .EQ. 'V' then N rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV .EQ. 'A' then MV rows of V are post-multipled 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 .GE. N.
          If JOBV = 'A', LDV .GE. 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)))) .GT. 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 .GE. 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.
Date
June 2016
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.

Contributors: Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

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 219 of file zgsvj0.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: