LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dgsvj1()

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

DGSVJ1 pre-processor for the routine dgesvj, applies Jacobi rotations targeting only particular pivots.

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

Purpose:
 DGSVJ1 is called from DGESVJ as a pre-processor and that is its main
 purpose. It applies Jacobi rotations in the same way as DGESVJ 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
 ~~~~~~~~~~~~~~~
 DGSVJ1 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 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 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 DOUBLE PRECISION 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 DABS(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 DOUBLE PRECISION 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 dgsvj1.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  DOUBLE PRECISION EPS, SFMIN, TOL
243  INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
244  CHARACTER*1 JOBV
245 * ..
246 * .. Array Arguments ..
247  DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
248  $ WORK( LWORK )
249 * ..
250 *
251 * =====================================================================
252 *
253 * .. Local Parameters ..
254  DOUBLE PRECISION ZERO, HALF, ONE
255  parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
256 * ..
257 * .. Local Scalars ..
258  DOUBLE PRECISION 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  DOUBLE PRECISION FASTR( 5 )
269 * ..
270 * .. Intrinsic Functions ..
271  INTRINSIC dabs, max, dble, min, dsign, dsqrt
272 * ..
273 * .. External Functions ..
274  DOUBLE PRECISION DDOT, DNRM2
275  INTEGER IDAMAX
276  LOGICAL LSAME
277  EXTERNAL idamax, lsame, ddot, dnrm2
278 * ..
279 * .. External Subroutines ..
280  EXTERNAL daxpy, dcopy, dlascl, dlassq, drotm, dswap,
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( 'DGSVJ1', -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 = dsqrt( eps )
328  rootsfmin = dsqrt( sfmin )
329  small = sfmin / eps
330  big = one / sfmin
331  rootbig = one / rootsfmin
332  large = big / dsqrt( dble( m*n ) )
333  bigtheta = one / rooteps
334  roottol = dsqrt( 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 = ( ddot( m, a( 1, p ), 1, a( 1,
427  $ q ), 1 )*d( p )*d( q ) / aaqq )
428  $ / aapp
429  ELSE
430  CALL dcopy( m, a( 1, p ), 1, work, 1 )
431  CALL dlascl( 'G', 0, 0, aapp, d( p ),
432  $ m, 1, work, lda, ierr )
433  aapq = ddot( 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 = ( ddot( m, a( 1, p ), 1, a( 1,
444  $ q ), 1 )*d( p )*d( q ) / aaqq )
445  $ / aapp
446  ELSE
447  CALL dcopy( m, a( 1, q ), 1, work, 1 )
448  CALL dlascl( 'G', 0, 0, aaqq, d( q ),
449  $ m, 1, work, lda, ierr )
450  aapq = ddot( m, work, 1, a( 1, p ),
451  $ 1 )*d( p ) / aapp
452  END IF
453  END IF
454 
455  mxaapq = max( mxaapq, dabs( aapq ) )
456 
457 * TO rotate or NOT to rotate, THAT is the question ...
458 *
459  IF( dabs( 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*dabs(aqoap-apoaq) / aapq
470  IF( aaqq.GT.aapp0 )theta = -theta
471 
472  IF( dabs( 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 drotm( m, a( 1, p ), 1,
477  $ a( 1, q ), 1, fastr )
478  IF( rsvec )CALL drotm( mvl,
479  $ v( 1, p ), 1,
480  $ v( 1, q ), 1,
481  $ fastr )
482  sva( q ) = aaqq*dsqrt( max( zero,
483  $ one+t*apoaq*aapq ) )
484  aapp = aapp*dsqrt( max( zero,
485  $ one-t*aqoap*aapq ) )
486  mxsinj = max( mxsinj, dabs( t ) )
487  ELSE
488 *
489 * .. choose correct signum for THETA and rotate
490 *
491  thsign = -dsign( one, aapq )
492  IF( aaqq.GT.aapp0 )thsign = -thsign
493  t = one / ( theta+thsign*
494  $ dsqrt( one+theta*theta ) )
495  cs = dsqrt( one / ( one+t*t ) )
496  sn = t*cs
497  mxsinj = max( mxsinj, dabs( sn ) )
498  sva( q ) = aaqq*dsqrt( max( zero,
499  $ one+t*apoaq*aapq ) )
500  aapp = aapp*dsqrt( 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 drotm( m, a( 1, p ), 1,
513  $ a( 1, q ), 1,
514  $ fastr )
515  IF( rsvec )CALL drotm( mvl,
516  $ v( 1, p ), 1, v( 1, q ),
517  $ 1, fastr )
518  ELSE
519  CALL daxpy( m, -t*aqoap,
520  $ a( 1, q ), 1,
521  $ a( 1, p ), 1 )
522  CALL daxpy( m, cs*sn*apoaq,
523  $ a( 1, p ), 1,
524  $ a( 1, q ), 1 )
525  IF( rsvec ) THEN
526  CALL daxpy( mvl, -t*aqoap,
527  $ v( 1, q ), 1,
528  $ v( 1, p ), 1 )
529  CALL daxpy( 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 daxpy( m, t*apoaq,
540  $ a( 1, p ), 1,
541  $ a( 1, q ), 1 )
542  CALL daxpy( m, -cs*sn*aqoap,
543  $ a( 1, q ), 1,
544  $ a( 1, p ), 1 )
545  IF( rsvec ) THEN
546  CALL daxpy( mvl, t*apoaq,
547  $ v( 1, p ), 1,
548  $ v( 1, q ), 1 )
549  CALL daxpy( 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 daxpy( m, -t*aqoap,
559  $ a( 1, q ), 1,
560  $ a( 1, p ), 1 )
561  CALL daxpy( 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 daxpy( mvl,
568  $ -t*aqoap,
569  $ v( 1, q ), 1,
570  $ v( 1, p ), 1 )
571  CALL daxpy( mvl,
572  $ cs*sn*apoaq,
573  $ v( 1, p ), 1,
574  $ v( 1, q ), 1 )
575  END IF
576  ELSE
577  CALL daxpy( m, t*apoaq,
578  $ a( 1, p ), 1,
579  $ a( 1, q ), 1 )
580  CALL daxpy( 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 daxpy( mvl,
588  $ t*apoaq, v( 1, p ),
589  $ 1, v( 1, q ), 1 )
590  CALL daxpy( 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 dcopy( m, a( 1, p ), 1, work,
603  $ 1 )
604  CALL dlascl( 'G', 0, 0, aapp, one,
605  $ m, 1, work, lda, ierr )
606  CALL dlascl( 'G', 0, 0, aaqq, one,
607  $ m, 1, a( 1, q ), lda,
608  $ ierr )
609  temp1 = -aapq*d( p ) / d( q )
610  CALL daxpy( m, temp1, work, 1,
611  $ a( 1, q ), 1 )
612  CALL dlascl( 'G', 0, 0, one, aaqq,
613  $ m, 1, a( 1, q ), lda,
614  $ ierr )
615  sva( q ) = aaqq*dsqrt( max( zero,
616  $ one-aapq*aapq ) )
617  mxsinj = max( mxsinj, sfmin )
618  ELSE
619  CALL dcopy( m, a( 1, q ), 1, work,
620  $ 1 )
621  CALL dlascl( 'G', 0, 0, aaqq, one,
622  $ m, 1, work, lda, ierr )
623  CALL dlascl( 'G', 0, 0, aapp, one,
624  $ m, 1, a( 1, p ), lda,
625  $ ierr )
626  temp1 = -aapq*d( q ) / d( p )
627  CALL daxpy( m, temp1, work, 1,
628  $ a( 1, p ), 1 )
629  CALL dlascl( 'G', 0, 0, one, aapp,
630  $ m, 1, a( 1, p ), lda,
631  $ ierr )
632  sva( p ) = aapp*dsqrt( 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 ) = dnrm2( m, a( 1, q ), 1 )*
646  $ d( q )
647  ELSE
648  t = zero
649  aaqq = one
650  CALL dlassq( m, a( 1, q ), 1, t,
651  $ aaqq )
652  sva( q ) = t*dsqrt( 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 = dnrm2( m, a( 1, p ), 1 )*
659  $ d( p )
660  ELSE
661  t = zero
662  aapp = one
663  CALL dlassq( m, a( 1, p ), 1, t,
664  $ aapp )
665  aapp = t*dsqrt( 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 ) = dabs( 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 ) = dnrm2( m, a( 1, n ), 1 )*d( n )
727  ELSE
728  t = zero
729  aapp = one
730  CALL dlassq( m, a( 1, n ), 1, t, aapp )
731  sva( n ) = t*dsqrt( 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.dble( n )*tol ) .AND.
740  $ ( dble( 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 = idamax( 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 dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
773  IF( rsvec )CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
774  END IF
775  5991 CONTINUE
776 *
777  RETURN
778 * ..
779 * .. END OF DGSVJ1
780 * ..
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f90:137
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
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 drotm(N, DX, INCX, DY, INCY, DPARAM)
DROTM
Definition: drotm.f:96
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: