LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 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).
[in]MV
          MV is INTEGER
          If JOBV = '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 = 'V' then N rows of V are post-multipled by a
                           sequence of Jacobi rotations.
          If JOBV = '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 >= 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 216 of file zgsvj0.f.

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