LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgsvj0.f
Go to the documentation of this file.
1*> \brief \b SGSVJ0 pre-processor for the routine sgesvj.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SGSVJ0 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgsvj0.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgsvj0.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgsvj0.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
20* SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
24* REAL EPS, SFMIN, TOL
25* CHARACTER*1 JOBV
26* ..
27* .. Array Arguments ..
28* REAL A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
29* $ WORK( LWORK )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SGSVJ0 is called from SGESVJ as a pre-processor and that is its main
39*> purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
40*> it does not check convergence (stopping criterion). Few tuning
41*> parameters (marked by [TP]) are available for the implementer.
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] JOBV
48*> \verbatim
49*> JOBV is CHARACTER*1
50*> Specifies whether the output from this procedure is used
51*> to compute the matrix V:
52*> = 'V': the product of the Jacobi rotations is accumulated
53*> by postmultiplying the N-by-N array V.
54*> (See the description of V.)
55*> = 'A': the product of the Jacobi rotations is accumulated
56*> by postmultiplying the MV-by-N array V.
57*> (See the descriptions of MV and V.)
58*> = 'N': the Jacobi rotations are not accumulated.
59*> \endverbatim
60*>
61*> \param[in] M
62*> \verbatim
63*> M is INTEGER
64*> The number of rows of the input matrix A. M >= 0.
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*> N is INTEGER
70*> The number of columns of the input matrix A.
71*> M >= N >= 0.
72*> \endverbatim
73*>
74*> \param[in,out] A
75*> \verbatim
76*> A is REAL array, dimension (LDA,N)
77*> On entry, M-by-N matrix A, such that A*diag(D) represents
78*> the input matrix.
79*> On exit,
80*> A_onexit * D_onexit represents the input matrix A*diag(D)
81*> post-multiplied by a sequence of Jacobi rotations, where the
82*> rotation threshold and the total number of sweeps are given in
83*> TOL and NSWEEP, respectively.
84*> (See the descriptions of D, TOL and NSWEEP.)
85*> \endverbatim
86*>
87*> \param[in] LDA
88*> \verbatim
89*> LDA is INTEGER
90*> The leading dimension of the array A. LDA >= max(1,M).
91*> \endverbatim
92*>
93*> \param[in,out] D
94*> \verbatim
95*> D is REAL array, dimension (N)
96*> The array D accumulates the scaling factors from the fast scaled
97*> Jacobi rotations.
98*> On entry, A*diag(D) represents the input matrix.
99*> On exit, A_onexit*diag(D_onexit) represents the input matrix
100*> post-multiplied by a sequence of Jacobi rotations, where the
101*> rotation threshold and the total number of sweeps are given in
102*> TOL and NSWEEP, respectively.
103*> (See the descriptions of A, TOL and NSWEEP.)
104*> \endverbatim
105*>
106*> \param[in,out] SVA
107*> \verbatim
108*> SVA is REAL array, dimension (N)
109*> On entry, SVA contains the Euclidean norms of the columns of
110*> the matrix A*diag(D).
111*> On exit, SVA contains the Euclidean norms of the columns of
112*> the matrix onexit*diag(D_onexit).
113*> \endverbatim
114*>
115*> \param[in] MV
116*> \verbatim
117*> MV is INTEGER
118*> If JOBV = 'A', then MV rows of V are post-multiplied by a
119*> sequence of Jacobi rotations.
120*> If JOBV = 'N', then MV is not referenced.
121*> \endverbatim
122*>
123*> \param[in,out] V
124*> \verbatim
125*> V is REAL array, dimension (LDV,N)
126*> If JOBV = 'V' then N rows of V are post-multiplied by a
127*> sequence of Jacobi rotations.
128*> If JOBV = 'A' then MV rows of V are post-multiplied by a
129*> sequence of Jacobi rotations.
130*> If JOBV = 'N', then V is not referenced.
131*> \endverbatim
132*>
133*> \param[in] LDV
134*> \verbatim
135*> LDV is INTEGER
136*> The leading dimension of the array V, LDV >= 1.
137*> If JOBV = 'V', LDV >= N.
138*> If JOBV = 'A', LDV >= MV.
139*> \endverbatim
140*>
141*> \param[in] EPS
142*> \verbatim
143*> EPS is REAL
144*> EPS = SLAMCH('Epsilon')
145*> \endverbatim
146*>
147*> \param[in] SFMIN
148*> \verbatim
149*> SFMIN is REAL
150*> SFMIN = SLAMCH('Safe Minimum')
151*> \endverbatim
152*>
153*> \param[in] TOL
154*> \verbatim
155*> TOL is REAL
156*> TOL is the threshold for Jacobi rotations. For a pair
157*> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
158*> applied only if ABS(COS(angle(A(:,p),A(:,q)))) > TOL.
159*> \endverbatim
160*>
161*> \param[in] NSWEEP
162*> \verbatim
163*> NSWEEP is INTEGER
164*> NSWEEP is the number of sweeps of Jacobi rotations to be
165*> performed.
166*> \endverbatim
167*>
168*> \param[out] WORK
169*> \verbatim
170*> WORK is REAL array, dimension (LWORK)
171*> \endverbatim
172*>
173*> \param[in] LWORK
174*> \verbatim
175*> LWORK is INTEGER
176*> LWORK is the dimension of WORK. LWORK >= M.
177*> \endverbatim
178*>
179*> \param[out] INFO
180*> \verbatim
181*> INFO is INTEGER
182*> = 0: successful exit.
183*> < 0: if INFO = -i, then the i-th argument had an illegal value
184*> \endverbatim
185*
186* Authors:
187* ========
188*
189*> \author Univ. of Tennessee
190*> \author Univ. of California Berkeley
191*> \author Univ. of Colorado Denver
192*> \author NAG Ltd.
193*
194*> \ingroup gsvj0
195*
196*> \par Further Details:
197* =====================
198*>
199*> SGSVJ0 is used just to enable SGESVJ to call a simplified version of
200*> itself to work on a submatrix of the original matrix.
201*>
202*> \par Contributors:
203* ==================
204*>
205*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
206*>
207*> \par Bugs, Examples and Comments:
208* =================================
209*>
210*> Please report all bugs and send interesting test examples and comments to
211*> drmac@math.hr. Thank you.
212*
213* =====================================================================
214 SUBROUTINE sgsvj0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
215 $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
216*
217* -- LAPACK computational routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
223 REAL EPS, SFMIN, TOL
224 CHARACTER*1 JOBV
225* ..
226* .. Array Arguments ..
227 REAL A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
228 $ work( lwork )
229* ..
230*
231* =====================================================================
232*
233* .. Local Parameters ..
234 REAL ZERO, HALF, ONE
235 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
236* ..
237* .. Local Scalars ..
238 REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
239 $ bigtheta, cs, mxaapq, mxsinj, rootbig, rooteps,
240 $ rootsfmin, roottol, small, sn, t, temp1, theta,
241 $ thsign
242 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
243 $ iswrot, jbc, jgl, kbl, lkahead, mvl, nbl,
244 $ notrot, p, pskipped, q, rowskip, swband
245 LOGICAL APPLV, ROTOK, RSVEC
246* ..
247* .. Local Arrays ..
248 REAL FASTR( 5 )
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC abs, max, float, min, sign, sqrt
252* ..
253* .. External Functions ..
254 REAL SDOT, SNRM2
255 INTEGER ISAMAX
256 LOGICAL LSAME
257 EXTERNAL isamax, lsame, sdot, snrm2
258* ..
259* .. External Subroutines ..
260 EXTERNAL saxpy, scopy, slascl, slassq, srotm,
261 $ sswap,
262 $ xerbla
263* ..
264* .. Executable Statements ..
265*
266* Test the input parameters.
267*
268 applv = lsame( jobv, 'A' )
269 rsvec = lsame( jobv, 'V' )
270 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
271 info = -1
272 ELSE IF( m.LT.0 ) THEN
273 info = -2
274 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
275 info = -3
276 ELSE IF( lda.LT.m ) THEN
277 info = -5
278 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
279 info = -8
280 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
281 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
282 info = -10
283 ELSE IF( tol.LE.eps ) THEN
284 info = -13
285 ELSE IF( nsweep.LT.0 ) THEN
286 info = -14
287 ELSE IF( lwork.LT.m ) THEN
288 info = -16
289 ELSE
290 info = 0
291 END IF
292*
293* #:(
294 IF( info.NE.0 ) THEN
295 CALL xerbla( 'SGSVJ0', -info )
296 RETURN
297 END IF
298*
299 IF( rsvec ) THEN
300 mvl = n
301 ELSE IF( applv ) THEN
302 mvl = mv
303 END IF
304 rsvec = rsvec .OR. applv
305
306 rooteps = sqrt( eps )
307 rootsfmin = sqrt( sfmin )
308 small = sfmin / eps
309 big = one / sfmin
310 rootbig = one / rootsfmin
311 bigtheta = one / rooteps
312 roottol = sqrt( tol )
313*
314* .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
315*
316 emptsw = ( n*( n-1 ) ) / 2
317 notrot = 0
318 fastr( 1 ) = zero
319*
320* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
321*
322
323 swband = 0
324*[TP] SWBAND is a tuning parameter. It is meaningful and effective
325* if SGESVJ is used as a computational routine in the preconditioned
326* Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
327* ......
328
329 kbl = min( 8, n )
330*[TP] KBL is a tuning parameter that defines the tile size in the
331* tiling of the p-q loops of pivot pairs. In general, an optimal
332* value of KBL depends on the matrix dimensions and on the
333* parameters of the computer's memory.
334*
335 nbl = n / kbl
336 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
337
338 blskip = ( kbl**2 ) + 1
339*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
340
341 rowskip = min( 5, kbl )
342*[TP] ROWSKIP is a tuning parameter.
343
344 lkahead = 1
345*[TP] LKAHEAD is a tuning parameter.
346 swband = 0
347 pskipped = 0
348*
349 DO 1993 i = 1, nsweep
350* .. go go go ...
351*
352 mxaapq = zero
353 mxsinj = zero
354 iswrot = 0
355*
356 notrot = 0
357 pskipped = 0
358*
359 DO 2000 ibr = 1, nbl
360
361 igl = ( ibr-1 )*kbl + 1
362*
363 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
364*
365 igl = igl + ir1*kbl
366*
367 DO 2001 p = igl, min( igl+kbl-1, n-1 )
368
369* .. de Rijk's pivoting
370 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
371 IF( p.NE.q ) THEN
372 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
373 IF( rsvec )CALL sswap( mvl, v( 1, p ), 1,
374 $ v( 1, q ), 1 )
375 temp1 = sva( p )
376 sva( p ) = sva( q )
377 sva( q ) = temp1
378 temp1 = d( p )
379 d( p ) = d( q )
380 d( q ) = temp1
381 END IF
382*
383 IF( ir1.EQ.0 ) THEN
384*
385* Column norms are periodically updated by explicit
386* norm computation.
387* Caveat:
388* Some BLAS implementations compute SNRM2(M,A(1,p),1)
389* as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may result in
390* overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and
391* underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
392* Hence, SNRM2 cannot be trusted, not even in the case when
393* the true norm is far from the under(over)flow boundaries.
394* If properly implemented SNRM2 is available, the IF-THEN-ELSE
395* below should read "AAPP = SNRM2( M, A(1,p), 1 ) * D(p)".
396*
397 IF( ( sva( p ).LT.rootbig ) .AND.
398 $ ( sva( p ).GT.rootsfmin ) ) THEN
399 sva( p ) = snrm2( m, a( 1, p ), 1 )*d( p )
400 ELSE
401 temp1 = zero
402 aapp = one
403 CALL slassq( m, a( 1, p ), 1, temp1, aapp )
404 sva( p ) = temp1*sqrt( aapp )*d( p )
405 END IF
406 aapp = sva( p )
407 ELSE
408 aapp = sva( p )
409 END IF
410
411*
412 IF( aapp.GT.zero ) THEN
413*
414 pskipped = 0
415*
416 DO 2002 q = p + 1, min( igl+kbl-1, n )
417*
418 aaqq = sva( q )
419
420 IF( aaqq.GT.zero ) THEN
421*
422 aapp0 = aapp
423 IF( aaqq.GE.one ) THEN
424 rotok = ( small*aapp ).LE.aaqq
425 IF( aapp.LT.( big / aaqq ) ) THEN
426 aapq = ( sdot( m, a( 1, p ), 1,
427 $ a( 1,
428 $ q ), 1 )*d( p )*d( q ) / aaqq )
429 $ / aapp
430 ELSE
431 CALL scopy( m, a( 1, p ), 1, work,
432 $ 1 )
433 CALL slascl( 'G', 0, 0, aapp,
434 $ d( p ),
435 $ m, 1, work, lda, ierr )
436 aapq = sdot( m, work, 1, a( 1, q ),
437 $ 1 )*d( q ) / aaqq
438 END IF
439 ELSE
440 rotok = aapp.LE.( aaqq / small )
441 IF( aapp.GT.( small / aaqq ) ) THEN
442 aapq = ( sdot( m, a( 1, p ), 1,
443 $ a( 1,
444 $ q ), 1 )*d( p )*d( q ) / aaqq )
445 $ / aapp
446 ELSE
447 CALL scopy( m, a( 1, q ), 1, work,
448 $ 1 )
449 CALL slascl( 'G', 0, 0, aaqq,
450 $ d( q ),
451 $ m, 1, work, lda, ierr )
452 aapq = sdot( m, work, 1, a( 1, p ),
453 $ 1 )*d( p ) / aapp
454 END IF
455 END IF
456*
457 mxaapq = max( mxaapq, abs( aapq ) )
458*
459* TO rotate or NOT to rotate, THAT is the question ...
460*
461 IF( abs( aapq ).GT.tol ) THEN
462*
463* .. rotate
464* ROTATED = ROTATED + ONE
465*
466 IF( ir1.EQ.0 ) THEN
467 notrot = 0
468 pskipped = 0
469 iswrot = iswrot + 1
470 END IF
471*
472 IF( rotok ) THEN
473*
474 aqoap = aaqq / aapp
475 apoaq = aapp / aaqq
476 theta = -half*abs( aqoap-apoaq ) / aapq
477*
478 IF( abs( theta ).GT.bigtheta ) THEN
479*
480 t = half / theta
481 fastr( 3 ) = t*d( p ) / d( q )
482 fastr( 4 ) = -t*d( q ) / d( p )
483 CALL srotm( m, a( 1, p ), 1,
484 $ a( 1, q ), 1, fastr )
485 IF( rsvec )CALL srotm( mvl,
486 $ v( 1, p ), 1,
487 $ v( 1, q ), 1,
488 $ fastr )
489 sva( q ) = aaqq*sqrt( max( zero,
490 $ one+t*apoaq*aapq ) )
491 aapp = aapp*sqrt( max( zero,
492 $ one-t*aqoap*aapq ) )
493 mxsinj = max( mxsinj, abs( t ) )
494*
495 ELSE
496*
497* .. choose correct signum for THETA and rotate
498*
499 thsign = -sign( one, aapq )
500 t = one / ( theta+thsign*
501 $ sqrt( one+theta*theta ) )
502 cs = sqrt( one / ( one+t*t ) )
503 sn = t*cs
504*
505 mxsinj = max( mxsinj, abs( sn ) )
506 sva( q ) = aaqq*sqrt( max( zero,
507 $ one+t*apoaq*aapq ) )
508 aapp = aapp*sqrt( max( zero,
509 $ one-t*aqoap*aapq ) )
510*
511 apoaq = d( p ) / d( q )
512 aqoap = d( q ) / d( p )
513 IF( d( p ).GE.one ) THEN
514 IF( d( q ).GE.one ) THEN
515 fastr( 3 ) = t*apoaq
516 fastr( 4 ) = -t*aqoap
517 d( p ) = d( p )*cs
518 d( q ) = d( q )*cs
519 CALL srotm( m, a( 1, p ),
520 $ 1,
521 $ a( 1, q ), 1,
522 $ fastr )
523 IF( rsvec )CALL srotm( mvl,
524 $ v( 1, p ), 1, v( 1, q ),
525 $ 1, fastr )
526 ELSE
527 CALL saxpy( m, -t*aqoap,
528 $ a( 1, q ), 1,
529 $ a( 1, p ), 1 )
530 CALL saxpy( m, cs*sn*apoaq,
531 $ a( 1, p ), 1,
532 $ a( 1, q ), 1 )
533 d( p ) = d( p )*cs
534 d( q ) = d( q ) / cs
535 IF( rsvec ) THEN
536 CALL saxpy( mvl,
537 $ -t*aqoap,
538 $ v( 1, q ), 1,
539 $ v( 1, p ), 1 )
540 CALL saxpy( mvl,
541 $ cs*sn*apoaq,
542 $ v( 1, p ), 1,
543 $ v( 1, q ), 1 )
544 END IF
545 END IF
546 ELSE
547 IF( d( q ).GE.one ) THEN
548 CALL saxpy( m, t*apoaq,
549 $ a( 1, p ), 1,
550 $ a( 1, q ), 1 )
551 CALL saxpy( m,
552 $ -cs*sn*aqoap,
553 $ a( 1, q ), 1,
554 $ a( 1, p ), 1 )
555 d( p ) = d( p ) / cs
556 d( q ) = d( q )*cs
557 IF( rsvec ) THEN
558 CALL saxpy( mvl,
559 $ t*apoaq,
560 $ v( 1, p ), 1,
561 $ v( 1, q ), 1 )
562 CALL saxpy( mvl,
563 $ -cs*sn*aqoap,
564 $ v( 1, q ), 1,
565 $ v( 1, p ), 1 )
566 END IF
567 ELSE
568 IF( d( p ).GE.d( q ) ) THEN
569 CALL saxpy( m, -t*aqoap,
570 $ a( 1, q ), 1,
571 $ a( 1, p ), 1 )
572 CALL saxpy( m,
573 $ cs*sn*apoaq,
574 $ a( 1, p ), 1,
575 $ a( 1, q ), 1 )
576 d( p ) = d( p )*cs
577 d( q ) = d( q ) / cs
578 IF( rsvec ) THEN
579 CALL saxpy( mvl,
580 $ -t*aqoap,
581 $ v( 1, q ), 1,
582 $ v( 1, p ), 1 )
583 CALL saxpy( mvl,
584 $ cs*sn*apoaq,
585 $ v( 1, p ), 1,
586 $ v( 1, q ), 1 )
587 END IF
588 ELSE
589 CALL saxpy( m, t*apoaq,
590 $ a( 1, p ), 1,
591 $ a( 1, q ), 1 )
592 CALL saxpy( m,
593 $ -cs*sn*aqoap,
594 $ a( 1, q ), 1,
595 $ a( 1, p ), 1 )
596 d( p ) = d( p ) / cs
597 d( q ) = d( q )*cs
598 IF( rsvec ) THEN
599 CALL saxpy( mvl,
600 $ t*apoaq, v( 1, p ),
601 $ 1, v( 1, q ), 1 )
602 CALL saxpy( mvl,
603 $ -cs*sn*aqoap,
604 $ v( 1, q ), 1,
605 $ v( 1, p ), 1 )
606 END IF
607 END IF
608 END IF
609 END IF
610 END IF
611*
612 ELSE
613* .. have to use modified Gram-Schmidt like transformation
614 CALL scopy( m, a( 1, p ), 1, work,
615 $ 1 )
616 CALL slascl( 'G', 0, 0, aapp, one,
617 $ m,
618 $ 1, work, lda, ierr )
619 CALL slascl( 'G', 0, 0, aaqq, one,
620 $ m,
621 $ 1, a( 1, q ), lda, ierr )
622 temp1 = -aapq*d( p ) / d( q )
623 CALL saxpy( m, temp1, work, 1,
624 $ a( 1, q ), 1 )
625 CALL slascl( 'G', 0, 0, one, aaqq,
626 $ m,
627 $ 1, a( 1, q ), lda, ierr )
628 sva( q ) = aaqq*sqrt( max( zero,
629 $ one-aapq*aapq ) )
630 mxsinj = max( mxsinj, sfmin )
631 END IF
632* END IF ROTOK THEN ... ELSE
633*
634* In the case of cancellation in updating SVA(q), SVA(p)
635* recompute SVA(q), SVA(p).
636 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
637 $ THEN
638 IF( ( aaqq.LT.rootbig ) .AND.
639 $ ( aaqq.GT.rootsfmin ) ) THEN
640 sva( q ) = snrm2( m, a( 1, q ),
641 $ 1 )*
642 $ d( q )
643 ELSE
644 t = zero
645 aaqq = one
646 CALL slassq( m, a( 1, q ), 1, t,
647 $ aaqq )
648 sva( q ) = t*sqrt( aaqq )*d( q )
649 END IF
650 END IF
651 IF( ( aapp / aapp0 ).LE.rooteps ) THEN
652 IF( ( aapp.LT.rootbig ) .AND.
653 $ ( aapp.GT.rootsfmin ) ) THEN
654 aapp = snrm2( m, a( 1, p ), 1 )*
655 $ d( p )
656 ELSE
657 t = zero
658 aapp = one
659 CALL slassq( m, a( 1, p ), 1, t,
660 $ aapp )
661 aapp = t*sqrt( aapp )*d( p )
662 END IF
663 sva( p ) = aapp
664 END IF
665*
666 ELSE
667* A(:,p) and A(:,q) already numerically orthogonal
668 IF( ir1.EQ.0 )notrot = notrot + 1
669 pskipped = pskipped + 1
670 END IF
671 ELSE
672* A(:,q) is zero column
673 IF( ir1.EQ.0 )notrot = notrot + 1
674 pskipped = pskipped + 1
675 END IF
676*
677 IF( ( i.LE.swband ) .AND.
678 $ ( pskipped.GT.rowskip ) ) THEN
679 IF( ir1.EQ.0 )aapp = -aapp
680 notrot = 0
681 GO TO 2103
682 END IF
683*
684 2002 CONTINUE
685* END q-LOOP
686*
687 2103 CONTINUE
688* bailed out of q-loop
689
690 sva( p ) = aapp
691
692 ELSE
693 sva( p ) = aapp
694 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
695 $ notrot = notrot + min( igl+kbl-1, n ) - p
696 END IF
697*
698 2001 CONTINUE
699* end of the p-loop
700* end of doing the block ( ibr, ibr )
701 1002 CONTINUE
702* end of ir1-loop
703*
704*........................................................
705* ... go to the off diagonal blocks
706*
707 igl = ( ibr-1 )*kbl + 1
708*
709 DO 2010 jbc = ibr + 1, nbl
710*
711 jgl = ( jbc-1 )*kbl + 1
712*
713* doing the block at ( ibr, jbc )
714*
715 ijblsk = 0
716 DO 2100 p = igl, min( igl+kbl-1, n )
717*
718 aapp = sva( p )
719*
720 IF( aapp.GT.zero ) THEN
721*
722 pskipped = 0
723*
724 DO 2200 q = jgl, min( jgl+kbl-1, n )
725*
726 aaqq = sva( q )
727*
728 IF( aaqq.GT.zero ) THEN
729 aapp0 = aapp
730*
731* .. M x 2 Jacobi SVD ..
732*
733* .. Safe Gram matrix computation ..
734*
735 IF( aaqq.GE.one ) THEN
736 IF( aapp.GE.aaqq ) THEN
737 rotok = ( small*aapp ).LE.aaqq
738 ELSE
739 rotok = ( small*aaqq ).LE.aapp
740 END IF
741 IF( aapp.LT.( big / aaqq ) ) THEN
742 aapq = ( sdot( m, a( 1, p ), 1,
743 $ a( 1,
744 $ q ), 1 )*d( p )*d( q ) / aaqq )
745 $ / aapp
746 ELSE
747 CALL scopy( m, a( 1, p ), 1, work,
748 $ 1 )
749 CALL slascl( 'G', 0, 0, aapp,
750 $ d( p ),
751 $ m, 1, work, lda, ierr )
752 aapq = sdot( m, work, 1, a( 1, q ),
753 $ 1 )*d( q ) / aaqq
754 END IF
755 ELSE
756 IF( aapp.GE.aaqq ) THEN
757 rotok = aapp.LE.( aaqq / small )
758 ELSE
759 rotok = aaqq.LE.( aapp / small )
760 END IF
761 IF( aapp.GT.( small / aaqq ) ) THEN
762 aapq = ( sdot( m, a( 1, p ), 1,
763 $ a( 1,
764 $ q ), 1 )*d( p )*d( q ) / aaqq )
765 $ / aapp
766 ELSE
767 CALL scopy( m, a( 1, q ), 1, work,
768 $ 1 )
769 CALL slascl( 'G', 0, 0, aaqq,
770 $ d( q ),
771 $ m, 1, work, lda, ierr )
772 aapq = sdot( m, work, 1, a( 1, p ),
773 $ 1 )*d( p ) / aapp
774 END IF
775 END IF
776*
777 mxaapq = max( mxaapq, abs( aapq ) )
778*
779* TO rotate or NOT to rotate, THAT is the question ...
780*
781 IF( abs( aapq ).GT.tol ) THEN
782 notrot = 0
783* ROTATED = ROTATED + 1
784 pskipped = 0
785 iswrot = iswrot + 1
786*
787 IF( rotok ) THEN
788*
789 aqoap = aaqq / aapp
790 apoaq = aapp / aaqq
791 theta = -half*abs( aqoap-apoaq ) / aapq
792 IF( aaqq.GT.aapp0 )theta = -theta
793*
794 IF( abs( theta ).GT.bigtheta ) THEN
795 t = half / theta
796 fastr( 3 ) = t*d( p ) / d( q )
797 fastr( 4 ) = -t*d( q ) / d( p )
798 CALL srotm( m, a( 1, p ), 1,
799 $ a( 1, q ), 1, fastr )
800 IF( rsvec )CALL srotm( mvl,
801 $ v( 1, p ), 1,
802 $ v( 1, q ), 1,
803 $ fastr )
804 sva( q ) = aaqq*sqrt( max( zero,
805 $ one+t*apoaq*aapq ) )
806 aapp = aapp*sqrt( max( zero,
807 $ one-t*aqoap*aapq ) )
808 mxsinj = max( mxsinj, abs( t ) )
809 ELSE
810*
811* .. choose correct signum for THETA and rotate
812*
813 thsign = -sign( one, aapq )
814 IF( aaqq.GT.aapp0 )thsign = -thsign
815 t = one / ( theta+thsign*
816 $ sqrt( one+theta*theta ) )
817 cs = sqrt( one / ( one+t*t ) )
818 sn = t*cs
819 mxsinj = max( mxsinj, abs( sn ) )
820 sva( q ) = aaqq*sqrt( max( zero,
821 $ one+t*apoaq*aapq ) )
822 aapp = aapp*sqrt( max( zero,
823 $ one-t*aqoap*aapq ) )
824*
825 apoaq = d( p ) / d( q )
826 aqoap = d( q ) / d( p )
827 IF( d( p ).GE.one ) THEN
828*
829 IF( d( q ).GE.one ) THEN
830 fastr( 3 ) = t*apoaq
831 fastr( 4 ) = -t*aqoap
832 d( p ) = d( p )*cs
833 d( q ) = d( q )*cs
834 CALL srotm( m, a( 1, p ),
835 $ 1,
836 $ a( 1, q ), 1,
837 $ fastr )
838 IF( rsvec )CALL srotm( mvl,
839 $ v( 1, p ), 1, v( 1, q ),
840 $ 1, fastr )
841 ELSE
842 CALL saxpy( m, -t*aqoap,
843 $ a( 1, q ), 1,
844 $ a( 1, p ), 1 )
845 CALL saxpy( m, cs*sn*apoaq,
846 $ a( 1, p ), 1,
847 $ a( 1, q ), 1 )
848 IF( rsvec ) THEN
849 CALL saxpy( mvl,
850 $ -t*aqoap,
851 $ v( 1, q ), 1,
852 $ v( 1, p ), 1 )
853 CALL saxpy( mvl,
854 $ cs*sn*apoaq,
855 $ v( 1, p ), 1,
856 $ v( 1, q ), 1 )
857 END IF
858 d( p ) = d( p )*cs
859 d( q ) = d( q ) / cs
860 END IF
861 ELSE
862 IF( d( q ).GE.one ) THEN
863 CALL saxpy( m, t*apoaq,
864 $ a( 1, p ), 1,
865 $ a( 1, q ), 1 )
866 CALL saxpy( m,
867 $ -cs*sn*aqoap,
868 $ a( 1, q ), 1,
869 $ a( 1, p ), 1 )
870 IF( rsvec ) THEN
871 CALL saxpy( mvl,
872 $ t*apoaq,
873 $ v( 1, p ), 1,
874 $ v( 1, q ), 1 )
875 CALL saxpy( mvl,
876 $ -cs*sn*aqoap,
877 $ v( 1, q ), 1,
878 $ v( 1, p ), 1 )
879 END IF
880 d( p ) = d( p ) / cs
881 d( q ) = d( q )*cs
882 ELSE
883 IF( d( p ).GE.d( q ) ) THEN
884 CALL saxpy( m, -t*aqoap,
885 $ a( 1, q ), 1,
886 $ a( 1, p ), 1 )
887 CALL saxpy( m,
888 $ cs*sn*apoaq,
889 $ a( 1, p ), 1,
890 $ a( 1, q ), 1 )
891 d( p ) = d( p )*cs
892 d( q ) = d( q ) / cs
893 IF( rsvec ) THEN
894 CALL saxpy( mvl,
895 $ -t*aqoap,
896 $ v( 1, q ), 1,
897 $ v( 1, p ), 1 )
898 CALL saxpy( mvl,
899 $ cs*sn*apoaq,
900 $ v( 1, p ), 1,
901 $ v( 1, q ), 1 )
902 END IF
903 ELSE
904 CALL saxpy( m, t*apoaq,
905 $ a( 1, p ), 1,
906 $ a( 1, q ), 1 )
907 CALL saxpy( m,
908 $ -cs*sn*aqoap,
909 $ a( 1, q ), 1,
910 $ a( 1, p ), 1 )
911 d( p ) = d( p ) / cs
912 d( q ) = d( q )*cs
913 IF( rsvec ) THEN
914 CALL saxpy( mvl,
915 $ t*apoaq, v( 1, p ),
916 $ 1, v( 1, q ), 1 )
917 CALL saxpy( mvl,
918 $ -cs*sn*aqoap,
919 $ v( 1, q ), 1,
920 $ v( 1, p ), 1 )
921 END IF
922 END IF
923 END IF
924 END IF
925 END IF
926*
927 ELSE
928 IF( aapp.GT.aaqq ) THEN
929 CALL scopy( m, a( 1, p ), 1,
930 $ work,
931 $ 1 )
932 CALL slascl( 'G', 0, 0, aapp,
933 $ one,
934 $ m, 1, work, lda, ierr )
935 CALL slascl( 'G', 0, 0, aaqq,
936 $ one,
937 $ m, 1, a( 1, q ), lda,
938 $ ierr )
939 temp1 = -aapq*d( p ) / d( q )
940 CALL saxpy( m, temp1, work, 1,
941 $ a( 1, q ), 1 )
942 CALL slascl( 'G', 0, 0, one,
943 $ aaqq,
944 $ m, 1, a( 1, q ), lda,
945 $ ierr )
946 sva( q ) = aaqq*sqrt( max( zero,
947 $ one-aapq*aapq ) )
948 mxsinj = max( mxsinj, sfmin )
949 ELSE
950 CALL scopy( m, a( 1, q ), 1,
951 $ work,
952 $ 1 )
953 CALL slascl( 'G', 0, 0, aaqq,
954 $ one,
955 $ m, 1, work, lda, ierr )
956 CALL slascl( 'G', 0, 0, aapp,
957 $ one,
958 $ m, 1, a( 1, p ), lda,
959 $ ierr )
960 temp1 = -aapq*d( q ) / d( p )
961 CALL saxpy( m, temp1, work, 1,
962 $ a( 1, p ), 1 )
963 CALL slascl( 'G', 0, 0, one,
964 $ aapp,
965 $ m, 1, a( 1, p ), lda,
966 $ ierr )
967 sva( p ) = aapp*sqrt( max( zero,
968 $ one-aapq*aapq ) )
969 mxsinj = max( mxsinj, sfmin )
970 END IF
971 END IF
972* END IF ROTOK THEN ... ELSE
973*
974* In the case of cancellation in updating SVA(q)
975* .. recompute SVA(q)
976 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
977 $ THEN
978 IF( ( aaqq.LT.rootbig ) .AND.
979 $ ( aaqq.GT.rootsfmin ) ) THEN
980 sva( q ) = snrm2( m, a( 1, q ),
981 $ 1 )*
982 $ d( q )
983 ELSE
984 t = zero
985 aaqq = one
986 CALL slassq( m, a( 1, q ), 1, t,
987 $ aaqq )
988 sva( q ) = t*sqrt( aaqq )*d( q )
989 END IF
990 END IF
991 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
992 IF( ( aapp.LT.rootbig ) .AND.
993 $ ( aapp.GT.rootsfmin ) ) THEN
994 aapp = snrm2( m, a( 1, p ), 1 )*
995 $ d( p )
996 ELSE
997 t = zero
998 aapp = one
999 CALL slassq( m, a( 1, p ), 1, t,
1000 $ aapp )
1001 aapp = t*sqrt( aapp )*d( p )
1002 END IF
1003 sva( p ) = aapp
1004 END IF
1005* end of OK rotation
1006 ELSE
1007 notrot = notrot + 1
1008 pskipped = pskipped + 1
1009 ijblsk = ijblsk + 1
1010 END IF
1011 ELSE
1012 notrot = notrot + 1
1013 pskipped = pskipped + 1
1014 ijblsk = ijblsk + 1
1015 END IF
1016*
1017 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
1018 $ THEN
1019 sva( p ) = aapp
1020 notrot = 0
1021 GO TO 2011
1022 END IF
1023 IF( ( i.LE.swband ) .AND.
1024 $ ( pskipped.GT.rowskip ) ) THEN
1025 aapp = -aapp
1026 notrot = 0
1027 GO TO 2203
1028 END IF
1029*
1030 2200 CONTINUE
1031* end of the q-loop
1032 2203 CONTINUE
1033*
1034 sva( p ) = aapp
1035*
1036 ELSE
1037 IF( aapp.EQ.zero )notrot = notrot +
1038 $ min( jgl+kbl-1, n ) - jgl + 1
1039 IF( aapp.LT.zero )notrot = 0
1040 END IF
1041
1042 2100 CONTINUE
1043* end of the p-loop
1044 2010 CONTINUE
1045* end of the jbc-loop
1046 2011 CONTINUE
1047*2011 bailed out of the jbc-loop
1048 DO 2012 p = igl, min( igl+kbl-1, n )
1049 sva( p ) = abs( sva( p ) )
1050 2012 CONTINUE
1051*
1052 2000 CONTINUE
1053*2000 :: end of the ibr-loop
1054*
1055* .. update SVA(N)
1056 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1057 $ THEN
1058 sva( n ) = snrm2( m, a( 1, n ), 1 )*d( n )
1059 ELSE
1060 t = zero
1061 aapp = one
1062 CALL slassq( m, a( 1, n ), 1, t, aapp )
1063 sva( n ) = t*sqrt( aapp )*d( n )
1064 END IF
1065*
1066* Additional steering devices
1067*
1068 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
1069 $ ( iswrot.LE.n ) ) )swband = i
1070*
1071 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.float( n )*tol ) .AND.
1072 $ ( float( n )*mxaapq*mxsinj.LT.tol ) ) THEN
1073 GO TO 1994
1074 END IF
1075*
1076 IF( notrot.GE.emptsw )GO TO 1994
1077
1078 1993 CONTINUE
1079* end i=1:NSWEEP loop
1080* #:) Reaching this point means that the procedure has completed the given
1081* number of iterations.
1082 info = nsweep - 1
1083 GO TO 1995
1084 1994 CONTINUE
1085* #:) Reaching this point means that during the i-th sweep all pivots were
1086* below the given tolerance, causing early exit.
1087*
1088 info = 0
1089* #:) INFO = 0 confirms successful iterations.
1090 1995 CONTINUE
1091*
1092* Sort the vector D.
1093 DO 5991 p = 1, n - 1
1094 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
1095 IF( p.NE.q ) THEN
1096 temp1 = sva( p )
1097 sva( p ) = sva( q )
1098 sva( q ) = temp1
1099 temp1 = d( p )
1100 d( p ) = d( q )
1101 d( q ) = temp1
1102 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1103 IF( rsvec )CALL sswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
1104 END IF
1105 5991 CONTINUE
1106*
1107 RETURN
1108* ..
1109* .. END OF SGSVJ0
1110* ..
1111 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgsvj0
subroutine sgsvj0(jobv, m, n, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
SGSVJ0 pre-processor for the routine sgesvj.
Definition sgsvj0.f:216
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:142
slassq
subroutine slassq(n, x, incx, scale, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:122
srotm
subroutine srotm(n, sx, incx, sy, incy, sparam)
SROTM
Definition srotm.f:97
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82