3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SGSVJ0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgsvj0.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgsvj0.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgsvj0.f">
21 * SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
22 * SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
26 * REAL EPS, SFMIN, TOL
29 * .. Array Arguments ..
30 * REAL A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
40 *> SGSVJ0 is called from SGESVJ as a pre-processor and that is its main
41 *> purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
42 *> it does not check convergence (stopping criterion). Few tuning
43 *> parameters (marked by [TP]) are available for the implementer.
51 *> JOBV is CHARACTER*1
52 *> Specifies whether the output from this procedure is used
53 *> to compute the matrix V:
54 *> = 'V': the product of the Jacobi rotations is accumulated
55 *> by postmulyiplying the N-by-N array V.
56 *> (See the description of V.)
57 *> = 'A': the product of the Jacobi rotations is accumulated
58 *> by postmulyiplying the MV-by-N array V.
59 *> (See the descriptions of MV and V.)
60 *> = 'N': the Jacobi rotations are not accumulated.
66 *> The number of rows of the input matrix A. M >= 0.
72 *> The number of columns of the input matrix A.
78 *> A is REAL array, dimension (LDA,N)
79 *> On entry, M-by-N matrix A, such that A*diag(D) represents
82 *> A_onexit * D_onexit represents the input matrix A*diag(D)
83 *> post-multiplied by a sequence of Jacobi rotations, where the
84 *> rotation threshold and the total number of sweeps are given in
85 *> TOL and NSWEEP, respectively.
86 *> (See the descriptions of D, TOL and NSWEEP.)
92 *> The leading dimension of the array A. LDA >= max(1,M).
97 *> D is REAL array, dimension (N)
98 *> The array D accumulates the scaling factors from the fast scaled
100 *> On entry, A*diag(D) represents the input matrix.
101 *> On exit, A_onexit*diag(D_onexit) represents the input matrix
102 *> post-multiplied by a sequence of Jacobi rotations, where the
103 *> rotation threshold and the total number of sweeps are given in
104 *> TOL and NSWEEP, respectively.
105 *> (See the descriptions of A, TOL and NSWEEP.)
108 *> \param[in,out] SVA
110 *> SVA is REAL array, dimension (N)
111 *> On entry, SVA contains the Euclidean norms of the columns of
112 *> the matrix A*diag(D).
113 *> On exit, SVA contains the Euclidean norms of the columns of
114 *> the matrix onexit*diag(D_onexit).
120 *> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
121 *> sequence of Jacobi rotations.
122 *> If JOBV = 'N', then MV is not referenced.
127 *> V is REAL array, dimension (LDV,N)
128 *> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
129 *> sequence of Jacobi rotations.
130 *> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
131 *> sequence of Jacobi rotations.
132 *> If JOBV = 'N', then V is not referenced.
138 *> The leading dimension of the array V, LDV >= 1.
139 *> If JOBV = 'V', LDV .GE. N.
140 *> If JOBV = 'A', LDV .GE. MV.
146 *> EPS = SLAMCH('Epsilon')
152 *> SFMIN = SLAMCH('Safe Minimum')
158 *> TOL is the threshold for Jacobi rotations. For a pair
159 *> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
160 *> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
166 *> NSWEEP is the number of sweeps of Jacobi rotations to be
172 *> WORK is REAL array, dimension LWORK.
178 *> LWORK is the dimension of WORK. LWORK .GE. M.
184 *> = 0 : successful exit.
185 *> < 0 : if INFO = -i, then the i-th argument had an illegal value
191 *> \author Univ. of Tennessee
192 *> \author Univ. of California Berkeley
193 *> \author Univ. of Colorado Denver
196 *> \date November 2011
198 *> \ingroup realOTHERcomputational
200 *> \par Further Details:
201 * =====================
203 *> SGSVJ0 is used just to enable SGESVJ to call a simplified version of
204 *> itself to work on a submatrix of the original matrix.
206 *> \par Contributors:
209 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
211 *> \par Bugs, Examples and Comments:
212 * =================================
214 *> Please report all bugs and send interesting test examples and comments to
215 *> drmac@math.hr. Thank you.
217 * =====================================================================
218 SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
219 $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
221 * -- LAPACK computational routine (version 1.23, October 23. 2008.) --
222 * -- LAPACK is a software package provided by Univ. of Tennessee, --
223 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
226 * .. Scalar Arguments ..
227 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
231 * .. Array Arguments ..
232 REAL A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
236 * =====================================================================
238 * .. Local Parameters ..
239 REAL ZERO, HALF, ONE, TWO
240 PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0,
243 * .. Local Scalars ..
244 REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
245 $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
246 $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
248 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
249 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
250 $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
251 LOGICAL APPLV, ROTOK, RSVEC
256 * .. Intrinsic Functions ..
257 INTRINSIC ABS, AMAX1, FLOAT, MIN0, SIGN, SQRT
259 * .. External Functions ..
263 EXTERNAL ISAMAX, LSAME, SDOT, SNRM2
265 * .. External Subroutines ..
266 EXTERNAL SAXPY, SCOPY, SLASCL, SLASSQ, SROTM, SSWAP
268 * .. Executable Statements ..
270 * Test the input parameters.
272 APPLV = LSAME( JOBV, 'A' )
273 RSVEC = LSAME( JOBV, 'V' )
274 IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
276 ELSE IF( M.LT.0 ) THEN
278 ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
280 ELSE IF( LDA.LT.M ) THEN
282 ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
284 ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
285 $ ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
287 ELSE IF( TOL.LE.EPS ) THEN
289 ELSE IF( NSWEEP.LT.0 ) THEN
291 ELSE IF( LWORK.LT.M ) THEN
299 CALL XERBLA( 'SGSVJ0', -INFO )
305 ELSE IF( APPLV ) THEN
308 RSVEC = RSVEC .OR. APPLV
310 ROOTEPS = SQRT( EPS )
311 ROOTSFMIN = SQRT( SFMIN )
314 ROOTBIG = ONE / ROOTSFMIN
315 BIGTHETA = ONE / ROOTEPS
316 ROOTTOL = SQRT( TOL )
318 * .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
320 EMPTSW = ( N*( N-1 ) ) / 2
324 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
328 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
329 * if SGESVJ is used as a computational routine in the preconditioned
330 * Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
334 *[TP] KBL is a tuning parameter that defines the tile size in the
335 * tiling of the p-q loops of pivot pairs. In general, an optimal
336 * value of KBL depends on the matrix dimensions and on the
337 * parameters of the computer's memory.
340 IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
342 BLSKIP = ( KBL**2 ) + 1
343 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
345 ROWSKIP = MIN0( 5, KBL )
346 *[TP] ROWSKIP is a tuning parameter.
349 *[TP] LKAHEAD is a tuning parameter.
353 DO 1993 i = 1, NSWEEP
365 igl = ( ibr-1 )*KBL + 1
367 DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
371 DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
373 * .. de Rijk's pivoting
374 q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
376 CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
377 IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1,
389 * Column norms are periodically updated by explicit
392 * Some BLAS implementations compute SNRM2(M,A(1,p),1)
393 * as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may result in
394 * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and
395 * undeflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
396 * Hence, SNRM2 cannot be trusted, not even in the case when
397 * the true norm is far from the under(over)flow boundaries.
398 * If properly implemented SNRM2 is available, the IF-THEN-ELSE
399 * below should read "AAPP = SNRM2( M, A(1,p), 1 ) * D(p)".
401 IF( ( SVA( p ).LT.ROOTBIG ) .AND.
402 $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN
403 SVA( p ) = SNRM2( M, A( 1, p ), 1 )*D( p )
407 CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
408 SVA( p ) = TEMP1*SQRT( AAPP )*D( p )
416 IF( AAPP.GT.ZERO ) THEN
420 DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
424 IF( AAQQ.GT.ZERO ) THEN
427 IF( AAQQ.GE.ONE ) THEN
428 ROTOK = ( SMALL*AAPP ).LE.AAQQ
429 IF( AAPP.LT.( BIG / AAQQ ) ) THEN
430 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
431 $ q ), 1 )*D( p )*D( q ) / AAQQ )
434 CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
435 CALL SLASCL( 'G', 0, 0, AAPP, D( p ),
436 $ M, 1, WORK, LDA, IERR )
437 AAPQ = SDOT( M, WORK, 1, A( 1, q ),
441 ROTOK = AAPP.LE.( AAQQ / SMALL )
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 )
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 ),
455 MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
457 * TO rotate or NOT to rotate, THAT is the question ...
459 IF( ABS( AAPQ ).GT.TOL ) THEN
462 * ROTATED = ROTATED + ONE
474 THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
476 IF( ABS( THETA ).GT.BIGTHETA ) THEN
479 FASTR( 3 ) = T*D( p ) / D( q )
480 FASTR( 4 ) = -T*D( q ) / D( p )
481 CALL SROTM( M, A( 1, p ), 1,
482 $ A( 1, q ), 1, FASTR )
483 IF( RSVEC )CALL SROTM( MVL,
487 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
488 $ ONE+T*APOAQ*AAPQ ) )
489 AAPP = AAPP*SQRT( AMAX1( ZERO,
490 $ ONE-T*AQOAP*AAPQ ) )
491 MXSINJ = AMAX1( MXSINJ, ABS( T ) )
495 * .. choose correct signum for THETA and rotate
497 THSIGN = -SIGN( ONE, AAPQ )
498 T = ONE / ( THETA+THSIGN*
499 $ SQRT( ONE+THETA*THETA ) )
500 CS = SQRT( ONE / ( ONE+T*T ) )
503 MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
504 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
505 $ ONE+T*APOAQ*AAPQ ) )
506 AAPP = AAPP*SQRT( AMAX1( ZERO,
507 $ ONE-T*AQOAP*AAPQ ) )
509 APOAQ = D( p ) / D( q )
510 AQOAP = D( q ) / D( p )
511 IF( D( p ).GE.ONE ) THEN
512 IF( D( q ).GE.ONE ) THEN
514 FASTR( 4 ) = -T*AQOAP
517 CALL SROTM( M, A( 1, p ), 1,
520 IF( RSVEC )CALL SROTM( MVL,
521 $ V( 1, p ), 1, V( 1, q ),
524 CALL SAXPY( M, -T*AQOAP,
527 CALL SAXPY( M, CS*SN*APOAQ,
533 CALL SAXPY( MVL, -T*AQOAP,
543 IF( D( q ).GE.ONE ) THEN
544 CALL SAXPY( M, T*APOAQ,
547 CALL SAXPY( M, -CS*SN*AQOAP,
553 CALL SAXPY( MVL, T*APOAQ,
562 IF( D( p ).GE.D( q ) ) THEN
563 CALL SAXPY( M, -T*AQOAP,
566 CALL SAXPY( M, CS*SN*APOAQ,
582 CALL SAXPY( M, T*APOAQ,
593 $ T*APOAQ, V( 1, p ),
606 * .. have to use modified Gram-Schmidt like transformation
607 CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
608 CALL SLASCL( 'G', 0, 0, AAPP, ONE, M,
609 $ 1, WORK, LDA, IERR )
610 CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M,
611 $ 1, A( 1, q ), LDA, IERR )
612 TEMP1 = -AAPQ*D( p ) / D( q )
613 CALL SAXPY( M, TEMP1, WORK, 1,
615 CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M,
616 $ 1, A( 1, q ), LDA, IERR )
617 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
619 MXSINJ = AMAX1( MXSINJ, SFMIN )
621 * END IF ROTOK THEN ... ELSE
623 * In the case of cancellation in updating SVA(q), SVA(p)
624 * recompute SVA(q), SVA(p).
625 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
627 IF( ( AAQQ.LT.ROOTBIG ) .AND.
628 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
629 SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
634 CALL SLASSQ( M, A( 1, q ), 1, T,
636 SVA( q ) = T*SQRT( AAQQ )*D( q )
639 IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
640 IF( ( AAPP.LT.ROOTBIG ) .AND.
641 $ ( AAPP.GT.ROOTSFMIN ) ) THEN
642 AAPP = SNRM2( M, A( 1, p ), 1 )*
647 CALL SLASSQ( M, A( 1, p ), 1, T,
649 AAPP = T*SQRT( AAPP )*D( p )
655 * A(:,p) and A(:,q) already numerically orthogonal
656 IF( ir1.EQ.0 )NOTROT = NOTROT + 1
657 PSKIPPED = PSKIPPED + 1
660 * A(:,q) is zero column
661 IF( ir1.EQ.0 )NOTROT = NOTROT + 1
662 PSKIPPED = PSKIPPED + 1
665 IF( ( i.LE.SWBAND ) .AND.
666 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
667 IF( ir1.EQ.0 )AAPP = -AAPP
676 * bailed out of q-loop
682 IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
683 $ NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
688 * end of doing the block ( ibr, ibr )
692 *........................................................
693 * ... go to the off diagonal blocks
695 igl = ( ibr-1 )*KBL + 1
697 DO 2010 jbc = ibr + 1, NBL
699 jgl = ( jbc-1 )*KBL + 1
701 * doing the block at ( ibr, jbc )
704 DO 2100 p = igl, MIN0( igl+KBL-1, N )
708 IF( AAPP.GT.ZERO ) THEN
712 DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
716 IF( AAQQ.GT.ZERO ) THEN
719 * .. M x 2 Jacobi SVD ..
721 * .. Safe Gram matrix computation ..
723 IF( AAQQ.GE.ONE ) THEN
724 IF( AAPP.GE.AAQQ ) THEN
725 ROTOK = ( SMALL*AAPP ).LE.AAQQ
727 ROTOK = ( SMALL*AAQQ ).LE.AAPP
729 IF( AAPP.LT.( BIG / AAQQ ) ) THEN
730 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
731 $ q ), 1 )*D( p )*D( q ) / AAQQ )
734 CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
735 CALL SLASCL( 'G', 0, 0, AAPP, D( p ),
736 $ M, 1, WORK, LDA, IERR )
737 AAPQ = SDOT( M, WORK, 1, A( 1, q ),
741 IF( AAPP.GE.AAQQ ) THEN
742 ROTOK = AAPP.LE.( AAQQ / SMALL )
744 ROTOK = AAQQ.LE.( AAPP / SMALL )
746 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
747 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
748 $ q ), 1 )*D( p )*D( q ) / AAQQ )
751 CALL SCOPY( M, A( 1, q ), 1, WORK, 1 )
752 CALL SLASCL( 'G', 0, 0, AAQQ, D( q ),
753 $ M, 1, WORK, LDA, IERR )
754 AAPQ = SDOT( M, WORK, 1, A( 1, p ),
759 MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
761 * TO rotate or NOT to rotate, THAT is the question ...
763 IF( ABS( AAPQ ).GT.TOL ) THEN
765 * ROTATED = ROTATED + 1
773 THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
774 IF( AAQQ.GT.AAPP0 )THETA = -THETA
776 IF( ABS( THETA ).GT.BIGTHETA ) THEN
778 FASTR( 3 ) = T*D( p ) / D( q )
779 FASTR( 4 ) = -T*D( q ) / D( p )
780 CALL SROTM( M, A( 1, p ), 1,
781 $ A( 1, q ), 1, FASTR )
782 IF( RSVEC )CALL SROTM( MVL,
786 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
787 $ ONE+T*APOAQ*AAPQ ) )
788 AAPP = AAPP*SQRT( AMAX1( ZERO,
789 $ ONE-T*AQOAP*AAPQ ) )
790 MXSINJ = AMAX1( MXSINJ, ABS( T ) )
793 * .. choose correct signum for THETA and rotate
795 THSIGN = -SIGN( ONE, AAPQ )
796 IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
797 T = ONE / ( THETA+THSIGN*
798 $ SQRT( ONE+THETA*THETA ) )
799 CS = SQRT( ONE / ( ONE+T*T ) )
801 MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
802 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
803 $ ONE+T*APOAQ*AAPQ ) )
804 AAPP = AAPP*SQRT( AMAX1( ZERO,
805 $ ONE-T*AQOAP*AAPQ ) )
807 APOAQ = D( p ) / D( q )
808 AQOAP = D( q ) / D( p )
809 IF( D( p ).GE.ONE ) THEN
811 IF( D( q ).GE.ONE ) THEN
813 FASTR( 4 ) = -T*AQOAP
816 CALL SROTM( M, A( 1, p ), 1,
819 IF( RSVEC )CALL SROTM( MVL,
820 $ V( 1, p ), 1, V( 1, q ),
823 CALL SAXPY( M, -T*AQOAP,
826 CALL SAXPY( M, CS*SN*APOAQ,
830 CALL SAXPY( MVL, -T*AQOAP,
842 IF( D( q ).GE.ONE ) THEN
843 CALL SAXPY( M, T*APOAQ,
846 CALL SAXPY( M, -CS*SN*AQOAP,
850 CALL SAXPY( MVL, T*APOAQ,
861 IF( D( p ).GE.D( q ) ) THEN
862 CALL SAXPY( M, -T*AQOAP,
865 CALL SAXPY( M, CS*SN*APOAQ,
881 CALL SAXPY( M, T*APOAQ,
892 $ T*APOAQ, V( 1, p ),
905 IF( AAPP.GT.AAQQ ) THEN
906 CALL SCOPY( M, A( 1, p ), 1, WORK,
908 CALL SLASCL( 'G', 0, 0, AAPP, ONE,
909 $ M, 1, WORK, LDA, IERR )
910 CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
911 $ M, 1, A( 1, q ), LDA,
913 TEMP1 = -AAPQ*D( p ) / D( q )
914 CALL SAXPY( M, TEMP1, WORK, 1,
916 CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
917 $ M, 1, A( 1, q ), LDA,
919 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
921 MXSINJ = AMAX1( MXSINJ, SFMIN )
923 CALL SCOPY( M, A( 1, q ), 1, WORK,
925 CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
926 $ M, 1, WORK, LDA, IERR )
927 CALL SLASCL( 'G', 0, 0, AAPP, ONE,
928 $ M, 1, A( 1, p ), LDA,
930 TEMP1 = -AAPQ*D( q ) / D( p )
931 CALL SAXPY( M, TEMP1, WORK, 1,
933 CALL SLASCL( 'G', 0, 0, ONE, AAPP,
934 $ M, 1, A( 1, p ), LDA,
936 SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
938 MXSINJ = AMAX1( MXSINJ, SFMIN )
941 * END IF ROTOK THEN ... ELSE
943 * In the case of cancellation in updating SVA(q)
944 * .. recompute SVA(q)
945 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
947 IF( ( AAQQ.LT.ROOTBIG ) .AND.
948 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
949 SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
954 CALL SLASSQ( M, A( 1, q ), 1, T,
956 SVA( q ) = T*SQRT( AAQQ )*D( q )
959 IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
960 IF( ( AAPP.LT.ROOTBIG ) .AND.
961 $ ( AAPP.GT.ROOTSFMIN ) ) THEN
962 AAPP = SNRM2( M, A( 1, p ), 1 )*
967 CALL SLASSQ( M, A( 1, p ), 1, T,
969 AAPP = T*SQRT( AAPP )*D( p )
976 PSKIPPED = PSKIPPED + 1
981 PSKIPPED = PSKIPPED + 1
985 IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
991 IF( ( i.LE.SWBAND ) .AND.
992 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
1005 IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1006 $ MIN0( jgl+KBL-1, N ) - jgl + 1
1007 IF( AAPP.LT.ZERO )NOTROT = 0
1013 * end of the jbc-loop
1015 *2011 bailed out of the jbc-loop
1016 DO 2012 p = igl, MIN0( igl+KBL-1, N )
1017 SVA( p ) = ABS( SVA( p ) )
1021 *2000 :: end of the ibr-loop
1024 IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
1026 SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
1030 CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
1031 SVA( N ) = T*SQRT( AAPP )*D( N )
1034 * Additional steering devices
1036 IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
1037 $ ( ISWROT.LE.N ) ) )SWBAND = i
1039 IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.FLOAT( N )*TOL ) .AND.
1040 $ ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1044 IF( NOTROT.GE.EMPTSW )GO TO 1994
1047 * end i=1:NSWEEP loop
1048 * #:) Reaching this point means that the procedure has comleted the given
1049 * number of iterations.
1053 * #:) Reaching this point means that during the i-th sweep all pivots were
1054 * below the given tolerance, causing early exit.
1057 * #:) INFO = 0 confirms successful iterations.
1060 * Sort the vector D.
1061 DO 5991 p = 1, N - 1
1062 q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
1070 CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
1071 IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )