LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 postmultiplying the N-by-N array V.
                (See the description of V.)
          = 'A': the product of the Jacobi rotations is accumulated
                 by postmultiplying 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-multiplied 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-multiplied by a
                           sequence of Jacobi rotations.
          If JOBV = 'A' then MV rows of V are post-multiplied 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 xerbla(srname, info)
Definition cblat2.f:3285
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
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
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
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
subroutine slassq(n, x, incx, scale, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:124
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
subroutine srotm(n, sx, incx, sy, incy, sparam)
SROTM
Definition srotm.f:97
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: