OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sgsvj1.f
Go to the documentation of this file.
1*> \brief \b SGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivots.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGSVJ1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgsvj1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgsvj1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgsvj1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
22* EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* REAL EPS, SFMIN, TOL
26* INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
27* CHARACTER*1 JOBV
28* ..
29* .. Array Arguments ..
30* REAL A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
31* $ WORK( LWORK )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> SGSVJ1 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 targets only particular pivots and it does not check convergence
43*> (stopping criterion). Few tuning parameters (marked by [TP]) are
44*> available for the implementer.
45*>
46*> Further Details
47*> ~~~~~~~~~~~~~~~
48*> SGSVJ1 applies few sweeps of Jacobi rotations in the column space of
49*> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
50*> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
51*> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
52*> [x]'s in the following scheme:
53*>
54*> | * * * [x] [x] [x]|
55*> | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
56*> | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
57*> |[x] [x] [x] * * * |
58*> |[x] [x] [x] * * * |
59*> |[x] [x] [x] * * * |
60*>
61*> In terms of the columns of A, the first N1 columns are rotated 'against'
62*> the remaining N-N1 columns, trying to increase the angle between the
63*> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
64*> tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter.
65*> The number of sweeps is given in NSWEEP and the orthogonality threshold
66*> is given in TOL.
67*> \endverbatim
68*
69* Arguments:
70* ==========
71*
72*> \param[in] JOBV
73*> \verbatim
74*> JOBV is CHARACTER*1
75*> Specifies whether the output from this procedure is used
76*> to compute the matrix V:
77*> = 'V': the product of the Jacobi rotations is accumulated
78*> by postmulyiplying the N-by-N array V.
79*> (See the description of V.)
80*> = 'A': the product of the Jacobi rotations is accumulated
81*> by postmulyiplying the MV-by-N array V.
82*> (See the descriptions of MV and V.)
83*> = 'N': the Jacobi rotations are not accumulated.
84*> \endverbatim
85*>
86*> \param[in] M
87*> \verbatim
88*> M is INTEGER
89*> The number of rows of the input matrix A. M >= 0.
90*> \endverbatim
91*>
92*> \param[in] N
93*> \verbatim
94*> N is INTEGER
95*> The number of columns of the input matrix A.
96*> M >= N >= 0.
97*> \endverbatim
98*>
99*> \param[in] N1
100*> \verbatim
101*> N1 is INTEGER
102*> N1 specifies the 2 x 2 block partition, the first N1 columns are
103*> rotated 'against' the remaining N-N1 columns of A.
104*> \endverbatim
105*>
106*> \param[in,out] A
107*> \verbatim
108*> A is REAL array, dimension (LDA,N)
109*> On entry, M-by-N matrix A, such that A*diag(D) represents
110*> the input matrix.
111*> On exit,
112*> A_onexit * D_onexit represents the input matrix A*diag(D)
113*> post-multiplied by a sequence of Jacobi rotations, where the
114*> rotation threshold and the total number of sweeps are given in
115*> TOL and NSWEEP, respectively.
116*> (See the descriptions of N1, D, TOL and NSWEEP.)
117*> \endverbatim
118*>
119*> \param[in] LDA
120*> \verbatim
121*> LDA is INTEGER
122*> The leading dimension of the array A. LDA >= max(1,M).
123*> \endverbatim
124*>
125*> \param[in,out] D
126*> \verbatim
127*> D is REAL array, dimension (N)
128*> The array D accumulates the scaling factors from the fast scaled
129*> Jacobi rotations.
130*> On entry, A*diag(D) represents the input matrix.
131*> On exit, A_onexit*diag(D_onexit) represents the input matrix
132*> post-multiplied by a sequence of Jacobi rotations, where the
133*> rotation threshold and the total number of sweeps are given in
134*> TOL and NSWEEP, respectively.
135*> (See the descriptions of N1, A, TOL and NSWEEP.)
136*> \endverbatim
137*>
138*> \param[in,out] SVA
139*> \verbatim
140*> SVA is REAL array, dimension (N)
141*> On entry, SVA contains the Euclidean norms of the columns of
142*> the matrix A*diag(D).
143*> On exit, SVA contains the Euclidean norms of the columns of
144*> the matrix onexit*diag(D_onexit).
145*> \endverbatim
146*>
147*> \param[in] MV
148*> \verbatim
149*> MV is INTEGER
150*> If JOBV = 'A', then MV rows of V are post-multipled by a
151*> sequence of Jacobi rotations.
152*> If JOBV = 'N', then MV is not referenced.
153*> \endverbatim
154*>
155*> \param[in,out] V
156*> \verbatim
157*> V is REAL array, dimension (LDV,N)
158*> If JOBV = 'V' then N rows of V are post-multipled by a
159*> sequence of Jacobi rotations.
160*> If JOBV = 'A' then MV rows of V are post-multipled by a
161*> sequence of Jacobi rotations.
162*> If JOBV = 'N', then V is not referenced.
163*> \endverbatim
164*>
165*> \param[in] LDV
166*> \verbatim
167*> LDV is INTEGER
168*> The leading dimension of the array V, LDV >= 1.
169*> If JOBV = 'V', LDV >= N.
170*> If JOBV = 'A', LDV >= MV.
171*> \endverbatim
172*>
173*> \param[in] EPS
174*> \verbatim
175*> EPS is REAL
176*> EPS = SLAMCH('Epsilon')
177*> \endverbatim
178*>
179*> \param[in] SFMIN
180*> \verbatim
181*> SFMIN is REAL
182*> SFMIN = SLAMCH('Safe Minimum')
183*> \endverbatim
184*>
185*> \param[in] TOL
186*> \verbatim
187*> TOL is REAL
188*> TOL is the threshold for Jacobi rotations. For a pair
189*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
190*> applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.
191*> \endverbatim
192*>
193*> \param[in] NSWEEP
194*> \verbatim
195*> NSWEEP is INTEGER
196*> NSWEEP is the number of sweeps of Jacobi rotations to be
197*> performed.
198*> \endverbatim
199*>
200*> \param[out] WORK
201*> \verbatim
202*> WORK is REAL array, dimension (LWORK)
203*> \endverbatim
204*>
205*> \param[in] LWORK
206*> \verbatim
207*> LWORK is INTEGER
208*> LWORK is the dimension of WORK. LWORK >= M.
209*> \endverbatim
210*>
211*> \param[out] INFO
212*> \verbatim
213*> INFO is INTEGER
214*> = 0: successful exit.
215*> < 0: if INFO = -i, then the i-th argument had an illegal value
216*> \endverbatim
217*
218* Authors:
219* ========
220*
221*> \author Univ. of Tennessee
222*> \author Univ. of California Berkeley
223*> \author Univ. of Colorado Denver
224*> \author NAG Ltd.
225*
226*> \ingroup realOTHERcomputational
227*
228*> \par Contributors:
229* ==================
230*>
231*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
232*
233* =====================================================================
234 SUBROUTINE sgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
235 $ EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
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.NOT..OR..OR. IF( ( RSVEC APPLV LSAME( JOBV, 'n' ) ) ) THEN
290 INFO = -1
291.LT. ELSE IF( M0 ) THEN
292 INFO = -2
293.LT..OR..GT. ELSE IF( ( N0 ) ( NM ) ) THEN
294 INFO = -3
295.LT. ELSE IF( N10 ) THEN
296 INFO = -4
297.LT. ELSE IF( LDAM ) THEN
298 INFO = -6
299.OR..AND..LT. ELSE IF( ( RSVECAPPLV ) ( MV0 ) ) THEN
300 INFO = -9
301.AND..LT..OR. ELSE IF( ( RSVEC( LDVN ) )
302.AND..LT. $ ( APPLV( LDVMV ) ) ) THEN
303 INFO = -11
304.LE. ELSE IF( TOLEPS ) THEN
305 INFO = -14
306.LT. ELSE IF( NSWEEP0 ) THEN
307 INFO = -15
308.LT. ELSE IF( LWORKM ) THEN
309 INFO = -17
310 ELSE
311 INFO = 0
312 END IF
313*
314* #:(
315.NE. IF( INFO0 ) 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.OR. RSVEC = RSVEC 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.NE. IF( ( NBLR*KBL )N1 )NBLR = NBLR + 1
349
350* .. the tiling is nblr-by-nblc [tiles]
351
352 NBLC = ( N-N1 ) / KBL
353.NE. IF( ( NBLC*KBL )( 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.GT. IF( AAPPZERO ) THEN
405
406 PSKIPPED = 0
407
408 DO 2200 q = jgl, MIN( jgl+KBL-1, N )
409*
410 AAQQ = SVA( q )
411
412.GT. IF( AAQQZERO ) THEN
413 AAPP0 = AAPP
414*
415* .. M x 2 Jacobi SVD ..
416*
417* .. Safe Gram matrix computation ..
418*
419.GE. IF( AAQQONE ) THEN
420.GE. IF( AAPPAAQQ ) THEN
421.LE. ROTOK = ( SMALL*AAPP )AAQQ
422 ELSE
423.LE. ROTOK = ( SMALL*AAQQ )AAPP
424 END IF
425.LT. IF( AAPP( 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.GE. IF( AAPPAAQQ ) THEN
438.LE. ROTOK = AAPP( AAQQ / SMALL )
439 ELSE
440.LE. ROTOK = AAQQ( AAPP / SMALL )
441 END IF
442.GT. IF( AAPP( 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.GT. IF( ABS( AAPQ )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.GT. IF( AAQQAAPP0 )THETA = -THETA
471
472.GT. IF( ABS( THETA )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.GT. IF( AAQQAAPP0 )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.GE. IF( D( p )ONE ) THEN
506*
507.GE. IF( D( q )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.GE. IF( D( q )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.GE. IF( D( p )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.GT. IF( AAPPAAQQ ) 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* ..
781 END
slassq
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:137
slascl
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
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
sgsvj1
subroutine sgsvj1(jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
SGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivot...
Definition sgsvj1.f:236
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
srotm
subroutine srotm(n, sx, incx, sy, incy, sparam)
SROTM
Definition srotm.f:97
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21