LAPACK  3.8.0
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 tunning 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 tunning parmeter.
 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 .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 DOUBLE PRECISION 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 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)))) .GT. 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 .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
Contributors:
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

Definition at line 238 of file dgsvj1.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 * June 2016
243 *
244 * .. Scalar Arguments ..
245  DOUBLE PRECISION eps, sfmin, tol
246  INTEGER info, lda, ldv, lwork, m, mv, n, n1, nsweep
247  CHARACTER*1 jobv
248 * ..
249 * .. Array Arguments ..
250  DOUBLE PRECISION a( lda, * ), d( n ), sva( n ), v( ldv, * ),
251  $ work( lwork )
252 * ..
253 *
254 * =====================================================================
255 *
256 * .. Local Parameters ..
257  DOUBLE PRECISION zero, half, one
258  parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
259 * ..
260 * .. Local Scalars ..
261  DOUBLE PRECISION 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  DOUBLE PRECISION fastr( 5 )
272 * ..
273 * .. Intrinsic Functions ..
274  INTRINSIC dabs, max, dble, min, dsign, dsqrt
275 * ..
276 * .. External Functions ..
277  DOUBLE PRECISION ddot, dnrm2
278  INTEGER idamax
279  LOGICAL lsame
280  EXTERNAL idamax, lsame, ddot, dnrm2
281 * ..
282 * .. External Subroutines ..
283  EXTERNAL daxpy, dcopy, dlascl, dlassq, drotm, dswap,
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( 'DGSVJ1', -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 = dsqrt( eps )
331  rootsfmin = dsqrt( sfmin )
332  small = sfmin / eps
333  big = one / sfmin
334  rootbig = one / rootsfmin
335  large = big / dsqrt( dble( m*n ) )
336  bigtheta = one / rooteps
337  roottol = dsqrt( 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 = ( ddot( m, a( 1, p ), 1, a( 1,
430  $ q ), 1 )*d( p )*d( q ) / aaqq )
431  $ / aapp
432  ELSE
433  CALL dcopy( m, a( 1, p ), 1, work, 1 )
434  CALL dlascl( 'G', 0, 0, aapp, d( p ),
435  $ m, 1, work, lda, ierr )
436  aapq = ddot( 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 = ( ddot( m, a( 1, p ), 1, a( 1,
447  $ q ), 1 )*d( p )*d( q ) / aaqq )
448  $ / aapp
449  ELSE
450  CALL dcopy( m, a( 1, q ), 1, work, 1 )
451  CALL dlascl( 'G', 0, 0, aaqq, d( q ),
452  $ m, 1, work, lda, ierr )
453  aapq = ddot( m, work, 1, a( 1, p ),
454  $ 1 )*d( p ) / aapp
455  END IF
456  END IF
457 
458  mxaapq = max( mxaapq, dabs( aapq ) )
459 
460 * TO rotate or NOT to rotate, THAT is the question ...
461 *
462  IF( dabs( 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*dabs(aqoap-apoaq) / aapq
473  IF( aaqq.GT.aapp0 )theta = -theta
474 
475  IF( dabs( 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 drotm( m, a( 1, p ), 1,
480  $ a( 1, q ), 1, fastr )
481  IF( rsvec )CALL drotm( mvl,
482  $ v( 1, p ), 1,
483  $ v( 1, q ), 1,
484  $ fastr )
485  sva( q ) = aaqq*dsqrt( max( zero,
486  $ one+t*apoaq*aapq ) )
487  aapp = aapp*dsqrt( max( zero,
488  $ one-t*aqoap*aapq ) )
489  mxsinj = max( mxsinj, dabs( t ) )
490  ELSE
491 *
492 * .. choose correct signum for THETA and rotate
493 *
494  thsign = -dsign( one, aapq )
495  IF( aaqq.GT.aapp0 )thsign = -thsign
496  t = one / ( theta+thsign*
497  $ dsqrt( one+theta*theta ) )
498  cs = dsqrt( one / ( one+t*t ) )
499  sn = t*cs
500  mxsinj = max( mxsinj, dabs( sn ) )
501  sva( q ) = aaqq*dsqrt( max( zero,
502  $ one+t*apoaq*aapq ) )
503  aapp = aapp*dsqrt( 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 drotm( m, a( 1, p ), 1,
516  $ a( 1, q ), 1,
517  $ fastr )
518  IF( rsvec )CALL drotm( mvl,
519  $ v( 1, p ), 1, v( 1, q ),
520  $ 1, fastr )
521  ELSE
522  CALL daxpy( m, -t*aqoap,
523  $ a( 1, q ), 1,
524  $ a( 1, p ), 1 )
525  CALL daxpy( m, cs*sn*apoaq,
526  $ a( 1, p ), 1,
527  $ a( 1, q ), 1 )
528  IF( rsvec ) THEN
529  CALL daxpy( mvl, -t*aqoap,
530  $ v( 1, q ), 1,
531  $ v( 1, p ), 1 )
532  CALL daxpy( 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 daxpy( m, t*apoaq,
543  $ a( 1, p ), 1,
544  $ a( 1, q ), 1 )
545  CALL daxpy( m, -cs*sn*aqoap,
546  $ a( 1, q ), 1,
547  $ a( 1, p ), 1 )
548  IF( rsvec ) THEN
549  CALL daxpy( mvl, t*apoaq,
550  $ v( 1, p ), 1,
551  $ v( 1, q ), 1 )
552  CALL daxpy( 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 daxpy( m, -t*aqoap,
562  $ a( 1, q ), 1,
563  $ a( 1, p ), 1 )
564  CALL daxpy( 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 daxpy( mvl,
571  $ -t*aqoap,
572  $ v( 1, q ), 1,
573  $ v( 1, p ), 1 )
574  CALL daxpy( mvl,
575  $ cs*sn*apoaq,
576  $ v( 1, p ), 1,
577  $ v( 1, q ), 1 )
578  END IF
579  ELSE
580  CALL daxpy( m, t*apoaq,
581  $ a( 1, p ), 1,
582  $ a( 1, q ), 1 )
583  CALL daxpy( 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 daxpy( mvl,
591  $ t*apoaq, v( 1, p ),
592  $ 1, v( 1, q ), 1 )
593  CALL daxpy( 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 dcopy( m, a( 1, p ), 1, work,
606  $ 1 )
607  CALL dlascl( 'G', 0, 0, aapp, one,
608  $ m, 1, work, lda, ierr )
609  CALL dlascl( 'G', 0, 0, aaqq, one,
610  $ m, 1, a( 1, q ), lda,
611  $ ierr )
612  temp1 = -aapq*d( p ) / d( q )
613  CALL daxpy( m, temp1, work, 1,
614  $ a( 1, q ), 1 )
615  CALL dlascl( 'G', 0, 0, one, aaqq,
616  $ m, 1, a( 1, q ), lda,
617  $ ierr )
618  sva( q ) = aaqq*dsqrt( max( zero,
619  $ one-aapq*aapq ) )
620  mxsinj = max( mxsinj, sfmin )
621  ELSE
622  CALL dcopy( m, a( 1, q ), 1, work,
623  $ 1 )
624  CALL dlascl( 'G', 0, 0, aaqq, one,
625  $ m, 1, work, lda, ierr )
626  CALL dlascl( 'G', 0, 0, aapp, one,
627  $ m, 1, a( 1, p ), lda,
628  $ ierr )
629  temp1 = -aapq*d( q ) / d( p )
630  CALL daxpy( m, temp1, work, 1,
631  $ a( 1, p ), 1 )
632  CALL dlascl( 'G', 0, 0, one, aapp,
633  $ m, 1, a( 1, p ), lda,
634  $ ierr )
635  sva( p ) = aapp*dsqrt( 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 ) = dnrm2( m, a( 1, q ), 1 )*
649  $ d( q )
650  ELSE
651  t = zero
652  aaqq = one
653  CALL dlassq( m, a( 1, q ), 1, t,
654  $ aaqq )
655  sva( q ) = t*dsqrt( 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 = dnrm2( m, a( 1, p ), 1 )*
662  $ d( p )
663  ELSE
664  t = zero
665  aapp = one
666  CALL dlassq( m, a( 1, p ), 1, t,
667  $ aapp )
668  aapp = t*dsqrt( 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 ) = dabs( 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 ) = dnrm2( m, a( 1, n ), 1 )*d( n )
730  ELSE
731  t = zero
732  aapp = one
733  CALL dlassq( m, a( 1, n ), 1, t, aapp )
734  sva( n ) = t*dsqrt( 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.dble( n )*tol ) .AND.
743  $ ( dble( 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 = idamax( 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 dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
776  IF( rsvec )CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
777  END IF
778  5991 CONTINUE
779 *
780  RETURN
781 * ..
782 * .. END OF DGSVJ1
783 * ..
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:73
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:84
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:145
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:84
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
subroutine drotm(N, DX, INCX, DY, INCY, DPARAM)
DROTM
Definition: drotm.f:98
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
Here is the call graph for this function:
Here is the caller graph for this function: