OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dgsvj1.f
Go to the documentation of this file.
1*> \brief \b DGSVJ1 pre-processor for the routine dgesvj, 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 DGSVJ1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
22* EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* DOUBLE PRECISION EPS, SFMIN, TOL
26* INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
27* CHARACTER*1 JOBV
28* ..
29* .. Array Arguments ..
30* DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
31* $ WORK( LWORK )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> DGSVJ1 is called from DGESVJ as a pre-processor and that is its main
41*> purpose. It applies Jacobi rotations in the same way as DGESVJ 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*> DGSVJ1 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
176*> EPS = DLAMCH('Epsilon')
177*> \endverbatim
178*>
179*> \param[in] SFMIN
180*> \verbatim
181*> SFMIN is DOUBLE PRECISION
182*> SFMIN = DLAMCH('Safe Minimum')
183*> \endverbatim
184*>
185*> \param[in] TOL
186*> \verbatim
187*> TOL is DOUBLE PRECISION
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 DABS(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 DOUBLE PRECISION 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 doubleOTHERcomputational
227*
228*> \par Contributors:
229* ==================
230*>
231*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
232*
233* =====================================================================
234 SUBROUTINE dgsvj1( 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 DOUBLE PRECISION EPS, SFMIN, TOL
243 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
244 CHARACTER*1 JOBV
245* ..
246* .. Array Arguments ..
247 DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
248 $ work( lwork )
249* ..
250*
251* =====================================================================
252*
253* .. Local Parameters ..
254 DOUBLE PRECISION ZERO, HALF, ONE
255 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
256* ..
257* .. Local Scalars ..
258 DOUBLE PRECISION 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 DOUBLE PRECISION FASTR( 5 )
269* ..
270* .. Intrinsic Functions ..
271 INTRINSIC dabs, max, dble, min, dsign, dsqrt
272* ..
273* .. External Functions ..
274 DOUBLE PRECISION DDOT, DNRM2
275 INTEGER IDAMAX
276 LOGICAL LSAME
277 EXTERNAL idamax, lsame, ddot, dnrm2
278* ..
279* .. External Subroutines ..
280 EXTERNAL daxpy, dcopy, dlascl, dlassq, drotm, dswap,
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( 'DGSVJ1', -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 = dsqrt( eps )
328 rootsfmin = dsqrt( sfmin )
329 small = sfmin / eps
330 big = one / sfmin
331 rootbig = one / rootsfmin
332 large = big / dsqrt( dble( m*n ) )
333 bigtheta = one / rooteps
334 roottol = dsqrt( 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 = ( ddot( m, a( 1, p ), 1, a( 1,
427 $ q ), 1 )*d( p )*d( q ) / aaqq )
428 $ / aapp
429 ELSE
430 CALL dcopy( m, a( 1, p ), 1, work, 1 )
431 CALL dlascl( 'g', 0, 0, AAPP, D( p ),
432 $ M, 1, WORK, LDA, IERR )
433 AAPQ = DDOT( 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 = ( DDOT( M, A( 1, p ), 1, A( 1,
444 $ q ), 1 )*D( p )*D( q ) / AAQQ )
445 $ / AAPP
446 ELSE
447 CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
448 CALL DLASCL( 'g', 0, 0, AAQQ, D( q ),
449 $ M, 1, WORK, LDA, IERR )
450 AAPQ = DDOT( M, WORK, 1, A( 1, p ),
451 $ 1 )*D( p ) / AAPP
452 END IF
453 END IF
454
455 MXAAPQ = MAX( MXAAPQ, DABS( AAPQ ) )
456
457* TO rotate or NOT to rotate, THAT is the question ...
458*
459.GT. IF( DABS( 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*DABS(AQOAP-APOAQ) / AAPQ
470.GT. IF( AAQQAAPP0 )THETA = -THETA
471
472.GT. IF( DABS( 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 DROTM( M, A( 1, p ), 1,
477 $ A( 1, q ), 1, FASTR )
478 IF( RSVEC )CALL DROTM( MVL,
479 $ V( 1, p ), 1,
480 $ V( 1, q ), 1,
481 $ FASTR )
482 SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
483 $ ONE+T*APOAQ*AAPQ ) )
484 AAPP = AAPP*DSQRT( MAX( ZERO,
485 $ ONE-T*AQOAP*AAPQ ) )
486 MXSINJ = MAX( MXSINJ, DABS( T ) )
487 ELSE
488*
489* .. choose correct signum for THETA and rotate
490*
491 THSIGN = -DSIGN( ONE, AAPQ )
492.GT. IF( AAQQAAPP0 )THSIGN = -THSIGN
493 T = ONE / ( THETA+THSIGN*
494 $ DSQRT( ONE+THETA*THETA ) )
495 CS = DSQRT( ONE / ( ONE+T*T ) )
496 SN = T*CS
497 MXSINJ = MAX( MXSINJ, DABS( SN ) )
498 SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
499 $ ONE+T*APOAQ*AAPQ ) )
500 AAPP = AAPP*DSQRT( 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 DROTM( M, A( 1, p ), 1,
513 $ A( 1, q ), 1,
514 $ FASTR )
515 IF( RSVEC )CALL DROTM( MVL,
516 $ V( 1, p ), 1, V( 1, q ),
517 $ 1, FASTR )
518 ELSE
519 CALL DAXPY( M, -T*AQOAP,
520 $ A( 1, q ), 1,
521 $ A( 1, p ), 1 )
522 CALL DAXPY( M, CS*SN*APOAQ,
523 $ A( 1, p ), 1,
524 $ A( 1, q ), 1 )
525 IF( RSVEC ) THEN
526 CALL DAXPY( MVL, -T*AQOAP,
527 $ V( 1, q ), 1,
528 $ V( 1, p ), 1 )
529 CALL DAXPY( 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 DAXPY( M, T*APOAQ,
540 $ A( 1, p ), 1,
541 $ A( 1, q ), 1 )
542 CALL DAXPY( M, -CS*SN*AQOAP,
543 $ A( 1, q ), 1,
544 $ A( 1, p ), 1 )
545 IF( RSVEC ) THEN
546 CALL DAXPY( MVL, T*APOAQ,
547 $ V( 1, p ), 1,
548 $ V( 1, q ), 1 )
549 CALL DAXPY( 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 DAXPY( M, -T*AQOAP,
559 $ A( 1, q ), 1,
560 $ A( 1, p ), 1 )
561 CALL DAXPY( 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 DAXPY( MVL,
568 $ -T*AQOAP,
569 $ V( 1, q ), 1,
570 $ V( 1, p ), 1 )
571 CALL DAXPY( MVL,
572 $ CS*SN*APOAQ,
573 $ V( 1, p ), 1,
574 $ V( 1, q ), 1 )
575 END IF
576 ELSE
577 CALL DAXPY( M, T*APOAQ,
578 $ A( 1, p ), 1,
579 $ A( 1, q ), 1 )
580 CALL DAXPY( 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 DAXPY( MVL,
588 $ T*APOAQ, V( 1, p ),
589 $ 1, V( 1, q ), 1 )
590 CALL DAXPY( 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 DCOPY( M, A( 1, p ), 1, WORK,
603 $ 1 )
604 CALL DLASCL( 'g', 0, 0, AAPP, ONE,
605 $ M, 1, WORK, LDA, IERR )
606 CALL DLASCL( 'g', 0, 0, AAQQ, ONE,
607 $ M, 1, A( 1, q ), LDA,
608 $ IERR )
609 TEMP1 = -AAPQ*D( p ) / D( q )
610 CALL DAXPY( M, TEMP1, WORK, 1,
611 $ A( 1, q ), 1 )
612 CALL DLASCL( 'g', 0, 0, ONE, AAQQ,
613 $ M, 1, A( 1, q ), LDA,
614 $ IERR )
615 SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
616 $ ONE-AAPQ*AAPQ ) )
617 MXSINJ = MAX( MXSINJ, SFMIN )
618 ELSE
619 CALL DCOPY( M, A( 1, q ), 1, WORK,
620 $ 1 )
621 CALL DLASCL( 'g', 0, 0, AAQQ, ONE,
622 $ M, 1, WORK, LDA, IERR )
623 CALL DLASCL( 'g', 0, 0, AAPP, ONE,
624 $ M, 1, A( 1, p ), LDA,
625 $ IERR )
626 TEMP1 = -AAPQ*D( q ) / D( p )
627 CALL DAXPY( M, TEMP1, WORK, 1,
628 $ A( 1, p ), 1 )
629 CALL DLASCL( 'g', 0, 0, ONE, AAPP,
630 $ M, 1, A( 1, p ), LDA,
631 $ IERR )
632 SVA( p ) = AAPP*DSQRT( 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.LE. IF( ( SVA( q ) / AAQQ )**2ROOTEPS )
642 $ THEN
643.LT..AND. IF( ( AAQQROOTBIG )
644.GT. $ ( AAQQROOTSFMIN ) ) THEN
645 SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
646 $ D( q )
647 ELSE
648 T = ZERO
649 AAQQ = ONE
650 CALL DLASSQ( M, A( 1, q ), 1, T,
651 $ AAQQ )
652 SVA( q ) = T*DSQRT( AAQQ )*D( q )
653 END IF
654 END IF
655.LE. IF( ( AAPP / AAPP0 )**2ROOTEPS ) THEN
656.LT..AND. IF( ( AAPPROOTBIG )
657.GT. $ ( AAPPROOTSFMIN ) ) THEN
658 AAPP = DNRM2( M, A( 1, p ), 1 )*
659 $ D( p )
660 ELSE
661 T = ZERO
662 AAPP = ONE
663 CALL DLASSQ( M, A( 1, p ), 1, T,
664 $ AAPP )
665 AAPP = T*DSQRT( 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.LE..AND..GE. IF( ( iSWBAND ) ( IJBLSKBLSKIP ) )
684 $ THEN
685 SVA( p ) = AAPP
686 NOTROT = 0
687 GO TO 2011
688 END IF
689.LE..AND. IF( ( iSWBAND )
690.GT. $ ( PSKIPPEDROWSKIP ) ) 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.EQ. IF( AAPPZERO )NOTROT = NOTROT +
705 $ MIN( jgl+KBL-1, N ) - jgl + 1
706.LT. IF( AAPPZERO )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 ) = DABS( 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.LT..AND..GT. IF( ( SVA( N )ROOTBIG ) ( SVA( N )ROOTSFMIN ) )
725 $ THEN
726 SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
727 ELSE
728 T = ZERO
729 AAPP = ONE
730 CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
731 SVA( N ) = T*DSQRT( AAPP )*D( N )
732 END IF
733*
734* Additional steering devices
735*
736.LT..AND..LE..OR. IF( ( iSWBAND ) ( ( MXAAPQROOTTOL )
737.LE. $ ( ISWROTN ) ) )SWBAND = i
738
739.GT..AND..LT..AND. IF( ( iSWBAND+1 ) ( MXAAPQDBLE( N )*TOL )
740.LT. $ ( DBLE( N )*MXAAPQ*MXSINJTOL ) ) THEN
741 GO TO 1994
742 END IF
743
744*
745.GE. IF( NOTROTEMPTSW )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 = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
765.NE. IF( pq ) 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 DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
773 IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
774 END IF
775 5991 CONTINUE
776*
777 RETURN
778* ..
779* .. END OF DGSVJ1
780* ..
781 END
dlassq
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:137
dlascl
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:143
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dgsvj1
subroutine dgsvj1(jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
DGSVJ1 pre-processor for the routine dgesvj, applies Jacobi rotations targeting only particular pivot...
Definition dgsvj1.f:236
drotm
subroutine drotm(n, dx, incx, dy, incy, dparam)
DROTM
Definition drotm.f:96
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21