LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.6.1) --
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 * ..
285 * .. Executable Statements ..
286 *
287 * Test the input parameters.
288 *
289  applv = lsame( jobv, 'A' )
290  rsvec = lsame( jobv, 'V' )
291  IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
292  info = -1
293  ELSE IF( m.LT.0 ) THEN
294  info = -2
295  ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
296  info = -3
297  ELSE IF( n1.LT.0 ) THEN
298  info = -4
299  ELSE IF( lda.LT.m ) THEN
300  info = -6
301  ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
302  info = -9
303  ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
304  $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
305  info = -11
306  ELSE IF( tol.LE.eps ) THEN
307  info = -14
308  ELSE IF( nsweep.LT.0 ) THEN
309  info = -15
310  ELSE IF( lwork.LT.m ) THEN
311  info = -17
312  ELSE
313  info = 0
314  END IF
315 *
316 * #:(
317  IF( info.NE.0 ) THEN
318  CALL xerbla( 'DGSVJ1', -info )
319  RETURN
320  END IF
321 *
322  IF( rsvec ) THEN
323  mvl = n
324  ELSE IF( applv ) THEN
325  mvl = mv
326  END IF
327  rsvec = rsvec .OR. applv
328 
329  rooteps = dsqrt( eps )
330  rootsfmin = dsqrt( sfmin )
331  small = sfmin / eps
332  big = one / sfmin
333  rootbig = one / rootsfmin
334  large = big / dsqrt( dble( m*n ) )
335  bigtheta = one / rooteps
336  roottol = dsqrt( tol )
337 *
338 * .. Initialize the right singular vector matrix ..
339 *
340 * RSVEC = LSAME( JOBV, 'Y' )
341 *
342  emptsw = n1*( n-n1 )
343  notrot = 0
344  fastr( 1 ) = zero
345 *
346 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
347 *
348  kbl = min( 8, n )
349  nblr = n1 / kbl
350  IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
351 
352 * .. the tiling is nblr-by-nblc [tiles]
353 
354  nblc = ( n-n1 ) / kbl
355  IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
356  blskip = ( kbl**2 ) + 1
357 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
358 
359  rowskip = min( 5, kbl )
360 *[TP] ROWSKIP is a tuning parameter.
361  swband = 0
362 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
363 * if SGESVJ is used as a computational routine in the preconditioned
364 * Jacobi SVD algorithm SGESVJ.
365 *
366 *
367 * | * * * [x] [x] [x]|
368 * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
369 * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
370 * |[x] [x] [x] * * * |
371 * |[x] [x] [x] * * * |
372 * |[x] [x] [x] * * * |
373 *
374 *
375  DO 1993 i = 1, nsweep
376 * .. go go go ...
377 *
378  mxaapq = zero
379  mxsinj = zero
380  iswrot = 0
381 *
382  notrot = 0
383  pskipped = 0
384 *
385  DO 2000 ibr = 1, nblr
386 
387  igl = ( ibr-1 )*kbl + 1
388 *
389 *
390 *........................................................
391 * ... go to the off diagonal blocks
392 
393  igl = ( ibr-1 )*kbl + 1
394 
395  DO 2010 jbc = 1, nblc
396 
397  jgl = n1 + ( jbc-1 )*kbl + 1
398 
399 * doing the block at ( ibr, jbc )
400 
401  ijblsk = 0
402  DO 2100 p = igl, min( igl+kbl-1, n1 )
403 
404  aapp = sva( p )
405 
406  IF( aapp.GT.zero ) THEN
407 
408  pskipped = 0
409 
410  DO 2200 q = jgl, min( jgl+kbl-1, n )
411 *
412  aaqq = sva( q )
413 
414  IF( aaqq.GT.zero ) THEN
415  aapp0 = aapp
416 *
417 * .. M x 2 Jacobi SVD ..
418 *
419 * .. Safe Gram matrix computation ..
420 *
421  IF( aaqq.GE.one ) THEN
422  IF( aapp.GE.aaqq ) THEN
423  rotok = ( small*aapp ).LE.aaqq
424  ELSE
425  rotok = ( small*aaqq ).LE.aapp
426  END IF
427  IF( aapp.LT.( big / aaqq ) ) THEN
428  aapq = ( ddot( m, a( 1, p ), 1, a( 1,
429  $ q ), 1 )*d( p )*d( q ) / aaqq )
430  $ / aapp
431  ELSE
432  CALL dcopy( m, a( 1, p ), 1, work, 1 )
433  CALL dlascl( 'G', 0, 0, aapp, d( p ),
434  $ m, 1, work, lda, ierr )
435  aapq = ddot( m, work, 1, a( 1, q ),
436  $ 1 )*d( q ) / aaqq
437  END IF
438  ELSE
439  IF( aapp.GE.aaqq ) THEN
440  rotok = aapp.LE.( aaqq / small )
441  ELSE
442  rotok = aaqq.LE.( aapp / small )
443  END IF
444  IF( aapp.GT.( small / aaqq ) ) THEN
445  aapq = ( ddot( m, a( 1, p ), 1, a( 1,
446  $ q ), 1 )*d( p )*d( q ) / aaqq )
447  $ / aapp
448  ELSE
449  CALL dcopy( m, a( 1, q ), 1, work, 1 )
450  CALL dlascl( 'G', 0, 0, aaqq, d( q ),
451  $ m, 1, work, lda, ierr )
452  aapq = ddot( m, work, 1, a( 1, p ),
453  $ 1 )*d( p ) / aapp
454  END IF
455  END IF
456 
457  mxaapq = max( mxaapq, dabs( aapq ) )
458 
459 * TO rotate or NOT to rotate, THAT is the question ...
460 *
461  IF( dabs( aapq ).GT.tol ) THEN
462  notrot = 0
463 * ROTATED = ROTATED + 1
464  pskipped = 0
465  iswrot = iswrot + 1
466 *
467  IF( rotok ) THEN
468 *
469  aqoap = aaqq / aapp
470  apoaq = aapp / aaqq
471  theta = -half*dabs(aqoap-apoaq) / aapq
472  IF( aaqq.GT.aapp0 )theta = -theta
473 
474  IF( dabs( theta ).GT.bigtheta ) THEN
475  t = half / theta
476  fastr( 3 ) = t*d( p ) / d( q )
477  fastr( 4 ) = -t*d( q ) / d( p )
478  CALL drotm( m, a( 1, p ), 1,
479  $ a( 1, q ), 1, fastr )
480  IF( rsvec )CALL drotm( mvl,
481  $ v( 1, p ), 1,
482  $ v( 1, q ), 1,
483  $ fastr )
484  sva( q ) = aaqq*dsqrt( max( zero,
485  $ one+t*apoaq*aapq ) )
486  aapp = aapp*dsqrt( max( zero,
487  $ one-t*aqoap*aapq ) )
488  mxsinj = max( mxsinj, dabs( t ) )
489  ELSE
490 *
491 * .. choose correct signum for THETA and rotate
492 *
493  thsign = -dsign( one, aapq )
494  IF( aaqq.GT.aapp0 )thsign = -thsign
495  t = one / ( theta+thsign*
496  $ dsqrt( one+theta*theta ) )
497  cs = dsqrt( one / ( one+t*t ) )
498  sn = t*cs
499  mxsinj = max( mxsinj, dabs( sn ) )
500  sva( q ) = aaqq*dsqrt( max( zero,
501  $ one+t*apoaq*aapq ) )
502  aapp = aapp*dsqrt( max( zero,
503  $ one-t*aqoap*aapq ) )
504 
505  apoaq = d( p ) / d( q )
506  aqoap = d( q ) / d( p )
507  IF( d( p ).GE.one ) THEN
508 *
509  IF( d( q ).GE.one ) THEN
510  fastr( 3 ) = t*apoaq
511  fastr( 4 ) = -t*aqoap
512  d( p ) = d( p )*cs
513  d( q ) = d( q )*cs
514  CALL drotm( m, a( 1, p ), 1,
515  $ a( 1, q ), 1,
516  $ fastr )
517  IF( rsvec )CALL drotm( mvl,
518  $ v( 1, p ), 1, v( 1, q ),
519  $ 1, fastr )
520  ELSE
521  CALL daxpy( m, -t*aqoap,
522  $ a( 1, q ), 1,
523  $ a( 1, p ), 1 )
524  CALL daxpy( m, cs*sn*apoaq,
525  $ a( 1, p ), 1,
526  $ a( 1, q ), 1 )
527  IF( rsvec ) THEN
528  CALL daxpy( mvl, -t*aqoap,
529  $ v( 1, q ), 1,
530  $ v( 1, p ), 1 )
531  CALL daxpy( mvl,
532  $ cs*sn*apoaq,
533  $ v( 1, p ), 1,
534  $ v( 1, q ), 1 )
535  END IF
536  d( p ) = d( p )*cs
537  d( q ) = d( q ) / cs
538  END IF
539  ELSE
540  IF( d( q ).GE.one ) THEN
541  CALL daxpy( m, t*apoaq,
542  $ a( 1, p ), 1,
543  $ a( 1, q ), 1 )
544  CALL daxpy( m, -cs*sn*aqoap,
545  $ a( 1, q ), 1,
546  $ a( 1, p ), 1 )
547  IF( rsvec ) THEN
548  CALL daxpy( mvl, t*apoaq,
549  $ v( 1, p ), 1,
550  $ v( 1, q ), 1 )
551  CALL daxpy( mvl,
552  $ -cs*sn*aqoap,
553  $ v( 1, q ), 1,
554  $ v( 1, p ), 1 )
555  END IF
556  d( p ) = d( p ) / cs
557  d( q ) = d( q )*cs
558  ELSE
559  IF( d( p ).GE.d( q ) ) THEN
560  CALL daxpy( m, -t*aqoap,
561  $ a( 1, q ), 1,
562  $ a( 1, p ), 1 )
563  CALL daxpy( m, cs*sn*apoaq,
564  $ a( 1, p ), 1,
565  $ a( 1, q ), 1 )
566  d( p ) = d( p )*cs
567  d( q ) = d( q ) / cs
568  IF( rsvec ) THEN
569  CALL daxpy( mvl,
570  $ -t*aqoap,
571  $ v( 1, q ), 1,
572  $ v( 1, p ), 1 )
573  CALL daxpy( mvl,
574  $ cs*sn*apoaq,
575  $ v( 1, p ), 1,
576  $ v( 1, q ), 1 )
577  END IF
578  ELSE
579  CALL daxpy( m, t*apoaq,
580  $ a( 1, p ), 1,
581  $ a( 1, q ), 1 )
582  CALL daxpy( m,
583  $ -cs*sn*aqoap,
584  $ a( 1, q ), 1,
585  $ a( 1, p ), 1 )
586  d( p ) = d( p ) / cs
587  d( q ) = d( q )*cs
588  IF( rsvec ) THEN
589  CALL daxpy( mvl,
590  $ t*apoaq, v( 1, p ),
591  $ 1, v( 1, q ), 1 )
592  CALL daxpy( mvl,
593  $ -cs*sn*aqoap,
594  $ v( 1, q ), 1,
595  $ v( 1, p ), 1 )
596  END IF
597  END IF
598  END IF
599  END IF
600  END IF
601 
602  ELSE
603  IF( aapp.GT.aaqq ) THEN
604  CALL dcopy( m, a( 1, p ), 1, work,
605  $ 1 )
606  CALL dlascl( 'G', 0, 0, aapp, one,
607  $ m, 1, work, lda, ierr )
608  CALL dlascl( 'G', 0, 0, aaqq, one,
609  $ m, 1, a( 1, q ), lda,
610  $ ierr )
611  temp1 = -aapq*d( p ) / d( q )
612  CALL daxpy( m, temp1, work, 1,
613  $ a( 1, q ), 1 )
614  CALL dlascl( 'G', 0, 0, one, aaqq,
615  $ m, 1, a( 1, q ), lda,
616  $ ierr )
617  sva( q ) = aaqq*dsqrt( max( zero,
618  $ one-aapq*aapq ) )
619  mxsinj = max( mxsinj, sfmin )
620  ELSE
621  CALL dcopy( m, a( 1, q ), 1, work,
622  $ 1 )
623  CALL dlascl( 'G', 0, 0, aaqq, one,
624  $ m, 1, work, lda, ierr )
625  CALL dlascl( 'G', 0, 0, aapp, one,
626  $ m, 1, a( 1, p ), lda,
627  $ ierr )
628  temp1 = -aapq*d( q ) / d( p )
629  CALL daxpy( m, temp1, work, 1,
630  $ a( 1, p ), 1 )
631  CALL dlascl( 'G', 0, 0, one, aapp,
632  $ m, 1, a( 1, p ), lda,
633  $ ierr )
634  sva( p ) = aapp*dsqrt( max( zero,
635  $ one-aapq*aapq ) )
636  mxsinj = max( mxsinj, sfmin )
637  END IF
638  END IF
639 * END IF ROTOK THEN ... ELSE
640 *
641 * In the case of cancellation in updating SVA(q)
642 * .. recompute SVA(q)
643  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
644  $ THEN
645  IF( ( aaqq.LT.rootbig ) .AND.
646  $ ( aaqq.GT.rootsfmin ) ) THEN
647  sva( q ) = dnrm2( m, a( 1, q ), 1 )*
648  $ d( q )
649  ELSE
650  t = zero
651  aaqq = one
652  CALL dlassq( m, a( 1, q ), 1, t,
653  $ aaqq )
654  sva( q ) = t*dsqrt( aaqq )*d( q )
655  END IF
656  END IF
657  IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
658  IF( ( aapp.LT.rootbig ) .AND.
659  $ ( aapp.GT.rootsfmin ) ) THEN
660  aapp = dnrm2( m, a( 1, p ), 1 )*
661  $ d( p )
662  ELSE
663  t = zero
664  aapp = one
665  CALL dlassq( m, a( 1, p ), 1, t,
666  $ aapp )
667  aapp = t*dsqrt( aapp )*d( p )
668  END IF
669  sva( p ) = aapp
670  END IF
671 * end of OK rotation
672  ELSE
673  notrot = notrot + 1
674 * SKIPPED = SKIPPED + 1
675  pskipped = pskipped + 1
676  ijblsk = ijblsk + 1
677  END IF
678  ELSE
679  notrot = notrot + 1
680  pskipped = pskipped + 1
681  ijblsk = ijblsk + 1
682  END IF
683 
684 * IF ( NOTROT .GE. EMPTSW ) GO TO 2011
685  IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
686  $ THEN
687  sva( p ) = aapp
688  notrot = 0
689  GO TO 2011
690  END IF
691  IF( ( i.LE.swband ) .AND.
692  $ ( pskipped.GT.rowskip ) ) THEN
693  aapp = -aapp
694  notrot = 0
695  GO TO 2203
696  END IF
697 
698 *
699  2200 CONTINUE
700 * end of the q-loop
701  2203 CONTINUE
702 
703  sva( p ) = aapp
704 *
705  ELSE
706  IF( aapp.EQ.zero )notrot = notrot +
707  $ min( jgl+kbl-1, n ) - jgl + 1
708  IF( aapp.LT.zero )notrot = 0
709 *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
710  END IF
711 
712  2100 CONTINUE
713 * end of the p-loop
714  2010 CONTINUE
715 * end of the jbc-loop
716  2011 CONTINUE
717 *2011 bailed out of the jbc-loop
718  DO 2012 p = igl, min( igl+kbl-1, n )
719  sva( p ) = dabs( sva( p ) )
720  2012 CONTINUE
721 *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
722  2000 CONTINUE
723 *2000 :: end of the ibr-loop
724 *
725 * .. update SVA(N)
726  IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
727  $ THEN
728  sva( n ) = dnrm2( m, a( 1, n ), 1 )*d( n )
729  ELSE
730  t = zero
731  aapp = one
732  CALL dlassq( m, a( 1, n ), 1, t, aapp )
733  sva( n ) = t*dsqrt( aapp )*d( n )
734  END IF
735 *
736 * Additional steering devices
737 *
738  IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
739  $ ( iswrot.LE.n ) ) )swband = i
740 
741  IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
742  $ ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
743  GO TO 1994
744  END IF
745 
746 *
747  IF( notrot.GE.emptsw )GO TO 1994
748 
749  1993 CONTINUE
750 * end i=1:NSWEEP loop
751 * #:) Reaching this point means that the procedure has completed the given
752 * number of sweeps.
753  info = nsweep - 1
754  GO TO 1995
755  1994 CONTINUE
756 * #:) Reaching this point means that during the i-th sweep all pivots were
757 * below the given threshold, causing early exit.
758 
759  info = 0
760 * #:) INFO = 0 confirms successful iterations.
761  1995 CONTINUE
762 *
763 * Sort the vector D
764 *
765  DO 5991 p = 1, n - 1
766  q = idamax( n-p+1, sva( p ), 1 ) + p - 1
767  IF( p.NE.q ) THEN
768  temp1 = sva( p )
769  sva( p ) = sva( q )
770  sva( q ) = temp1
771  temp1 = d( p )
772  d( p ) = d( q )
773  d( q ) = temp1
774  CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
775  IF( rsvec )CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
776  END IF
777  5991 CONTINUE
778 *
779  RETURN
780 * ..
781 * .. END OF DGSVJ1
782 * ..
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
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:54
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
subroutine drotm(N, DX, INCX, DY, INCY, DPARAM)
DROTM
Definition: drotm.f:100
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
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: