LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
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*> Download DGSVJ1 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj1.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj1.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj1.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
20* EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* DOUBLE PRECISION EPS, SFMIN, TOL
24* INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
25* CHARACTER*1 JOBV
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
29* $ WORK( LWORK )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DGSVJ1 is called from DGESVJ as a pre-processor and that is its main
39*> purpose. It applies Jacobi rotations in the same way as DGESVJ does, but
40*> it targets only particular pivots and it does not check convergence
41*> (stopping criterion). Few tuning parameters (marked by [TP]) are
42*> available for the implementer.
43*>
44*> Further Details
45*> ~~~~~~~~~~~~~~~
46*> DGSVJ1 applies few sweeps of Jacobi rotations in the column space of
47*> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
48*> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
49*> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
50*> [x]'s in the following scheme:
51*>
52*> | * * * [x] [x] [x]|
53*> | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
54*> | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
55*> |[x] [x] [x] * * * |
56*> |[x] [x] [x] * * * |
57*> |[x] [x] [x] * * * |
58*>
59*> In terms of the columns of A, the first N1 columns are rotated 'against'
60*> the remaining N-N1 columns, trying to increase the angle between the
61*> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
62*> tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter.
63*> The number of sweeps is given in NSWEEP and the orthogonality threshold
64*> is given in TOL.
65*> \endverbatim
66*
67* Arguments:
68* ==========
69*
70*> \param[in] JOBV
71*> \verbatim
72*> JOBV is CHARACTER*1
73*> Specifies whether the output from this procedure is used
74*> to compute the matrix V:
75*> = 'V': the product of the Jacobi rotations is accumulated
76*> by postmultiplying the N-by-N array V.
77*> (See the description of V.)
78*> = 'A': the product of the Jacobi rotations is accumulated
79*> by postmultiplying the MV-by-N array V.
80*> (See the descriptions of MV and V.)
81*> = 'N': the Jacobi rotations are not accumulated.
82*> \endverbatim
83*>
84*> \param[in] M
85*> \verbatim
86*> M is INTEGER
87*> The number of rows of the input matrix A. M >= 0.
88*> \endverbatim
89*>
90*> \param[in] N
91*> \verbatim
92*> N is INTEGER
93*> The number of columns of the input matrix A.
94*> M >= N >= 0.
95*> \endverbatim
96*>
97*> \param[in] N1
98*> \verbatim
99*> N1 is INTEGER
100*> N1 specifies the 2 x 2 block partition, the first N1 columns are
101*> rotated 'against' the remaining N-N1 columns of A.
102*> \endverbatim
103*>
104*> \param[in,out] A
105*> \verbatim
106*> A is DOUBLE PRECISION array, dimension (LDA,N)
107*> On entry, M-by-N matrix A, such that A*diag(D) represents
108*> the input matrix.
109*> On exit,
110*> A_onexit * D_onexit represents the input matrix A*diag(D)
111*> post-multiplied by a sequence of Jacobi rotations, where the
112*> rotation threshold and the total number of sweeps are given in
113*> TOL and NSWEEP, respectively.
114*> (See the descriptions of N1, D, TOL and NSWEEP.)
115*> \endverbatim
116*>
117*> \param[in] LDA
118*> \verbatim
119*> LDA is INTEGER
120*> The leading dimension of the array A. LDA >= max(1,M).
121*> \endverbatim
122*>
123*> \param[in,out] D
124*> \verbatim
125*> D is DOUBLE PRECISION array, dimension (N)
126*> The array D accumulates the scaling factors from the fast scaled
127*> Jacobi rotations.
128*> On entry, A*diag(D) represents the input matrix.
129*> On exit, A_onexit*diag(D_onexit) represents the input matrix
130*> post-multiplied by a sequence of Jacobi rotations, where the
131*> rotation threshold and the total number of sweeps are given in
132*> TOL and NSWEEP, respectively.
133*> (See the descriptions of N1, A, TOL and NSWEEP.)
134*> \endverbatim
135*>
136*> \param[in,out] SVA
137*> \verbatim
138*> SVA is DOUBLE PRECISION array, dimension (N)
139*> On entry, SVA contains the Euclidean norms of the columns of
140*> the matrix A*diag(D).
141*> On exit, SVA contains the Euclidean norms of the columns of
142*> the matrix onexit*diag(D_onexit).
143*> \endverbatim
144*>
145*> \param[in] MV
146*> \verbatim
147*> MV is INTEGER
148*> If JOBV = 'A', then MV rows of V are post-multiplied by a
149*> sequence of Jacobi rotations.
150*> If JOBV = 'N', then MV is not referenced.
151*> \endverbatim
152*>
153*> \param[in,out] V
154*> \verbatim
155*> V is DOUBLE PRECISION array, dimension (LDV,N)
156*> If JOBV = 'V', then N rows of V are post-multiplied by a
157*> sequence of Jacobi rotations.
158*> If JOBV = 'A', then MV rows of V are post-multiplied by a
159*> sequence of Jacobi rotations.
160*> If JOBV = 'N', then V is not referenced.
161*> \endverbatim
162*>
163*> \param[in] LDV
164*> \verbatim
165*> LDV is INTEGER
166*> The leading dimension of the array V, LDV >= 1.
167*> If JOBV = 'V', LDV >= N.
168*> If JOBV = 'A', LDV >= MV.
169*> \endverbatim
170*>
171*> \param[in] EPS
172*> \verbatim
173*> EPS is DOUBLE PRECISION
174*> EPS = DLAMCH('Epsilon')
175*> \endverbatim
176*>
177*> \param[in] SFMIN
178*> \verbatim
179*> SFMIN is DOUBLE PRECISION
180*> SFMIN = DLAMCH('Safe Minimum')
181*> \endverbatim
182*>
183*> \param[in] TOL
184*> \verbatim
185*> TOL is DOUBLE PRECISION
186*> TOL is the threshold for Jacobi rotations. For a pair
187*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
188*> applied only if DABS(COS(angle(A(:,p),A(:,q)))) > TOL.
189*> \endverbatim
190*>
191*> \param[in] NSWEEP
192*> \verbatim
193*> NSWEEP is INTEGER
194*> NSWEEP is the number of sweeps of Jacobi rotations to be
195*> performed.
196*> \endverbatim
197*>
198*> \param[out] WORK
199*> \verbatim
200*> WORK is DOUBLE PRECISION array, dimension (LWORK)
201*> \endverbatim
202*>
203*> \param[in] LWORK
204*> \verbatim
205*> LWORK is INTEGER
206*> LWORK is the dimension of WORK. LWORK >= M.
207*> \endverbatim
208*>
209*> \param[out] INFO
210*> \verbatim
211*> INFO is INTEGER
212*> = 0: successful exit.
213*> < 0: if INFO = -i, then the i-th argument had an illegal value
214*> \endverbatim
215*
216* Authors:
217* ========
218*
219*> \author Univ. of Tennessee
220*> \author Univ. of California Berkeley
221*> \author Univ. of Colorado Denver
222*> \author NAG Ltd.
223*
224*> \ingroup gsvj1
225*
226*> \par Contributors:
227* ==================
228*>
229*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
230*
231* =====================================================================
232 SUBROUTINE dgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
233 $ EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
234*
235* -- LAPACK computational routine --
236* -- LAPACK is a software package provided by Univ. of Tennessee, --
237* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
238*
239* .. Scalar Arguments ..
240 DOUBLE PRECISION EPS, SFMIN, TOL
241 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
242 CHARACTER*1 JOBV
243* ..
244* .. Array Arguments ..
245 DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
246 $ work( lwork )
247* ..
248*
249* =====================================================================
250*
251* .. Local Parameters ..
252 DOUBLE PRECISION ZERO, HALF, ONE
253 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
254* ..
255* .. Local Scalars ..
256 DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
257 $ bigtheta, cs, large, mxaapq, mxsinj, rootbig,
258 $ rooteps, rootsfmin, roottol, small, sn, t,
259 $ temp1, theta, thsign
260 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
261 $ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
262 $ p, pskipped, q, rowskip, swband
263 LOGICAL APPLV, ROTOK, RSVEC
264* ..
265* .. Local Arrays ..
266 DOUBLE PRECISION FASTR( 5 )
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC dabs, max, dble, min, dsign, dsqrt
270* ..
271* .. External Functions ..
272 DOUBLE PRECISION DDOT, DNRM2
273 INTEGER IDAMAX
274 LOGICAL LSAME
275 EXTERNAL idamax, lsame, ddot, dnrm2
276* ..
277* .. External Subroutines ..
278 EXTERNAL daxpy, dcopy, dlascl, dlassq, drotm,
279 $ dswap,
280 $ xerbla
281* ..
282* .. Executable Statements ..
283*
284* Test the input parameters.
285*
286 applv = lsame( jobv, 'A' )
287 rsvec = lsame( jobv, 'V' )
288 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
289 info = -1
290 ELSE IF( m.LT.0 ) THEN
291 info = -2
292 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
293 info = -3
294 ELSE IF( n1.LT.0 ) THEN
295 info = -4
296 ELSE IF( lda.LT.m ) THEN
297 info = -6
298 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
299 info = -9
300 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
301 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
302 info = -11
303 ELSE IF( tol.LE.eps ) THEN
304 info = -14
305 ELSE IF( nsweep.LT.0 ) THEN
306 info = -15
307 ELSE IF( lwork.LT.m ) THEN
308 info = -17
309 ELSE
310 info = 0
311 END IF
312*
313* #:(
314 IF( info.NE.0 ) THEN
315 CALL xerbla( 'DGSVJ1', -info )
316 RETURN
317 END IF
318*
319 IF( rsvec ) THEN
320 mvl = n
321 ELSE IF( applv ) THEN
322 mvl = mv
323 END IF
324 rsvec = rsvec .OR. applv
325
326 rooteps = dsqrt( eps )
327 rootsfmin = dsqrt( sfmin )
328 small = sfmin / eps
329 big = one / sfmin
330 rootbig = one / rootsfmin
331 large = big / dsqrt( dble( m*n ) )
332 bigtheta = one / rooteps
333 roottol = dsqrt( tol )
334*
335* .. Initialize the right singular vector matrix ..
336*
337* RSVEC = LSAME( JOBV, 'Y' )
338*
339 emptsw = n1*( n-n1 )
340 notrot = 0
341 fastr( 1 ) = zero
342*
343* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
344*
345 kbl = min( 8, n )
346 nblr = n1 / kbl
347 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
348
349* .. the tiling is nblr-by-nblc [tiles]
350
351 nblc = ( n-n1 ) / kbl
352 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
353 blskip = ( kbl**2 ) + 1
354*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
355
356 rowskip = min( 5, kbl )
357*[TP] ROWSKIP is a tuning parameter.
358 swband = 0
359*[TP] SWBAND is a tuning parameter. It is meaningful and effective
360* if SGESVJ is used as a computational routine in the preconditioned
361* Jacobi SVD algorithm SGESVJ.
362*
363*
364* | * * * [x] [x] [x]|
365* | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
366* | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
367* |[x] [x] [x] * * * |
368* |[x] [x] [x] * * * |
369* |[x] [x] [x] * * * |
370*
371*
372 DO 1993 i = 1, nsweep
373* .. go go go ...
374*
375 mxaapq = zero
376 mxsinj = zero
377 iswrot = 0
378*
379 notrot = 0
380 pskipped = 0
381*
382 DO 2000 ibr = 1, nblr
383
384 igl = ( ibr-1 )*kbl + 1
385*
386*
387*........................................................
388* ... go to the off diagonal blocks
389
390 igl = ( ibr-1 )*kbl + 1
391
392 DO 2010 jbc = 1, nblc
393
394 jgl = n1 + ( jbc-1 )*kbl + 1
395
396* doing the block at ( ibr, jbc )
397
398 ijblsk = 0
399 DO 2100 p = igl, min( igl+kbl-1, n1 )
400
401 aapp = sva( p )
402
403 IF( aapp.GT.zero ) THEN
404
405 pskipped = 0
406
407 DO 2200 q = jgl, min( jgl+kbl-1, n )
408*
409 aaqq = sva( q )
410
411 IF( aaqq.GT.zero ) THEN
412 aapp0 = aapp
413*
414* .. M x 2 Jacobi SVD ..
415*
416* .. Safe Gram matrix computation ..
417*
418 IF( aaqq.GE.one ) THEN
419 IF( aapp.GE.aaqq ) THEN
420 rotok = ( small*aapp ).LE.aaqq
421 ELSE
422 rotok = ( small*aaqq ).LE.aapp
423 END IF
424 IF( aapp.LT.( big / aaqq ) ) THEN
425 aapq = ( ddot( m, a( 1, p ), 1,
426 $ a( 1,
427 $ q ), 1 )*d( p )*d( q ) / aaqq )
428 $ / aapp
429 ELSE
430 CALL dcopy( m, a( 1, p ), 1, work,
431 $ 1 )
432 CALL dlascl( 'G', 0, 0, aapp,
433 $ d( p ),
434 $ m, 1, work, lda, ierr )
435 aapq = ddot( m, work, 1, a( 1, q ),
436 $ 1 )*d( q ) / aaqq
437 END IF
438 ELSE
439 IF( aapp.GE.aaqq ) THEN
440 rotok = aapp.LE.( aaqq / small )
441 ELSE
442 rotok = aaqq.LE.( aapp / small )
443 END IF
444 IF( aapp.GT.( small / aaqq ) ) THEN
445 aapq = ( ddot( m, a( 1, p ), 1,
446 $ a( 1,
447 $ q ), 1 )*d( p )*d( q ) / aaqq )
448 $ / aapp
449 ELSE
450 CALL dcopy( m, a( 1, q ), 1, work,
451 $ 1 )
452 CALL dlascl( 'G', 0, 0, aaqq,
453 $ d( q ),
454 $ m, 1, work, lda, ierr )
455 aapq = ddot( m, work, 1, a( 1, p ),
456 $ 1 )*d( p ) / aapp
457 END IF
458 END IF
459
460 mxaapq = max( mxaapq, dabs( aapq ) )
461
462* TO rotate or NOT to rotate, THAT is the question ...
463*
464 IF( dabs( aapq ).GT.tol ) THEN
465 notrot = 0
466* ROTATED = ROTATED + 1
467 pskipped = 0
468 iswrot = iswrot + 1
469*
470 IF( rotok ) THEN
471*
472 aqoap = aaqq / aapp
473 apoaq = aapp / aaqq
474 theta = -half*dabs(aqoap-apoaq) / aapq
475 IF( aaqq.GT.aapp0 )theta = -theta
476
477 IF( dabs( theta ).GT.bigtheta ) THEN
478 t = half / theta
479 fastr( 3 ) = t*d( p ) / d( q )
480 fastr( 4 ) = -t*d( q ) / d( p )
481 CALL drotm( m, a( 1, p ), 1,
482 $ a( 1, q ), 1, fastr )
483 IF( rsvec )CALL drotm( mvl,
484 $ v( 1, p ), 1,
485 $ v( 1, q ), 1,
486 $ fastr )
487 sva( q ) = aaqq*dsqrt( max( zero,
488 $ one+t*apoaq*aapq ) )
489 aapp = aapp*dsqrt( max( zero,
490 $ one-t*aqoap*aapq ) )
491 mxsinj = max( mxsinj, dabs( t ) )
492 ELSE
493*
494* .. choose correct signum for THETA and rotate
495*
496 thsign = -dsign( one, aapq )
497 IF( aaqq.GT.aapp0 )thsign = -thsign
498 t = one / ( theta+thsign*
499 $ dsqrt( one+theta*theta ) )
500 cs = dsqrt( one / ( one+t*t ) )
501 sn = t*cs
502 mxsinj = max( mxsinj, dabs( sn ) )
503 sva( q ) = aaqq*dsqrt( max( zero,
504 $ one+t*apoaq*aapq ) )
505 aapp = aapp*dsqrt( max( zero,
506 $ one-t*aqoap*aapq ) )
507
508 apoaq = d( p ) / d( q )
509 aqoap = d( q ) / d( p )
510 IF( d( p ).GE.one ) THEN
511*
512 IF( d( q ).GE.one ) THEN
513 fastr( 3 ) = t*apoaq
514 fastr( 4 ) = -t*aqoap
515 d( p ) = d( p )*cs
516 d( q ) = d( q )*cs
517 CALL drotm( m, a( 1, p ),
518 $ 1,
519 $ a( 1, q ), 1,
520 $ fastr )
521 IF( rsvec )CALL drotm( mvl,
522 $ v( 1, p ), 1, v( 1, q ),
523 $ 1, fastr )
524 ELSE
525 CALL daxpy( m, -t*aqoap,
526 $ a( 1, q ), 1,
527 $ a( 1, p ), 1 )
528 CALL daxpy( m, cs*sn*apoaq,
529 $ a( 1, p ), 1,
530 $ a( 1, q ), 1 )
531 IF( rsvec ) THEN
532 CALL daxpy( mvl,
533 $ -t*aqoap,
534 $ v( 1, q ), 1,
535 $ v( 1, p ), 1 )
536 CALL daxpy( mvl,
537 $ cs*sn*apoaq,
538 $ v( 1, p ), 1,
539 $ v( 1, q ), 1 )
540 END IF
541 d( p ) = d( p )*cs
542 d( q ) = d( q ) / cs
543 END IF
544 ELSE
545 IF( d( q ).GE.one ) THEN
546 CALL daxpy( m, t*apoaq,
547 $ a( 1, p ), 1,
548 $ a( 1, q ), 1 )
549 CALL daxpy( m,
550 $ -cs*sn*aqoap,
551 $ a( 1, q ), 1,
552 $ a( 1, p ), 1 )
553 IF( rsvec ) THEN
554 CALL daxpy( mvl,
555 $ t*apoaq,
556 $ v( 1, p ), 1,
557 $ v( 1, q ), 1 )
558 CALL daxpy( mvl,
559 $ -cs*sn*aqoap,
560 $ v( 1, q ), 1,
561 $ v( 1, p ), 1 )
562 END IF
563 d( p ) = d( p ) / cs
564 d( q ) = d( q )*cs
565 ELSE
566 IF( d( p ).GE.d( q ) ) THEN
567 CALL daxpy( m, -t*aqoap,
568 $ a( 1, q ), 1,
569 $ a( 1, p ), 1 )
570 CALL daxpy( m,
571 $ cs*sn*apoaq,
572 $ a( 1, p ), 1,
573 $ a( 1, q ), 1 )
574 d( p ) = d( p )*cs
575 d( q ) = d( q ) / cs
576 IF( rsvec ) THEN
577 CALL daxpy( mvl,
578 $ -t*aqoap,
579 $ v( 1, q ), 1,
580 $ v( 1, p ), 1 )
581 CALL daxpy( mvl,
582 $ cs*sn*apoaq,
583 $ v( 1, p ), 1,
584 $ v( 1, q ), 1 )
585 END IF
586 ELSE
587 CALL daxpy( m, t*apoaq,
588 $ a( 1, p ), 1,
589 $ a( 1, q ), 1 )
590 CALL daxpy( m,
591 $ -cs*sn*aqoap,
592 $ a( 1, q ), 1,
593 $ a( 1, p ), 1 )
594 d( p ) = d( p ) / cs
595 d( q ) = d( q )*cs
596 IF( rsvec ) THEN
597 CALL daxpy( mvl,
598 $ t*apoaq, v( 1, p ),
599 $ 1, v( 1, q ), 1 )
600 CALL daxpy( mvl,
601 $ -cs*sn*aqoap,
602 $ v( 1, q ), 1,
603 $ v( 1, p ), 1 )
604 END IF
605 END IF
606 END IF
607 END IF
608 END IF
609
610 ELSE
611 IF( aapp.GT.aaqq ) THEN
612 CALL dcopy( m, a( 1, p ), 1,
613 $ work,
614 $ 1 )
615 CALL dlascl( 'G', 0, 0, aapp,
616 $ one,
617 $ m, 1, work, lda, ierr )
618 CALL dlascl( 'G', 0, 0, aaqq,
619 $ one,
620 $ m, 1, a( 1, q ), lda,
621 $ ierr )
622 temp1 = -aapq*d( p ) / d( q )
623 CALL daxpy( m, temp1, work, 1,
624 $ a( 1, q ), 1 )
625 CALL dlascl( 'G', 0, 0, one,
626 $ aaqq,
627 $ m, 1, a( 1, q ), lda,
628 $ ierr )
629 sva( q ) = aaqq*dsqrt( max( zero,
630 $ one-aapq*aapq ) )
631 mxsinj = max( mxsinj, sfmin )
632 ELSE
633 CALL dcopy( m, a( 1, q ), 1,
634 $ work,
635 $ 1 )
636 CALL dlascl( 'G', 0, 0, aaqq,
637 $ one,
638 $ m, 1, work, lda, ierr )
639 CALL dlascl( 'G', 0, 0, aapp,
640 $ one,
641 $ m, 1, a( 1, p ), lda,
642 $ ierr )
643 temp1 = -aapq*d( q ) / d( p )
644 CALL daxpy( m, temp1, work, 1,
645 $ a( 1, p ), 1 )
646 CALL dlascl( 'G', 0, 0, one,
647 $ aapp,
648 $ m, 1, a( 1, p ), lda,
649 $ ierr )
650 sva( p ) = aapp*dsqrt( max( zero,
651 $ one-aapq*aapq ) )
652 mxsinj = max( mxsinj, sfmin )
653 END IF
654 END IF
655* END IF ROTOK THEN ... ELSE
656*
657* In the case of cancellation in updating SVA(q)
658* .. recompute SVA(q)
659 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
660 $ THEN
661 IF( ( aaqq.LT.rootbig ) .AND.
662 $ ( aaqq.GT.rootsfmin ) ) THEN
663 sva( q ) = dnrm2( m, a( 1, q ),
664 $ 1 )*
665 $ d( q )
666 ELSE
667 t = zero
668 aaqq = one
669 CALL dlassq( m, a( 1, q ), 1, t,
670 $ aaqq )
671 sva( q ) = t*dsqrt( aaqq )*d( q )
672 END IF
673 END IF
674 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
675 IF( ( aapp.LT.rootbig ) .AND.
676 $ ( aapp.GT.rootsfmin ) ) THEN
677 aapp = dnrm2( m, a( 1, p ), 1 )*
678 $ d( p )
679 ELSE
680 t = zero
681 aapp = one
682 CALL dlassq( m, a( 1, p ), 1, t,
683 $ aapp )
684 aapp = t*dsqrt( aapp )*d( p )
685 END IF
686 sva( p ) = aapp
687 END IF
688* end of OK rotation
689 ELSE
690 notrot = notrot + 1
691* SKIPPED = SKIPPED + 1
692 pskipped = pskipped + 1
693 ijblsk = ijblsk + 1
694 END IF
695 ELSE
696 notrot = notrot + 1
697 pskipped = pskipped + 1
698 ijblsk = ijblsk + 1
699 END IF
700
701* IF ( NOTROT .GE. EMPTSW ) GO TO 2011
702 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
703 $ THEN
704 sva( p ) = aapp
705 notrot = 0
706 GO TO 2011
707 END IF
708 IF( ( i.LE.swband ) .AND.
709 $ ( pskipped.GT.rowskip ) ) THEN
710 aapp = -aapp
711 notrot = 0
712 GO TO 2203
713 END IF
714
715*
716 2200 CONTINUE
717* end of the q-loop
718 2203 CONTINUE
719
720 sva( p ) = aapp
721*
722 ELSE
723 IF( aapp.EQ.zero )notrot = notrot +
724 $ min( jgl+kbl-1, n ) - jgl + 1
725 IF( aapp.LT.zero )notrot = 0
726*** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
727 END IF
728
729 2100 CONTINUE
730* end of the p-loop
731 2010 CONTINUE
732* end of the jbc-loop
733 2011 CONTINUE
734*2011 bailed out of the jbc-loop
735 DO 2012 p = igl, min( igl+kbl-1, n )
736 sva( p ) = dabs( sva( p ) )
737 2012 CONTINUE
738*** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
739 2000 CONTINUE
740*2000 :: end of the ibr-loop
741*
742* .. update SVA(N)
743 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
744 $ THEN
745 sva( n ) = dnrm2( m, a( 1, n ), 1 )*d( n )
746 ELSE
747 t = zero
748 aapp = one
749 CALL dlassq( m, a( 1, n ), 1, t, aapp )
750 sva( n ) = t*dsqrt( aapp )*d( n )
751 END IF
752*
753* Additional steering devices
754*
755 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
756 $ ( iswrot.LE.n ) ) )swband = i
757
758 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
759 $ ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
760 GO TO 1994
761 END IF
762
763*
764 IF( notrot.GE.emptsw )GO TO 1994
765
766 1993 CONTINUE
767* end i=1:NSWEEP loop
768* #:) Reaching this point means that the procedure has completed the given
769* number of sweeps.
770 info = nsweep - 1
771 GO TO 1995
772 1994 CONTINUE
773* #:) Reaching this point means that during the i-th sweep all pivots were
774* below the given threshold, causing early exit.
775
776 info = 0
777* #:) INFO = 0 confirms successful iterations.
778 1995 CONTINUE
779*
780* Sort the vector D
781*
782 DO 5991 p = 1, n - 1
783 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
784 IF( p.NE.q ) THEN
785 temp1 = sva( p )
786 sva( p ) = sva( q )
787 sva( q ) = temp1
788 temp1 = d( p )
789 d( p ) = d( q )
790 d( q ) = temp1
791 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
792 IF( rsvec )CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
793 END IF
794 5991 CONTINUE
795*
796 RETURN
797* ..
798* .. END OF DGSVJ1
799* ..
800 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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
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:234
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:142
dlassq
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:122
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