LAPACK  3.10.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 tuning 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 tuning 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.
Contributors:
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

Definition at line 234 of file sgsvj1.f.

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