LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ sgsvj1()

subroutine sgsvj1 ( character*1  JOBV,
integer  M,
integer  N,
integer  N1,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( n )  D,
real, dimension( n )  SVA,
integer  MV,
real, dimension( ldv, * )  V,
integer  LDV,
real  EPS,
real  SFMIN,
real  TOL,
integer  NSWEEP,
real, dimension( lwork )  WORK,
integer  LWORK,
integer  INFO 
)

SGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivots.

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

Purpose:
 SGSVJ1 is called from SGESVJ as a pre-processor and that is its main
 purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
 it targets only particular pivots and it does not check convergence
 (stopping criterion). Few tunning parameters (marked by [TP]) are
 available for the implementer.

 Further Details
 ~~~~~~~~~~~~~~~
 SGSVJ1 applies few sweeps of Jacobi rotations in the column space of
 the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
 off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
 block-entries (tiles) of the (1,2) off-diagonal block are marked by the
 [x]'s in the following scheme:

    | *  *  * [x] [x] [x]|
    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
    |[x] [x] [x] *  *  * |
    |[x] [x] [x] *  *  * |
    |[x] [x] [x] *  *  * |

 In terms of the columns of A, the first N1 columns are rotated 'against'
 the remaining N-N1 columns, trying to increase the angle between the
 corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
 tiled using quadratic tiles of side KBL. Here, KBL is a tunning parameter.
 The number of sweeps is given in NSWEEP and the orthogonality threshold
 is given in TOL.
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]N1
          N1 is INTEGER
          N1 specifies the 2 x 2 block partition, the first N1 columns are
          rotated 'against' the remaining N-N1 columns of A.
[in,out]A
          A is REAL array, dimension (LDA,N)
          On entry, M-by-N matrix A, such that A*diag(D) represents
          the input matrix.
          On exit,
          A_onexit * 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 N1, 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 REAL array, dimension (N)
          The array D accumulates the scaling factors from the fast 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 N1, 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 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 REAL 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 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)))) > TOL.
[in]NSWEEP
          NSWEEP is INTEGER
          NSWEEP is the number of sweeps of Jacobi rotations to be
          performed.
[out]WORK
          WORK is REAL 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.
Date
November 2017
Contributors:
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

Definition at line 238 of file sgsvj1.f.

238 *
239 * -- LAPACK computational routine (version 3.8.0) --
240 * -- LAPACK is a software package provided by Univ. of Tennessee, --
241 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242 * November 2017
243 *
244 * .. Scalar Arguments ..
245  REAL EPS, SFMIN, TOL
246  INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
247  CHARACTER*1 JOBV
248 * ..
249 * .. Array Arguments ..
250  REAL A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
251  $ WORK( LWORK )
252 * ..
253 *
254 * =====================================================================
255 *
256 * .. Local Parameters ..
257  REAL ZERO, HALF, ONE
258  parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
259 * ..
260 * .. Local Scalars ..
261  REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
262  $ BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
263  $ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
264  $ TEMP1, THETA, THSIGN
265  INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
266  $ ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
267  $ p, PSKIPPED, q, ROWSKIP, SWBAND
268  LOGICAL APPLV, ROTOK, RSVEC
269 * ..
270 * .. Local Arrays ..
271  REAL FASTR( 5 )
272 * ..
273 * .. Intrinsic Functions ..
274  INTRINSIC abs, max, float, min, sign, sqrt
275 * ..
276 * .. External Functions ..
277  REAL SDOT, SNRM2
278  INTEGER ISAMAX
279  LOGICAL LSAME
280  EXTERNAL isamax, lsame, sdot, snrm2
281 * ..
282 * .. External Subroutines ..
283  EXTERNAL saxpy, scopy, slascl, slassq, srotm, sswap,
284  $ xerbla
285 * ..
286 * .. Executable Statements ..
287 *
288 * Test the input parameters.
289 *
290  applv = lsame( jobv, 'A' )
291  rsvec = lsame( jobv, 'V' )
292  IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
293  info = -1
294  ELSE IF( m.LT.0 ) THEN
295  info = -2
296  ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
297  info = -3
298  ELSE IF( n1.LT.0 ) THEN
299  info = -4
300  ELSE IF( lda.LT.m ) THEN
301  info = -6
302  ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
303  info = -9
304  ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
305  $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
306  info = -11
307  ELSE IF( tol.LE.eps ) THEN
308  info = -14
309  ELSE IF( nsweep.LT.0 ) THEN
310  info = -15
311  ELSE IF( lwork.LT.m ) THEN
312  info = -17
313  ELSE
314  info = 0
315  END IF
316 *
317 * #:(
318  IF( info.NE.0 ) THEN
319  CALL xerbla( 'SGSVJ1', -info )
320  RETURN
321  END IF
322 *
323  IF( rsvec ) THEN
324  mvl = n
325  ELSE IF( applv ) THEN
326  mvl = mv
327  END IF
328  rsvec = rsvec .OR. applv
329 
330  rooteps = sqrt( eps )
331  rootsfmin = sqrt( sfmin )
332  small = sfmin / eps
333  big = one / sfmin
334  rootbig = one / rootsfmin
335  large = big / sqrt( float( m*n ) )
336  bigtheta = one / rooteps
337  roottol = sqrt( tol )
338 *
339 * .. Initialize the right singular vector matrix ..
340 *
341 * RSVEC = LSAME( JOBV, 'Y' )
342 *
343  emptsw = n1*( n-n1 )
344  notrot = 0
345  fastr( 1 ) = zero
346 *
347 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
348 *
349  kbl = min( 8, n )
350  nblr = n1 / kbl
351  IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
352 
353 * .. the tiling is nblr-by-nblc [tiles]
354 
355  nblc = ( n-n1 ) / kbl
356  IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
357  blskip = ( kbl**2 ) + 1
358 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
359 
360  rowskip = min( 5, kbl )
361 *[TP] ROWSKIP is a tuning parameter.
362  swband = 0
363 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
364 * if SGESVJ is used as a computational routine in the preconditioned
365 * Jacobi SVD algorithm SGESVJ.
366 *
367 *
368 * | * * * [x] [x] [x]|
369 * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
370 * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
371 * |[x] [x] [x] * * * |
372 * |[x] [x] [x] * * * |
373 * |[x] [x] [x] * * * |
374 *
375 *
376  DO 1993 i = 1, nsweep
377 * .. go go go ...
378 *
379  mxaapq = zero
380  mxsinj = zero
381  iswrot = 0
382 *
383  notrot = 0
384  pskipped = 0
385 *
386  DO 2000 ibr = 1, nblr
387 
388  igl = ( ibr-1 )*kbl + 1
389 *
390 *
391 *........................................................
392 * ... go to the off diagonal blocks
393 
394  igl = ( ibr-1 )*kbl + 1
395 
396  DO 2010 jbc = 1, nblc
397 
398  jgl = n1 + ( jbc-1 )*kbl + 1
399 
400 * doing the block at ( ibr, jbc )
401 
402  ijblsk = 0
403  DO 2100 p = igl, min( igl+kbl-1, n1 )
404 
405  aapp = sva( p )
406 
407  IF( aapp.GT.zero ) THEN
408 
409  pskipped = 0
410 
411  DO 2200 q = jgl, min( jgl+kbl-1, n )
412 *
413  aaqq = sva( q )
414 
415  IF( aaqq.GT.zero ) THEN
416  aapp0 = aapp
417 *
418 * .. M x 2 Jacobi SVD ..
419 *
420 * .. Safe Gram matrix computation ..
421 *
422  IF( aaqq.GE.one ) THEN
423  IF( aapp.GE.aaqq ) THEN
424  rotok = ( small*aapp ).LE.aaqq
425  ELSE
426  rotok = ( small*aaqq ).LE.aapp
427  END IF
428  IF( aapp.LT.( big / aaqq ) ) THEN
429  aapq = ( sdot( m, a( 1, p ), 1, a( 1,
430  $ q ), 1 )*d( p )*d( q ) / aaqq )
431  $ / aapp
432  ELSE
433  CALL scopy( m, a( 1, p ), 1, work, 1 )
434  CALL slascl( 'G', 0, 0, aapp, d( p ),
435  $ m, 1, work, lda, ierr )
436  aapq = sdot( m, work, 1, a( 1, q ),
437  $ 1 )*d( q ) / aaqq
438  END IF
439  ELSE
440  IF( aapp.GE.aaqq ) THEN
441  rotok = aapp.LE.( aaqq / small )
442  ELSE
443  rotok = aaqq.LE.( aapp / small )
444  END IF
445  IF( aapp.GT.( small / aaqq ) ) THEN
446  aapq = ( sdot( m, a( 1, p ), 1, a( 1,
447  $ q ), 1 )*d( p )*d( q ) / aaqq )
448  $ / aapp
449  ELSE
450  CALL scopy( m, a( 1, q ), 1, work, 1 )
451  CALL slascl( 'G', 0, 0, aaqq, d( q ),
452  $ m, 1, work, lda, ierr )
453  aapq = sdot( m, work, 1, a( 1, p ),
454  $ 1 )*d( p ) / aapp
455  END IF
456  END IF
457 
458  mxaapq = max( mxaapq, abs( aapq ) )
459 
460 * TO rotate or NOT to rotate, THAT is the question ...
461 *
462  IF( abs( aapq ).GT.tol ) THEN
463  notrot = 0
464 * ROTATED = ROTATED + 1
465  pskipped = 0
466  iswrot = iswrot + 1
467 *
468  IF( rotok ) THEN
469 *
470  aqoap = aaqq / aapp
471  apoaq = aapp / aaqq
472  theta = -half*abs( aqoap-apoaq ) / aapq
473  IF( aaqq.GT.aapp0 )theta = -theta
474 
475  IF( abs( theta ).GT.bigtheta ) THEN
476  t = half / theta
477  fastr( 3 ) = t*d( p ) / d( q )
478  fastr( 4 ) = -t*d( q ) / d( p )
479  CALL srotm( m, a( 1, p ), 1,
480  $ a( 1, q ), 1, fastr )
481  IF( rsvec )CALL srotm( mvl,
482  $ v( 1, p ), 1,
483  $ v( 1, q ), 1,
484  $ fastr )
485  sva( q ) = aaqq*sqrt( max( zero,
486  $ one+t*apoaq*aapq ) )
487  aapp = aapp*sqrt( max( zero,
488  $ one-t*aqoap*aapq ) )
489  mxsinj = max( mxsinj, abs( t ) )
490  ELSE
491 *
492 * .. choose correct signum for THETA and rotate
493 *
494  thsign = -sign( one, aapq )
495  IF( aaqq.GT.aapp0 )thsign = -thsign
496  t = one / ( theta+thsign*
497  $ sqrt( one+theta*theta ) )
498  cs = sqrt( one / ( one+t*t ) )
499  sn = t*cs
500  mxsinj = max( mxsinj, abs( sn ) )
501  sva( q ) = aaqq*sqrt( max( zero,
502  $ one+t*apoaq*aapq ) )
503  aapp = aapp*sqrt( max( zero,
504  $ one-t*aqoap*aapq ) )
505 
506  apoaq = d( p ) / d( q )
507  aqoap = d( q ) / d( p )
508  IF( d( p ).GE.one ) THEN
509 *
510  IF( d( q ).GE.one ) THEN
511  fastr( 3 ) = t*apoaq
512  fastr( 4 ) = -t*aqoap
513  d( p ) = d( p )*cs
514  d( q ) = d( q )*cs
515  CALL srotm( m, a( 1, p ), 1,
516  $ a( 1, q ), 1,
517  $ fastr )
518  IF( rsvec )CALL srotm( mvl,
519  $ v( 1, p ), 1, v( 1, q ),
520  $ 1, fastr )
521  ELSE
522  CALL saxpy( m, -t*aqoap,
523  $ a( 1, q ), 1,
524  $ a( 1, p ), 1 )
525  CALL saxpy( m, cs*sn*apoaq,
526  $ a( 1, p ), 1,
527  $ a( 1, q ), 1 )
528  IF( rsvec ) THEN
529  CALL saxpy( mvl, -t*aqoap,
530  $ v( 1, q ), 1,
531  $ v( 1, p ), 1 )
532  CALL saxpy( mvl,
533  $ cs*sn*apoaq,
534  $ v( 1, p ), 1,
535  $ v( 1, q ), 1 )
536  END IF
537  d( p ) = d( p )*cs
538  d( q ) = d( q ) / cs
539  END IF
540  ELSE
541  IF( d( q ).GE.one ) THEN
542  CALL saxpy( m, t*apoaq,
543  $ a( 1, p ), 1,
544  $ a( 1, q ), 1 )
545  CALL saxpy( m, -cs*sn*aqoap,
546  $ a( 1, q ), 1,
547  $ a( 1, p ), 1 )
548  IF( rsvec ) THEN
549  CALL saxpy( mvl, t*apoaq,
550  $ v( 1, p ), 1,
551  $ v( 1, q ), 1 )
552  CALL saxpy( mvl,
553  $ -cs*sn*aqoap,
554  $ v( 1, q ), 1,
555  $ v( 1, p ), 1 )
556  END IF
557  d( p ) = d( p ) / cs
558  d( q ) = d( q )*cs
559  ELSE
560  IF( d( p ).GE.d( q ) ) THEN
561  CALL saxpy( m, -t*aqoap,
562  $ a( 1, q ), 1,
563  $ a( 1, p ), 1 )
564  CALL saxpy( m, cs*sn*apoaq,
565  $ a( 1, p ), 1,
566  $ a( 1, q ), 1 )
567  d( p ) = d( p )*cs
568  d( q ) = d( q ) / cs
569  IF( rsvec ) THEN
570  CALL saxpy( mvl,
571  $ -t*aqoap,
572  $ v( 1, q ), 1,
573  $ v( 1, p ), 1 )
574  CALL saxpy( mvl,
575  $ cs*sn*apoaq,
576  $ v( 1, p ), 1,
577  $ v( 1, q ), 1 )
578  END IF
579  ELSE
580  CALL saxpy( m, t*apoaq,
581  $ a( 1, p ), 1,
582  $ a( 1, q ), 1 )
583  CALL saxpy( m,
584  $ -cs*sn*aqoap,
585  $ a( 1, q ), 1,
586  $ a( 1, p ), 1 )
587  d( p ) = d( p ) / cs
588  d( q ) = d( q )*cs
589  IF( rsvec ) THEN
590  CALL saxpy( mvl,
591  $ t*apoaq, v( 1, p ),
592  $ 1, v( 1, q ), 1 )
593  CALL saxpy( mvl,
594  $ -cs*sn*aqoap,
595  $ v( 1, q ), 1,
596  $ v( 1, p ), 1 )
597  END IF
598  END IF
599  END IF
600  END IF
601  END IF
602 
603  ELSE
604  IF( aapp.GT.aaqq ) THEN
605  CALL scopy( m, a( 1, p ), 1, work,
606  $ 1 )
607  CALL slascl( 'G', 0, 0, aapp, one,
608  $ m, 1, work, lda, ierr )
609  CALL slascl( 'G', 0, 0, aaqq, one,
610  $ m, 1, a( 1, q ), lda,
611  $ ierr )
612  temp1 = -aapq*d( p ) / d( q )
613  CALL saxpy( m, temp1, work, 1,
614  $ a( 1, q ), 1 )
615  CALL slascl( 'G', 0, 0, one, aaqq,
616  $ m, 1, a( 1, q ), lda,
617  $ ierr )
618  sva( q ) = aaqq*sqrt( max( zero,
619  $ one-aapq*aapq ) )
620  mxsinj = max( mxsinj, sfmin )
621  ELSE
622  CALL scopy( m, a( 1, q ), 1, work,
623  $ 1 )
624  CALL slascl( 'G', 0, 0, aaqq, one,
625  $ m, 1, work, lda, ierr )
626  CALL slascl( 'G', 0, 0, aapp, one,
627  $ m, 1, a( 1, p ), lda,
628  $ ierr )
629  temp1 = -aapq*d( q ) / d( p )
630  CALL saxpy( m, temp1, work, 1,
631  $ a( 1, p ), 1 )
632  CALL slascl( 'G', 0, 0, one, aapp,
633  $ m, 1, a( 1, p ), lda,
634  $ ierr )
635  sva( p ) = aapp*sqrt( max( zero,
636  $ one-aapq*aapq ) )
637  mxsinj = max( mxsinj, sfmin )
638  END IF
639  END IF
640 * END IF ROTOK THEN ... ELSE
641 *
642 * In the case of cancellation in updating SVA(q)
643 * .. recompute SVA(q)
644  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
645  $ THEN
646  IF( ( aaqq.LT.rootbig ) .AND.
647  $ ( aaqq.GT.rootsfmin ) ) THEN
648  sva( q ) = snrm2( m, a( 1, q ), 1 )*
649  $ d( q )
650  ELSE
651  t = zero
652  aaqq = one
653  CALL slassq( m, a( 1, q ), 1, t,
654  $ aaqq )
655  sva( q ) = t*sqrt( aaqq )*d( q )
656  END IF
657  END IF
658  IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
659  IF( ( aapp.LT.rootbig ) .AND.
660  $ ( aapp.GT.rootsfmin ) ) THEN
661  aapp = snrm2( m, a( 1, p ), 1 )*
662  $ d( p )
663  ELSE
664  t = zero
665  aapp = one
666  CALL slassq( m, a( 1, p ), 1, t,
667  $ aapp )
668  aapp = t*sqrt( aapp )*d( p )
669  END IF
670  sva( p ) = aapp
671  END IF
672 * end of OK rotation
673  ELSE
674  notrot = notrot + 1
675 * SKIPPED = SKIPPED + 1
676  pskipped = pskipped + 1
677  ijblsk = ijblsk + 1
678  END IF
679  ELSE
680  notrot = notrot + 1
681  pskipped = pskipped + 1
682  ijblsk = ijblsk + 1
683  END IF
684 
685 * IF ( NOTROT .GE. EMPTSW ) GO TO 2011
686  IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
687  $ THEN
688  sva( p ) = aapp
689  notrot = 0
690  GO TO 2011
691  END IF
692  IF( ( i.LE.swband ) .AND.
693  $ ( pskipped.GT.rowskip ) ) THEN
694  aapp = -aapp
695  notrot = 0
696  GO TO 2203
697  END IF
698 
699 *
700  2200 CONTINUE
701 * end of the q-loop
702  2203 CONTINUE
703 
704  sva( p ) = aapp
705 *
706  ELSE
707  IF( aapp.EQ.zero )notrot = notrot +
708  $ min( jgl+kbl-1, n ) - jgl + 1
709  IF( aapp.LT.zero )notrot = 0
710 *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
711  END IF
712 
713  2100 CONTINUE
714 * end of the p-loop
715  2010 CONTINUE
716 * end of the jbc-loop
717  2011 CONTINUE
718 *2011 bailed out of the jbc-loop
719  DO 2012 p = igl, min( igl+kbl-1, n )
720  sva( p ) = abs( sva( p ) )
721  2012 CONTINUE
722 *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
723  2000 CONTINUE
724 *2000 :: end of the ibr-loop
725 *
726 * .. update SVA(N)
727  IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
728  $ THEN
729  sva( n ) = snrm2( m, a( 1, n ), 1 )*d( n )
730  ELSE
731  t = zero
732  aapp = one
733  CALL slassq( m, a( 1, n ), 1, t, aapp )
734  sva( n ) = t*sqrt( aapp )*d( n )
735  END IF
736 *
737 * Additional steering devices
738 *
739  IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
740  $ ( iswrot.LE.n ) ) )swband = i
741 
742  IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.float( n )*tol ) .AND.
743  $ ( float( n )*mxaapq*mxsinj.LT.tol ) ) THEN
744  GO TO 1994
745  END IF
746 
747 *
748  IF( notrot.GE.emptsw )GO TO 1994
749 
750  1993 CONTINUE
751 * end i=1:NSWEEP loop
752 * #:) Reaching this point means that the procedure has completed the given
753 * number of sweeps.
754  info = nsweep - 1
755  GO TO 1995
756  1994 CONTINUE
757 * #:) Reaching this point means that during the i-th sweep all pivots were
758 * below the given threshold, causing early exit.
759 
760  info = 0
761 * #:) INFO = 0 confirms successful iterations.
762  1995 CONTINUE
763 *
764 * Sort the vector D
765 *
766  DO 5991 p = 1, n - 1
767  q = isamax( n-p+1, sva( p ), 1 ) + p - 1
768  IF( p.NE.q ) THEN
769  temp1 = sva( p )
770  sva( p ) = sva( q )
771  sva( q ) = temp1
772  temp1 = d( p )
773  d( p ) = d( q )
774  d( q ) = temp1
775  CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
776  IF( rsvec )CALL sswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
777  END IF
778  5991 CONTINUE
779 *
780  RETURN
781 * ..
782 * .. END OF SGSVJ1
783 * ..
Here is the call graph for this function:
Here is the caller graph for this function:
snrm2
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:76
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
srotm
subroutine srotm(N, SX, INCX, SY, INCY, SPARAM)
SROTM
Definition: srotm.f:99
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:73
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
slassq
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f:105
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
sdot
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:84
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:91