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 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 = ( 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 IF( dabs( 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*dabs(aqoap-apoaq) / aapq
470 IF( aaqq.GT.aapp0 )theta = -theta
471
472 IF( dabs( 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 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 IF( aaqq.GT.aapp0 )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 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 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 IF( d( q ).GE.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 IF( d( p ).GE.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 IF( aapp.GT.aaqq ) 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 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
642 $ THEN
643 IF( ( aaqq.LT.rootbig ) .AND.
644 $ ( aaqq.GT.rootsfmin ) ) 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 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
656 IF( ( aapp.LT.rootbig ) .AND.
657 $ ( aapp.GT.rootsfmin ) ) 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 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 ) = 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 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.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 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.dble( n )*tol ) .AND.
740 $ ( dble( 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 = idamax( 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 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