LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cgsvj0 ( character*1  JOBV,
integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( n )  D,
real, dimension( n )  SVA,
integer  MV,
complex, dimension( ldv, * )  V,
integer  LDV,
real  EPS,
real  SFMIN,
real  TOL,
integer  NSWEEP,
complex, dimension( lwork )  WORK,
integer  LWORK,
integer  INFO 
)

CGSVJ0 pre-processor for the routine cgesvj.

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

Purpose: 
 CGSVJ0 is called from CGESVJ as a pre-processor and that is its main
 purpose. It applies Jacobi rotations in the same way as CGESVJ 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 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 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 REAL 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 .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 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 REAL
          EPS = SLAMCH('Epsilon')
[in]SFMIN
          SFMIN is REAL
          SFMIN = SLAMCH('Safe Minimum')
[in]TOL
          TOL is REAL
          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 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:
CGSVJ0 is used just to enable CGESVJ 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 220 of file cgsvj0.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: