LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgesvj.f
Go to the documentation of this file.
1*> \brief \b SGESVJ
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGESVJ + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgesvj.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgesvj.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgesvj.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
22* LDV, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDV, LWORK, M, MV, N
26* CHARACTER*1 JOBA, JOBU, JOBV
27* ..
28* .. Array Arguments ..
29* REAL A( LDA, * ), SVA( N ), V( LDV, * ),
30* $ WORK( LWORK )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SGESVJ computes the singular value decomposition (SVD) of a real
40*> M-by-N matrix A, where M >= N. The SVD of A is written as
41*> [++] [xx] [x0] [xx]
42*> A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx]
43*> [++] [xx]
44*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
45*> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
46*> of SIGMA are the singular values of A. The columns of U and V are the
47*> left and the right singular vectors of A, respectively.
48*> SGESVJ can sometimes compute tiny singular values and their singular vectors much
49*> more accurately than other SVD routines, see below under Further Details.
50*> \endverbatim
51*
52* Arguments:
53* ==========
54*
55*> \param[in] JOBA
56*> \verbatim
57*> JOBA is CHARACTER*1
58*> Specifies the structure of A.
59*> = 'L': The input matrix A is lower triangular;
60*> = 'U': The input matrix A is upper triangular;
61*> = 'G': The input matrix A is general M-by-N matrix, M >= N.
62*> \endverbatim
63*>
64*> \param[in] JOBU
65*> \verbatim
66*> JOBU is CHARACTER*1
67*> Specifies whether to compute the left singular vectors
68*> (columns of U):
69*> = 'U': The left singular vectors corresponding to the nonzero
70*> singular values are computed and returned in the leading
71*> columns of A. See more details in the description of A.
72*> The default numerical orthogonality threshold is set to
73*> approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=SLAMCH('E').
74*> = 'C': Analogous to JOBU='U', except that user can control the
75*> level of numerical orthogonality of the computed left
76*> singular vectors. TOL can be set to TOL = CTOL*EPS, where
77*> CTOL is given on input in the array WORK.
78*> No CTOL smaller than ONE is allowed. CTOL greater
79*> than 1 / EPS is meaningless. The option 'C'
80*> can be used if M*EPS is satisfactory orthogonality
81*> of the computed left singular vectors, so CTOL=M could
82*> save few sweeps of Jacobi rotations.
83*> See the descriptions of A and WORK(1).
84*> = 'N': The matrix U is not computed. However, see the
85*> description of A.
86*> \endverbatim
87*>
88*> \param[in] JOBV
89*> \verbatim
90*> JOBV is CHARACTER*1
91*> Specifies whether to compute the right singular vectors, that
92*> is, the matrix V:
93*> = 'V': the matrix V is computed and returned in the array V
94*> = 'A': the Jacobi rotations are applied to the MV-by-N
95*> array V. In other words, the right singular vector
96*> matrix V is not computed explicitly; instead it is
97*> applied to an MV-by-N matrix initially stored in the
98*> first MV rows of V.
99*> = 'N': the matrix V is not computed and the array V is not
100*> referenced
101*> \endverbatim
102*>
103*> \param[in] M
104*> \verbatim
105*> M is INTEGER
106*> The number of rows of the input matrix A. 1/SLAMCH('E') > M >= 0.
107*> \endverbatim
108*>
109*> \param[in] N
110*> \verbatim
111*> N is INTEGER
112*> The number of columns of the input matrix A.
113*> M >= N >= 0.
114*> \endverbatim
115*>
116*> \param[in,out] A
117*> \verbatim
118*> A is REAL array, dimension (LDA,N)
119*> On entry, the M-by-N matrix A.
120*> On exit,
121*> If JOBU = 'U' .OR. JOBU = 'C':
122*> If INFO = 0:
123*> RANKA orthonormal columns of U are returned in the
124*> leading RANKA columns of the array A. Here RANKA <= N
125*> is the number of computed singular values of A that are
126*> above the underflow threshold SLAMCH('S'). The singular
127*> vectors corresponding to underflowed or zero singular
128*> values are not computed. The value of RANKA is returned
129*> in the array WORK as RANKA=NINT(WORK(2)). Also see the
130*> descriptions of SVA and WORK. The computed columns of U
131*> are mutually numerically orthogonal up to approximately
132*> TOL=SQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
133*> see the description of JOBU.
134*> If INFO > 0,
135*> the procedure SGESVJ did not converge in the given number
136*> of iterations (sweeps). In that case, the computed
137*> columns of U may not be orthogonal up to TOL. The output
138*> U (stored in A), SIGMA (given by the computed singular
139*> values in SVA(1:N)) and V is still a decomposition of the
140*> input matrix A in the sense that the residual
141*> ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
142*> If JOBU = 'N':
143*> If INFO = 0:
144*> Note that the left singular vectors are 'for free' in the
145*> one-sided Jacobi SVD algorithm. However, if only the
146*> singular values are needed, the level of numerical
147*> orthogonality of U is not an issue and iterations are
148*> stopped when the columns of the iterated matrix are
149*> numerically orthogonal up to approximately M*EPS. Thus,
150*> on exit, A contains the columns of U scaled with the
151*> corresponding singular values.
152*> If INFO > 0:
153*> the procedure SGESVJ did not converge in the given number
154*> of iterations (sweeps).
155*> \endverbatim
156*>
157*> \param[in] LDA
158*> \verbatim
159*> LDA is INTEGER
160*> The leading dimension of the array A. LDA >= max(1,M).
161*> \endverbatim
162*>
163*> \param[out] SVA
164*> \verbatim
165*> SVA is REAL array, dimension (N)
166*> On exit,
167*> If INFO = 0 :
168*> depending on the value SCALE = WORK(1), we have:
169*> If SCALE = ONE:
170*> SVA(1:N) contains the computed singular values of A.
171*> During the computation SVA contains the Euclidean column
172*> norms of the iterated matrices in the array A.
173*> If SCALE .NE. ONE:
174*> The singular values of A are SCALE*SVA(1:N), and this
175*> factored representation is due to the fact that some of the
176*> singular values of A might underflow or overflow.
177*>
178*> If INFO > 0 :
179*> the procedure SGESVJ did not converge in the given number of
180*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
181*> \endverbatim
182*>
183*> \param[in] MV
184*> \verbatim
185*> MV is INTEGER
186*> If JOBV = 'A', then the product of Jacobi rotations in SGESVJ
187*> is applied to the first MV rows of V. See the description of JOBV.
188*> \endverbatim
189*>
190*> \param[in,out] V
191*> \verbatim
192*> V is REAL array, dimension (LDV,N)
193*> If JOBV = 'V', then V contains on exit the N-by-N matrix of
194*> the right singular vectors;
195*> If JOBV = 'A', then V contains the product of the computed right
196*> singular vector matrix and the initial matrix in
197*> the array V.
198*> If JOBV = 'N', then V is not referenced.
199*> \endverbatim
200*>
201*> \param[in] LDV
202*> \verbatim
203*> LDV is INTEGER
204*> The leading dimension of the array V, LDV >= 1.
205*> If JOBV = 'V', then LDV >= max(1,N).
206*> If JOBV = 'A', then LDV >= max(1,MV) .
207*> \endverbatim
208*>
209*> \param[in,out] WORK
210*> \verbatim
211*> WORK is REAL array, dimension (LWORK)
212*> On entry,
213*> If JOBU = 'C' :
214*> WORK(1) = CTOL, where CTOL defines the threshold for convergence.
215*> The process stops if all columns of A are mutually
216*> orthogonal up to CTOL*EPS, EPS=SLAMCH('E').
217*> It is required that CTOL >= ONE, i.e. it is not
218*> allowed to force the routine to obtain orthogonality
219*> below EPSILON.
220*> On exit,
221*> WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
222*> are the computed singular vcalues of A.
223*> (See description of SVA().)
224*> WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
225*> singular values.
226*> WORK(3) = NINT(WORK(3)) is the number of the computed singular
227*> values that are larger than the underflow threshold.
228*> WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
229*> rotations needed for numerical convergence.
230*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
231*> This is useful information in cases when SGESVJ did
232*> not converge, as it can be used to estimate whether
233*> the output is still useful and for post festum analysis.
234*> WORK(6) = the largest absolute value over all sines of the
235*> Jacobi rotation angles in the last sweep. It can be
236*> useful for a post festum analysis.
237*> \endverbatim
238*>
239*> \param[in] LWORK
240*> \verbatim
241*> LWORK is INTEGER
242*> length of WORK, WORK >= MAX(6,M+N)
243*> \endverbatim
244*>
245*> \param[out] INFO
246*> \verbatim
247*> INFO is INTEGER
248*> = 0: successful exit.
249*> < 0: if INFO = -i, then the i-th argument had an illegal value
250*> > 0: SGESVJ did not converge in the maximal allowed number (30)
251*> of sweeps. The output may still be useful. See the
252*> description of WORK.
253*> \endverbatim
254*
255* Authors:
256* ========
257*
258*> \author Univ. of Tennessee
259*> \author Univ. of California Berkeley
260*> \author Univ. of Colorado Denver
261*> \author NAG Ltd.
262*
263*> \ingroup realGEcomputational
264*
265*> \par Further Details:
266* =====================
267*>
268*> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
269*> rotations. The rotations are implemented as fast scaled rotations of
270*> Anda and Park [1]. In the case of underflow of the Jacobi angle, a
271*> modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
272*> column interchanges of de Rijk [2]. The relative accuracy of the computed
273*> singular values and the accuracy of the computed singular vectors (in
274*> angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
275*> The condition number that determines the accuracy in the full rank case
276*> is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
277*> spectral condition number. The best performance of this Jacobi SVD
278*> procedure is achieved if used in an accelerated version of Drmac and
279*> Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
280*> Some tuning parameters (marked with [TP]) are available for the
281*> implementer. \n
282*> The computational range for the nonzero singular values is the machine
283*> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
284*> denormalized singular values can be computed with the corresponding
285*> gradual loss of accurate digits.
286*>
287*> \par Contributors:
288* ==================
289*>
290*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
291*>
292*> \par References:
293* ================
294*>
295*> [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. \n
296*> SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. \n\n
297*> [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
298*> singular value decomposition on a vector computer. \n
299*> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. \n\n
300*> [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. \n
301*> [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
302*> value computation in floating point arithmetic. \n
303*> SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. \n\n
304*> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. \n
305*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. \n
306*> LAPACK Working note 169. \n\n
307*> [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. \n
308*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. \n
309*> LAPACK Working note 170. \n\n
310*> [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
311*> QSVD, (H,K)-SVD computations.\n
312*> Department of Mathematics, University of Zagreb, 2008.
313*>
314*> \par Bugs, Examples and Comments:
315* =================================
316*>
317*> Please report all bugs and send interesting test examples and comments to
318*> drmac@math.hr. Thank you.
319*
320* =====================================================================
321 SUBROUTINE sgesvj( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
322 $ LDV, WORK, LWORK, INFO )
323*
324* -- LAPACK computational routine --
325* -- LAPACK is a software package provided by Univ. of Tennessee, --
326* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
327*
328* .. Scalar Arguments ..
329 INTEGER INFO, LDA, LDV, LWORK, M, MV, N
330 CHARACTER*1 JOBA, JOBU, JOBV
331* ..
332* .. Array Arguments ..
333 REAL A( LDA, * ), SVA( N ), V( LDV, * ),
334 $ work( lwork )
335* ..
336*
337* =====================================================================
338*
339* .. Local Parameters ..
340 REAL ZERO, HALF, ONE
341 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
342 INTEGER NSWEEP
343 parameter( nsweep = 30 )
344* ..
345* .. Local Scalars ..
346 REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
347 $ bigtheta, cs, ctol, epsln, large, mxaapq,
348 $ mxsinj, rootbig, rooteps, rootsfmin, roottol,
349 $ skl, sfmin, small, sn, t, temp1, theta,
350 $ thsign, tol
351 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
352 $ iswrot, jbc, jgl, kbl, lkahead, mvl, n2, n34,
353 $ n4, nbl, notrot, p, pskipped, q, rowskip,
354 $ swband
355 LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
356 $ rsvec, uctol, upper
357* ..
358* .. Local Arrays ..
359 REAL FASTR( 5 )
360* ..
361* .. Intrinsic Functions ..
362 INTRINSIC abs, max, min, float, sign, sqrt
363* ..
364* .. External Functions ..
365* ..
366* from BLAS
367 REAL SDOT, SNRM2
368 EXTERNAL sdot, snrm2
369 INTEGER ISAMAX
370 EXTERNAL isamax
371* from LAPACK
372 REAL SLAMCH
373 EXTERNAL slamch
374 LOGICAL LSAME
375 EXTERNAL lsame
376* ..
377* .. External Subroutines ..
378* ..
379* from BLAS
380 EXTERNAL saxpy, scopy, srotm, sscal, sswap
381* from LAPACK
382 EXTERNAL slascl, slaset, slassq, xerbla
383*
384 EXTERNAL sgsvj0, sgsvj1
385* ..
386* .. Executable Statements ..
387*
388* Test the input arguments
389*
390 lsvec = lsame( jobu, 'U' )
391 uctol = lsame( jobu, 'C' )
392 rsvec = lsame( jobv, 'V' )
393 applv = lsame( jobv, 'A' )
394 upper = lsame( joba, 'U' )
395 lower = lsame( joba, 'L' )
396*
397 IF( .NOT.( upper .OR. lower .OR. lsame( joba, 'G' ) ) ) THEN
398 info = -1
399 ELSE IF( .NOT.( lsvec .OR. uctol .OR. lsame( jobu, 'N' ) ) ) THEN
400 info = -2
401 ELSE IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
402 info = -3
403 ELSE IF( m.LT.0 ) THEN
404 info = -4
405 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
406 info = -5
407 ELSE IF( lda.LT.m ) THEN
408 info = -7
409 ELSE IF( mv.LT.0 ) THEN
410 info = -9
411 ELSE IF( ( rsvec .AND. ( ldv.LT.n ) ) .OR.
412 $ ( applv .AND. ( ldv.LT.mv ) ) ) THEN
413 info = -11
414 ELSE IF( uctol .AND. ( work( 1 ).LE.one ) ) THEN
415 info = -12
416 ELSE IF( lwork.LT.max( m+n, 6 ) ) THEN
417 info = -13
418 ELSE
419 info = 0
420 END IF
421*
422* #:(
423 IF( info.NE.0 ) THEN
424 CALL xerbla( 'SGESVJ', -info )
425 RETURN
426 END IF
427*
428* #:) Quick return for void matrix
429*
430 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )RETURN
431*
432* Set numerical parameters
433* The stopping criterion for Jacobi rotations is
434*
435* max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
436*
437* where EPS is the round-off and CTOL is defined as follows:
438*
439 IF( uctol ) THEN
440* ... user controlled
441 ctol = work( 1 )
442 ELSE
443* ... default
444 IF( lsvec .OR. rsvec .OR. applv ) THEN
445 ctol = sqrt( float( m ) )
446 ELSE
447 ctol = float( m )
448 END IF
449 END IF
450* ... and the machine dependent parameters are
451*[!] (Make sure that SLAMCH() works properly on the target machine.)
452*
453 epsln = slamch( 'Epsilon' )
454 rooteps = sqrt( epsln )
455 sfmin = slamch( 'SafeMinimum' )
456 rootsfmin = sqrt( sfmin )
457 small = sfmin / epsln
458 big = slamch( 'Overflow' )
459* BIG = ONE / SFMIN
460 rootbig = one / rootsfmin
461 large = big / sqrt( float( m*n ) )
462 bigtheta = one / rooteps
463*
464 tol = ctol*epsln
465 roottol = sqrt( tol )
466*
467 IF( float( m )*epsln.GE.one ) THEN
468 info = -4
469 CALL xerbla( 'SGESVJ', -info )
470 RETURN
471 END IF
472*
473* Initialize the right singular vector matrix.
474*
475 IF( rsvec ) THEN
476 mvl = n
477 CALL slaset( 'A', mvl, n, zero, one, v, ldv )
478 ELSE IF( applv ) THEN
479 mvl = mv
480 END IF
481 rsvec = rsvec .OR. applv
482*
483* Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
484*(!) If necessary, scale A to protect the largest singular value
485* from overflow. It is possible that saving the largest singular
486* value destroys the information about the small ones.
487* This initial scaling is almost minimal in the sense that the
488* goal is to make sure that no column norm overflows, and that
489* SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
490* in A are detected, the procedure returns with INFO=-6.
491*
492 skl = one / sqrt( float( m )*float( n ) )
493 noscale = .true.
494 goscale = .true.
495*
496 IF( lower ) THEN
497* the input matrix is M-by-N lower triangular (trapezoidal)
498 DO 1874 p = 1, n
499 aapp = zero
500 aaqq = one
501 CALL slassq( m-p+1, a( p, p ), 1, aapp, aaqq )
502 IF( aapp.GT.big ) THEN
503 info = -6
504 CALL xerbla( 'SGESVJ', -info )
505 RETURN
506 END IF
507 aaqq = sqrt( aaqq )
508 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale ) THEN
509 sva( p ) = aapp*aaqq
510 ELSE
511 noscale = .false.
512 sva( p ) = aapp*( aaqq*skl )
513 IF( goscale ) THEN
514 goscale = .false.
515 DO 1873 q = 1, p - 1
516 sva( q ) = sva( q )*skl
517 1873 CONTINUE
518 END IF
519 END IF
520 1874 CONTINUE
521 ELSE IF( upper ) THEN
522* the input matrix is M-by-N upper triangular (trapezoidal)
523 DO 2874 p = 1, n
524 aapp = zero
525 aaqq = one
526 CALL slassq( p, a( 1, p ), 1, aapp, aaqq )
527 IF( aapp.GT.big ) THEN
528 info = -6
529 CALL xerbla( 'SGESVJ', -info )
530 RETURN
531 END IF
532 aaqq = sqrt( aaqq )
533 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale ) THEN
534 sva( p ) = aapp*aaqq
535 ELSE
536 noscale = .false.
537 sva( p ) = aapp*( aaqq*skl )
538 IF( goscale ) THEN
539 goscale = .false.
540 DO 2873 q = 1, p - 1
541 sva( q ) = sva( q )*skl
542 2873 CONTINUE
543 END IF
544 END IF
545 2874 CONTINUE
546 ELSE
547* the input matrix is M-by-N general dense
548 DO 3874 p = 1, n
549 aapp = zero
550 aaqq = one
551 CALL slassq( m, a( 1, p ), 1, aapp, aaqq )
552 IF( aapp.GT.big ) THEN
553 info = -6
554 CALL xerbla( 'SGESVJ', -info )
555 RETURN
556 END IF
557 aaqq = sqrt( aaqq )
558 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale ) THEN
559 sva( p ) = aapp*aaqq
560 ELSE
561 noscale = .false.
562 sva( p ) = aapp*( aaqq*skl )
563 IF( goscale ) THEN
564 goscale = .false.
565 DO 3873 q = 1, p - 1
566 sva( q ) = sva( q )*skl
567 3873 CONTINUE
568 END IF
569 END IF
570 3874 CONTINUE
571 END IF
572*
573 IF( noscale )skl = one
574*
575* Move the smaller part of the spectrum from the underflow threshold
576*(!) Start by determining the position of the nonzero entries of the
577* array SVA() relative to ( SFMIN, BIG ).
578*
579 aapp = zero
580 aaqq = big
581 DO 4781 p = 1, n
582 IF( sva( p ).NE.zero )aaqq = min( aaqq, sva( p ) )
583 aapp = max( aapp, sva( p ) )
584 4781 CONTINUE
585*
586* #:) Quick return for zero matrix
587*
588 IF( aapp.EQ.zero ) THEN
589 IF( lsvec )CALL slaset( 'G', m, n, zero, one, a, lda )
590 work( 1 ) = one
591 work( 2 ) = zero
592 work( 3 ) = zero
593 work( 4 ) = zero
594 work( 5 ) = zero
595 work( 6 ) = zero
596 RETURN
597 END IF
598*
599* #:) Quick return for one-column matrix
600*
601 IF( n.EQ.1 ) THEN
602 IF( lsvec )CALL slascl( 'G', 0, 0, sva( 1 ), skl, m, 1,
603 $ a( 1, 1 ), lda, ierr )
604 work( 1 ) = one / skl
605 IF( sva( 1 ).GE.sfmin ) THEN
606 work( 2 ) = one
607 ELSE
608 work( 2 ) = zero
609 END IF
610 work( 3 ) = zero
611 work( 4 ) = zero
612 work( 5 ) = zero
613 work( 6 ) = zero
614 RETURN
615 END IF
616*
617* Protect small singular values from underflow, and try to
618* avoid underflows/overflows in computing Jacobi rotations.
619*
620 sn = sqrt( sfmin / epsln )
621 temp1 = sqrt( big / float( n ) )
622 IF( ( aapp.LE.sn ) .OR. ( aaqq.GE.temp1 ) .OR.
623 $ ( ( sn.LE.aaqq ) .AND. ( aapp.LE.temp1 ) ) ) THEN
624 temp1 = min( big, temp1 / aapp )
625* AAQQ = AAQQ*TEMP1
626* AAPP = AAPP*TEMP1
627 ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.LE.temp1 ) ) THEN
628 temp1 = min( sn / aaqq, big / ( aapp*sqrt( float( n ) ) ) )
629* AAQQ = AAQQ*TEMP1
630* AAPP = AAPP*TEMP1
631 ELSE IF( ( aaqq.GE.sn ) .AND. ( aapp.GE.temp1 ) ) THEN
632 temp1 = max( sn / aaqq, temp1 / aapp )
633* AAQQ = AAQQ*TEMP1
634* AAPP = AAPP*TEMP1
635 ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.GE.temp1 ) ) THEN
636 temp1 = min( sn / aaqq, big / ( sqrt( float( n ) )*aapp ) )
637* AAQQ = AAQQ*TEMP1
638* AAPP = AAPP*TEMP1
639 ELSE
640 temp1 = one
641 END IF
642*
643* Scale, if necessary
644*
645 IF( temp1.NE.one ) THEN
646 CALL slascl( 'G', 0, 0, one, temp1, n, 1, sva, n, ierr )
647 END IF
648 skl = temp1*skl
649 IF( skl.NE.one ) THEN
650 CALL slascl( joba, 0, 0, one, skl, m, n, a, lda, ierr )
651 skl = one / skl
652 END IF
653*
654* Row-cyclic Jacobi SVD algorithm with column pivoting
655*
656 emptsw = ( n*( n-1 ) ) / 2
657 notrot = 0
658 fastr( 1 ) = zero
659*
660* A is represented in factored form A = A * diag(WORK), where diag(WORK)
661* is initialized to identity. WORK is updated during fast scaled
662* rotations.
663*
664 DO 1868 q = 1, n
665 work( q ) = one
666 1868 CONTINUE
667*
668*
669 swband = 3
670*[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
671* if SGESVJ is used as a computational routine in the preconditioned
672* Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
673* works on pivots inside a band-like region around the diagonal.
674* The boundaries are determined dynamically, based on the number of
675* pivots above a threshold.
676*
677 kbl = min( 8, n )
678*[TP] KBL is a tuning parameter that defines the tile size in the
679* tiling of the p-q loops of pivot pairs. In general, an optimal
680* value of KBL depends on the matrix dimensions and on the
681* parameters of the computer's memory.
682*
683 nbl = n / kbl
684 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
685*
686 blskip = kbl**2
687*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
688*
689 rowskip = min( 5, kbl )
690*[TP] ROWSKIP is a tuning parameter.
691*
692 lkahead = 1
693*[TP] LKAHEAD is a tuning parameter.
694*
695* Quasi block transformations, using the lower (upper) triangular
696* structure of the input matrix. The quasi-block-cycling usually
697* invokes cubic convergence. Big part of this cycle is done inside
698* canonical subspaces of dimensions less than M.
699*
700 IF( ( lower .OR. upper ) .AND. ( n.GT.max( 64, 4*kbl ) ) ) THEN
701*[TP] The number of partition levels and the actual partition are
702* tuning parameters.
703 n4 = n / 4
704 n2 = n / 2
705 n34 = 3*n4
706 IF( applv ) THEN
707 q = 0
708 ELSE
709 q = 1
710 END IF
711*
712 IF( lower ) THEN
713*
714* This works very well on lower triangular matrices, in particular
715* in the framework of the preconditioned Jacobi SVD (xGEJSV).
716* The idea is simple:
717* [+ 0 0 0] Note that Jacobi transformations of [0 0]
718* [+ + 0 0] [0 0]
719* [+ + x 0] actually work on [x 0] [x 0]
720* [+ + x x] [x x]. [x x]
721*
722 CALL sgsvj0( jobv, m-n34, n-n34, a( n34+1, n34+1 ), lda,
723 $ work( n34+1 ), sva( n34+1 ), mvl,
724 $ v( n34*q+1, n34+1 ), ldv, epsln, sfmin, tol,
725 $ 2, work( n+1 ), lwork-n, ierr )
726*
727 CALL sgsvj0( jobv, m-n2, n34-n2, a( n2+1, n2+1 ), lda,
728 $ work( n2+1 ), sva( n2+1 ), mvl,
729 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 2,
730 $ work( n+1 ), lwork-n, ierr )
731*
732 CALL sgsvj1( jobv, m-n2, n-n2, n4, a( n2+1, n2+1 ), lda,
733 $ work( n2+1 ), sva( n2+1 ), mvl,
734 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
735 $ work( n+1 ), lwork-n, ierr )
736*
737 CALL sgsvj0( jobv, m-n4, n2-n4, a( n4+1, n4+1 ), lda,
738 $ work( n4+1 ), sva( n4+1 ), mvl,
739 $ v( n4*q+1, n4+1 ), ldv, epsln, sfmin, tol, 1,
740 $ work( n+1 ), lwork-n, ierr )
741*
742 CALL sgsvj0( jobv, m, n4, a, lda, work, sva, mvl, v, ldv,
743 $ epsln, sfmin, tol, 1, work( n+1 ), lwork-n,
744 $ ierr )
745*
746 CALL sgsvj1( jobv, m, n2, n4, a, lda, work, sva, mvl, v,
747 $ ldv, epsln, sfmin, tol, 1, work( n+1 ),
748 $ lwork-n, ierr )
749*
750*
751 ELSE IF( upper ) THEN
752*
753*
754 CALL sgsvj0( jobv, n4, n4, a, lda, work, sva, mvl, v, ldv,
755 $ epsln, sfmin, tol, 2, work( n+1 ), lwork-n,
756 $ ierr )
757*
758 CALL sgsvj0( jobv, n2, n4, a( 1, n4+1 ), lda, work( n4+1 ),
759 $ sva( n4+1 ), mvl, v( n4*q+1, n4+1 ), ldv,
760 $ epsln, sfmin, tol, 1, work( n+1 ), lwork-n,
761 $ ierr )
762*
763 CALL sgsvj1( jobv, n2, n2, n4, a, lda, work, sva, mvl, v,
764 $ ldv, epsln, sfmin, tol, 1, work( n+1 ),
765 $ lwork-n, ierr )
766*
767 CALL sgsvj0( jobv, n2+n4, n4, a( 1, n2+1 ), lda,
768 $ work( n2+1 ), sva( n2+1 ), mvl,
769 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
770 $ work( n+1 ), lwork-n, ierr )
771
772 END IF
773*
774 END IF
775*
776* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
777*
778 DO 1993 i = 1, nsweep
779*
780* .. go go go ...
781*
782 mxaapq = zero
783 mxsinj = zero
784 iswrot = 0
785*
786 notrot = 0
787 pskipped = 0
788*
789* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
790* 1 <= p < q <= N. This is the first step toward a blocked implementation
791* of the rotations. New implementation, based on block transformations,
792* is under development.
793*
794 DO 2000 ibr = 1, nbl
795*
796 igl = ( ibr-1 )*kbl + 1
797*
798 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
799*
800 igl = igl + ir1*kbl
801*
802 DO 2001 p = igl, min( igl+kbl-1, n-1 )
803*
804* .. de Rijk's pivoting
805*
806 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
807 IF( p.NE.q ) THEN
808 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
809 IF( rsvec )CALL sswap( mvl, v( 1, p ), 1,
810 $ v( 1, q ), 1 )
811 temp1 = sva( p )
812 sva( p ) = sva( q )
813 sva( q ) = temp1
814 temp1 = work( p )
815 work( p ) = work( q )
816 work( q ) = temp1
817 END IF
818*
819 IF( ir1.EQ.0 ) THEN
820*
821* Column norms are periodically updated by explicit
822* norm computation.
823* Caveat:
824* Unfortunately, some BLAS implementations compute SNRM2(M,A(1,p),1)
825* as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
826* overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
827* underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
828* Hence, SNRM2 cannot be trusted, not even in the case when
829* the true norm is far from the under(over)flow boundaries.
830* If properly implemented SNRM2 is available, the IF-THEN-ELSE
831* below should read "AAPP = SNRM2( M, A(1,p), 1 ) * WORK(p)".
832*
833 IF( ( sva( p ).LT.rootbig ) .AND.
834 $ ( sva( p ).GT.rootsfmin ) ) THEN
835 sva( p ) = snrm2( m, a( 1, p ), 1 )*work( p )
836 ELSE
837 temp1 = zero
838 aapp = one
839 CALL slassq( m, a( 1, p ), 1, temp1, aapp )
840 sva( p ) = temp1*sqrt( aapp )*work( p )
841 END IF
842 aapp = sva( p )
843 ELSE
844 aapp = sva( p )
845 END IF
846*
847 IF( aapp.GT.zero ) THEN
848*
849 pskipped = 0
850*
851 DO 2002 q = p + 1, min( igl+kbl-1, n )
852*
853 aaqq = sva( q )
854*
855 IF( aaqq.GT.zero ) THEN
856*
857 aapp0 = aapp
858 IF( aaqq.GE.one ) THEN
859 rotok = ( small*aapp ).LE.aaqq
860 IF( aapp.LT.( big / aaqq ) ) THEN
861 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
862 $ q ), 1 )*work( p )*work( q ) /
863 $ aaqq ) / aapp
864 ELSE
865 CALL scopy( m, a( 1, p ), 1,
866 $ work( n+1 ), 1 )
867 CALL slascl( 'G', 0, 0, aapp,
868 $ work( p ), m, 1,
869 $ work( n+1 ), lda, ierr )
870 aapq = sdot( m, work( n+1 ), 1,
871 $ a( 1, q ), 1 )*work( q ) / aaqq
872 END IF
873 ELSE
874 rotok = aapp.LE.( aaqq / small )
875 IF( aapp.GT.( small / aaqq ) ) THEN
876 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
877 $ q ), 1 )*work( p )*work( q ) /
878 $ aaqq ) / aapp
879 ELSE
880 CALL scopy( m, a( 1, q ), 1,
881 $ work( n+1 ), 1 )
882 CALL slascl( 'G', 0, 0, aaqq,
883 $ work( q ), m, 1,
884 $ work( n+1 ), lda, ierr )
885 aapq = sdot( m, work( n+1 ), 1,
886 $ a( 1, p ), 1 )*work( p ) / aapp
887 END IF
888 END IF
889*
890 mxaapq = max( mxaapq, abs( aapq ) )
891*
892* TO rotate or NOT to rotate, THAT is the question ...
893*
894 IF( abs( aapq ).GT.tol ) THEN
895*
896* .. rotate
897*[RTD] ROTATED = ROTATED + ONE
898*
899 IF( ir1.EQ.0 ) THEN
900 notrot = 0
901 pskipped = 0
902 iswrot = iswrot + 1
903 END IF
904*
905 IF( rotok ) THEN
906*
907 aqoap = aaqq / aapp
908 apoaq = aapp / aaqq
909 theta = -half*abs( aqoap-apoaq ) / aapq
910*
911 IF( abs( theta ).GT.bigtheta ) THEN
912*
913 t = half / theta
914 fastr( 3 ) = t*work( p ) / work( q )
915 fastr( 4 ) = -t*work( q ) /
916 $ work( p )
917 CALL srotm( m, a( 1, p ), 1,
918 $ a( 1, q ), 1, fastr )
919 IF( rsvec )CALL srotm( mvl,
920 $ v( 1, p ), 1,
921 $ v( 1, q ), 1,
922 $ fastr )
923 sva( q ) = aaqq*sqrt( max( zero,
924 $ one+t*apoaq*aapq ) )
925 aapp = aapp*sqrt( max( zero,
926 $ one-t*aqoap*aapq ) )
927 mxsinj = max( mxsinj, abs( t ) )
928*
929 ELSE
930*
931* .. choose correct signum for THETA and rotate
932*
933 thsign = -sign( one, aapq )
934 t = one / ( theta+thsign*
935 $ sqrt( one+theta*theta ) )
936 cs = sqrt( one / ( one+t*t ) )
937 sn = t*cs
938*
939 mxsinj = max( mxsinj, abs( sn ) )
940 sva( q ) = aaqq*sqrt( max( zero,
941 $ one+t*apoaq*aapq ) )
942 aapp = aapp*sqrt( max( zero,
943 $ one-t*aqoap*aapq ) )
944*
945 apoaq = work( p ) / work( q )
946 aqoap = work( q ) / work( p )
947 IF( work( p ).GE.one ) THEN
948 IF( work( q ).GE.one ) THEN
949 fastr( 3 ) = t*apoaq
950 fastr( 4 ) = -t*aqoap
951 work( p ) = work( p )*cs
952 work( q ) = work( q )*cs
953 CALL srotm( m, a( 1, p ), 1,
954 $ a( 1, q ), 1,
955 $ fastr )
956 IF( rsvec )CALL srotm( mvl,
957 $ v( 1, p ), 1, v( 1, q ),
958 $ 1, fastr )
959 ELSE
960 CALL saxpy( m, -t*aqoap,
961 $ a( 1, q ), 1,
962 $ a( 1, p ), 1 )
963 CALL saxpy( m, cs*sn*apoaq,
964 $ a( 1, p ), 1,
965 $ a( 1, q ), 1 )
966 work( p ) = work( p )*cs
967 work( q ) = work( q ) / cs
968 IF( rsvec ) THEN
969 CALL saxpy( mvl, -t*aqoap,
970 $ v( 1, q ), 1,
971 $ v( 1, p ), 1 )
972 CALL saxpy( mvl,
973 $ cs*sn*apoaq,
974 $ v( 1, p ), 1,
975 $ v( 1, q ), 1 )
976 END IF
977 END IF
978 ELSE
979 IF( work( q ).GE.one ) THEN
980 CALL saxpy( m, t*apoaq,
981 $ a( 1, p ), 1,
982 $ a( 1, q ), 1 )
983 CALL saxpy( m, -cs*sn*aqoap,
984 $ a( 1, q ), 1,
985 $ a( 1, p ), 1 )
986 work( p ) = work( p ) / cs
987 work( q ) = work( q )*cs
988 IF( rsvec ) THEN
989 CALL saxpy( mvl, t*apoaq,
990 $ v( 1, p ), 1,
991 $ v( 1, q ), 1 )
992 CALL saxpy( mvl,
993 $ -cs*sn*aqoap,
994 $ v( 1, q ), 1,
995 $ v( 1, p ), 1 )
996 END IF
997 ELSE
998 IF( work( p ).GE.work( q ) )
999 $ THEN
1000 CALL saxpy( m, -t*aqoap,
1001 $ a( 1, q ), 1,
1002 $ a( 1, p ), 1 )
1003 CALL saxpy( m, cs*sn*apoaq,
1004 $ a( 1, p ), 1,
1005 $ a( 1, q ), 1 )
1006 work( p ) = work( p )*cs
1007 work( q ) = work( q ) / cs
1008 IF( rsvec ) THEN
1009 CALL saxpy( mvl,
1010 $ -t*aqoap,
1011 $ v( 1, q ), 1,
1012 $ v( 1, p ), 1 )
1013 CALL saxpy( mvl,
1014 $ cs*sn*apoaq,
1015 $ v( 1, p ), 1,
1016 $ v( 1, q ), 1 )
1017 END IF
1018 ELSE
1019 CALL saxpy( m, t*apoaq,
1020 $ a( 1, p ), 1,
1021 $ a( 1, q ), 1 )
1022 CALL saxpy( m,
1023 $ -cs*sn*aqoap,
1024 $ a( 1, q ), 1,
1025 $ a( 1, p ), 1 )
1026 work( p ) = work( p ) / cs
1027 work( q ) = work( q )*cs
1028 IF( rsvec ) THEN
1029 CALL saxpy( mvl,
1030 $ t*apoaq, v( 1, p ),
1031 $ 1, v( 1, q ), 1 )
1032 CALL saxpy( mvl,
1033 $ -cs*sn*aqoap,
1034 $ v( 1, q ), 1,
1035 $ v( 1, p ), 1 )
1036 END IF
1037 END IF
1038 END IF
1039 END IF
1040 END IF
1041*
1042 ELSE
1043* .. have to use modified Gram-Schmidt like transformation
1044 CALL scopy( m, a( 1, p ), 1,
1045 $ work( n+1 ), 1 )
1046 CALL slascl( 'G', 0, 0, aapp, one, m,
1047 $ 1, work( n+1 ), lda,
1048 $ ierr )
1049 CALL slascl( 'G', 0, 0, aaqq, one, m,
1050 $ 1, a( 1, q ), lda, ierr )
1051 temp1 = -aapq*work( p ) / work( q )
1052 CALL saxpy( m, temp1, work( n+1 ), 1,
1053 $ a( 1, q ), 1 )
1054 CALL slascl( 'G', 0, 0, one, aaqq, m,
1055 $ 1, a( 1, q ), lda, ierr )
1056 sva( q ) = aaqq*sqrt( max( zero,
1057 $ one-aapq*aapq ) )
1058 mxsinj = max( mxsinj, sfmin )
1059 END IF
1060* END IF ROTOK THEN ... ELSE
1061*
1062* In the case of cancellation in updating SVA(q), SVA(p)
1063* recompute SVA(q), SVA(p).
1064*
1065 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
1066 $ THEN
1067 IF( ( aaqq.LT.rootbig ) .AND.
1068 $ ( aaqq.GT.rootsfmin ) ) THEN
1069 sva( q ) = snrm2( m, a( 1, q ), 1 )*
1070 $ work( q )
1071 ELSE
1072 t = zero
1073 aaqq = one
1074 CALL slassq( m, a( 1, q ), 1, t,
1075 $ aaqq )
1076 sva( q ) = t*sqrt( aaqq )*work( q )
1077 END IF
1078 END IF
1079 IF( ( aapp / aapp0 ).LE.rooteps ) THEN
1080 IF( ( aapp.LT.rootbig ) .AND.
1081 $ ( aapp.GT.rootsfmin ) ) THEN
1082 aapp = snrm2( m, a( 1, p ), 1 )*
1083 $ work( p )
1084 ELSE
1085 t = zero
1086 aapp = one
1087 CALL slassq( m, a( 1, p ), 1, t,
1088 $ aapp )
1089 aapp = t*sqrt( aapp )*work( p )
1090 END IF
1091 sva( p ) = aapp
1092 END IF
1093*
1094 ELSE
1095* A(:,p) and A(:,q) already numerically orthogonal
1096 IF( ir1.EQ.0 )notrot = notrot + 1
1097*[RTD] SKIPPED = SKIPPED + 1
1098 pskipped = pskipped + 1
1099 END IF
1100 ELSE
1101* A(:,q) is zero column
1102 IF( ir1.EQ.0 )notrot = notrot + 1
1103 pskipped = pskipped + 1
1104 END IF
1105*
1106 IF( ( i.LE.swband ) .AND.
1107 $ ( pskipped.GT.rowskip ) ) THEN
1108 IF( ir1.EQ.0 )aapp = -aapp
1109 notrot = 0
1110 GO TO 2103
1111 END IF
1112*
1113 2002 CONTINUE
1114* END q-LOOP
1115*
1116 2103 CONTINUE
1117* bailed out of q-loop
1118*
1119 sva( p ) = aapp
1120*
1121 ELSE
1122 sva( p ) = aapp
1123 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
1124 $ notrot = notrot + min( igl+kbl-1, n ) - p
1125 END IF
1126*
1127 2001 CONTINUE
1128* end of the p-loop
1129* end of doing the block ( ibr, ibr )
1130 1002 CONTINUE
1131* end of ir1-loop
1132*
1133* ... go to the off diagonal blocks
1134*
1135 igl = ( ibr-1 )*kbl + 1
1136*
1137 DO 2010 jbc = ibr + 1, nbl
1138*
1139 jgl = ( jbc-1 )*kbl + 1
1140*
1141* doing the block at ( ibr, jbc )
1142*
1143 ijblsk = 0
1144 DO 2100 p = igl, min( igl+kbl-1, n )
1145*
1146 aapp = sva( p )
1147 IF( aapp.GT.zero ) THEN
1148*
1149 pskipped = 0
1150*
1151 DO 2200 q = jgl, min( jgl+kbl-1, n )
1152*
1153 aaqq = sva( q )
1154 IF( aaqq.GT.zero ) THEN
1155 aapp0 = aapp
1156*
1157* .. M x 2 Jacobi SVD ..
1158*
1159* Safe Gram matrix computation
1160*
1161 IF( aaqq.GE.one ) THEN
1162 IF( aapp.GE.aaqq ) THEN
1163 rotok = ( small*aapp ).LE.aaqq
1164 ELSE
1165 rotok = ( small*aaqq ).LE.aapp
1166 END IF
1167 IF( aapp.LT.( big / aaqq ) ) THEN
1168 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
1169 $ q ), 1 )*work( p )*work( q ) /
1170 $ aaqq ) / aapp
1171 ELSE
1172 CALL scopy( m, a( 1, p ), 1,
1173 $ work( n+1 ), 1 )
1174 CALL slascl( 'G', 0, 0, aapp,
1175 $ work( p ), m, 1,
1176 $ work( n+1 ), lda, ierr )
1177 aapq = sdot( m, work( n+1 ), 1,
1178 $ a( 1, q ), 1 )*work( q ) / aaqq
1179 END IF
1180 ELSE
1181 IF( aapp.GE.aaqq ) THEN
1182 rotok = aapp.LE.( aaqq / small )
1183 ELSE
1184 rotok = aaqq.LE.( aapp / small )
1185 END IF
1186 IF( aapp.GT.( small / aaqq ) ) THEN
1187 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
1188 $ q ), 1 )*work( p )*work( q ) /
1189 $ aaqq ) / aapp
1190 ELSE
1191 CALL scopy( m, a( 1, q ), 1,
1192 $ work( n+1 ), 1 )
1193 CALL slascl( 'G', 0, 0, aaqq,
1194 $ work( q ), m, 1,
1195 $ work( n+1 ), lda, ierr )
1196 aapq = sdot( m, work( n+1 ), 1,
1197 $ a( 1, p ), 1 )*work( p ) / aapp
1198 END IF
1199 END IF
1200*
1201 mxaapq = max( mxaapq, abs( aapq ) )
1202*
1203* TO rotate or NOT to rotate, THAT is the question ...
1204*
1205 IF( abs( aapq ).GT.tol ) THEN
1206 notrot = 0
1207*[RTD] ROTATED = ROTATED + 1
1208 pskipped = 0
1209 iswrot = iswrot + 1
1210*
1211 IF( rotok ) THEN
1212*
1213 aqoap = aaqq / aapp
1214 apoaq = aapp / aaqq
1215 theta = -half*abs( aqoap-apoaq ) / aapq
1216 IF( aaqq.GT.aapp0 )theta = -theta
1217*
1218 IF( abs( theta ).GT.bigtheta ) THEN
1219 t = half / theta
1220 fastr( 3 ) = t*work( p ) / work( q )
1221 fastr( 4 ) = -t*work( q ) /
1222 $ work( p )
1223 CALL srotm( m, a( 1, p ), 1,
1224 $ a( 1, q ), 1, fastr )
1225 IF( rsvec )CALL srotm( mvl,
1226 $ v( 1, p ), 1,
1227 $ v( 1, q ), 1,
1228 $ fastr )
1229 sva( q ) = aaqq*sqrt( max( zero,
1230 $ one+t*apoaq*aapq ) )
1231 aapp = aapp*sqrt( max( zero,
1232 $ one-t*aqoap*aapq ) )
1233 mxsinj = max( mxsinj, abs( t ) )
1234 ELSE
1235*
1236* .. choose correct signum for THETA and rotate
1237*
1238 thsign = -sign( one, aapq )
1239 IF( aaqq.GT.aapp0 )thsign = -thsign
1240 t = one / ( theta+thsign*
1241 $ sqrt( one+theta*theta ) )
1242 cs = sqrt( one / ( one+t*t ) )
1243 sn = t*cs
1244 mxsinj = max( mxsinj, abs( sn ) )
1245 sva( q ) = aaqq*sqrt( max( zero,
1246 $ one+t*apoaq*aapq ) )
1247 aapp = aapp*sqrt( max( zero,
1248 $ one-t*aqoap*aapq ) )
1249*
1250 apoaq = work( p ) / work( q )
1251 aqoap = work( q ) / work( p )
1252 IF( work( p ).GE.one ) THEN
1253*
1254 IF( work( q ).GE.one ) THEN
1255 fastr( 3 ) = t*apoaq
1256 fastr( 4 ) = -t*aqoap
1257 work( p ) = work( p )*cs
1258 work( q ) = work( q )*cs
1259 CALL srotm( m, a( 1, p ), 1,
1260 $ a( 1, q ), 1,
1261 $ fastr )
1262 IF( rsvec )CALL srotm( mvl,
1263 $ v( 1, p ), 1, v( 1, q ),
1264 $ 1, fastr )
1265 ELSE
1266 CALL saxpy( m, -t*aqoap,
1267 $ a( 1, q ), 1,
1268 $ a( 1, p ), 1 )
1269 CALL saxpy( m, cs*sn*apoaq,
1270 $ a( 1, p ), 1,
1271 $ a( 1, q ), 1 )
1272 IF( rsvec ) THEN
1273 CALL saxpy( mvl, -t*aqoap,
1274 $ v( 1, q ), 1,
1275 $ v( 1, p ), 1 )
1276 CALL saxpy( mvl,
1277 $ cs*sn*apoaq,
1278 $ v( 1, p ), 1,
1279 $ v( 1, q ), 1 )
1280 END IF
1281 work( p ) = work( p )*cs
1282 work( q ) = work( q ) / cs
1283 END IF
1284 ELSE
1285 IF( work( q ).GE.one ) THEN
1286 CALL saxpy( m, t*apoaq,
1287 $ a( 1, p ), 1,
1288 $ a( 1, q ), 1 )
1289 CALL saxpy( m, -cs*sn*aqoap,
1290 $ a( 1, q ), 1,
1291 $ a( 1, p ), 1 )
1292 IF( rsvec ) THEN
1293 CALL saxpy( mvl, t*apoaq,
1294 $ v( 1, p ), 1,
1295 $ v( 1, q ), 1 )
1296 CALL saxpy( mvl,
1297 $ -cs*sn*aqoap,
1298 $ v( 1, q ), 1,
1299 $ v( 1, p ), 1 )
1300 END IF
1301 work( p ) = work( p ) / cs
1302 work( q ) = work( q )*cs
1303 ELSE
1304 IF( work( p ).GE.work( q ) )
1305 $ THEN
1306 CALL saxpy( m, -t*aqoap,
1307 $ a( 1, q ), 1,
1308 $ a( 1, p ), 1 )
1309 CALL saxpy( m, cs*sn*apoaq,
1310 $ a( 1, p ), 1,
1311 $ a( 1, q ), 1 )
1312 work( p ) = work( p )*cs
1313 work( q ) = work( q ) / cs
1314 IF( rsvec ) THEN
1315 CALL saxpy( mvl,
1316 $ -t*aqoap,
1317 $ v( 1, q ), 1,
1318 $ v( 1, p ), 1 )
1319 CALL saxpy( mvl,
1320 $ cs*sn*apoaq,
1321 $ v( 1, p ), 1,
1322 $ v( 1, q ), 1 )
1323 END IF
1324 ELSE
1325 CALL saxpy( m, t*apoaq,
1326 $ a( 1, p ), 1,
1327 $ a( 1, q ), 1 )
1328 CALL saxpy( m,
1329 $ -cs*sn*aqoap,
1330 $ a( 1, q ), 1,
1331 $ a( 1, p ), 1 )
1332 work( p ) = work( p ) / cs
1333 work( q ) = work( q )*cs
1334 IF( rsvec ) THEN
1335 CALL saxpy( mvl,
1336 $ t*apoaq, v( 1, p ),
1337 $ 1, v( 1, q ), 1 )
1338 CALL saxpy( mvl,
1339 $ -cs*sn*aqoap,
1340 $ v( 1, q ), 1,
1341 $ v( 1, p ), 1 )
1342 END IF
1343 END IF
1344 END IF
1345 END IF
1346 END IF
1347*
1348 ELSE
1349 IF( aapp.GT.aaqq ) THEN
1350 CALL scopy( m, a( 1, p ), 1,
1351 $ work( n+1 ), 1 )
1352 CALL slascl( 'G', 0, 0, aapp, one,
1353 $ m, 1, work( n+1 ), lda,
1354 $ ierr )
1355 CALL slascl( 'G', 0, 0, aaqq, one,
1356 $ m, 1, a( 1, q ), lda,
1357 $ ierr )
1358 temp1 = -aapq*work( p ) / work( q )
1359 CALL saxpy( m, temp1, work( n+1 ),
1360 $ 1, a( 1, q ), 1 )
1361 CALL slascl( 'G', 0, 0, one, aaqq,
1362 $ m, 1, a( 1, q ), lda,
1363 $ ierr )
1364 sva( q ) = aaqq*sqrt( max( zero,
1365 $ one-aapq*aapq ) )
1366 mxsinj = max( mxsinj, sfmin )
1367 ELSE
1368 CALL scopy( m, a( 1, q ), 1,
1369 $ work( n+1 ), 1 )
1370 CALL slascl( 'G', 0, 0, aaqq, one,
1371 $ m, 1, work( n+1 ), lda,
1372 $ ierr )
1373 CALL slascl( 'G', 0, 0, aapp, one,
1374 $ m, 1, a( 1, p ), lda,
1375 $ ierr )
1376 temp1 = -aapq*work( q ) / work( p )
1377 CALL saxpy( m, temp1, work( n+1 ),
1378 $ 1, a( 1, p ), 1 )
1379 CALL slascl( 'G', 0, 0, one, aapp,
1380 $ m, 1, a( 1, p ), lda,
1381 $ ierr )
1382 sva( p ) = aapp*sqrt( max( zero,
1383 $ one-aapq*aapq ) )
1384 mxsinj = max( mxsinj, sfmin )
1385 END IF
1386 END IF
1387* END IF ROTOK THEN ... ELSE
1388*
1389* In the case of cancellation in updating SVA(q)
1390* .. recompute SVA(q)
1391 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
1392 $ THEN
1393 IF( ( aaqq.LT.rootbig ) .AND.
1394 $ ( aaqq.GT.rootsfmin ) ) THEN
1395 sva( q ) = snrm2( m, a( 1, q ), 1 )*
1396 $ work( q )
1397 ELSE
1398 t = zero
1399 aaqq = one
1400 CALL slassq( m, a( 1, q ), 1, t,
1401 $ aaqq )
1402 sva( q ) = t*sqrt( aaqq )*work( q )
1403 END IF
1404 END IF
1405 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
1406 IF( ( aapp.LT.rootbig ) .AND.
1407 $ ( aapp.GT.rootsfmin ) ) THEN
1408 aapp = snrm2( m, a( 1, p ), 1 )*
1409 $ work( p )
1410 ELSE
1411 t = zero
1412 aapp = one
1413 CALL slassq( m, a( 1, p ), 1, t,
1414 $ aapp )
1415 aapp = t*sqrt( aapp )*work( p )
1416 END IF
1417 sva( p ) = aapp
1418 END IF
1419* end of OK rotation
1420 ELSE
1421 notrot = notrot + 1
1422*[RTD] SKIPPED = SKIPPED + 1
1423 pskipped = pskipped + 1
1424 ijblsk = ijblsk + 1
1425 END IF
1426 ELSE
1427 notrot = notrot + 1
1428 pskipped = pskipped + 1
1429 ijblsk = ijblsk + 1
1430 END IF
1431*
1432 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
1433 $ THEN
1434 sva( p ) = aapp
1435 notrot = 0
1436 GO TO 2011
1437 END IF
1438 IF( ( i.LE.swband ) .AND.
1439 $ ( pskipped.GT.rowskip ) ) THEN
1440 aapp = -aapp
1441 notrot = 0
1442 GO TO 2203
1443 END IF
1444*
1445 2200 CONTINUE
1446* end of the q-loop
1447 2203 CONTINUE
1448*
1449 sva( p ) = aapp
1450*
1451 ELSE
1452*
1453 IF( aapp.EQ.zero )notrot = notrot +
1454 $ min( jgl+kbl-1, n ) - jgl + 1
1455 IF( aapp.LT.zero )notrot = 0
1456*
1457 END IF
1458*
1459 2100 CONTINUE
1460* end of the p-loop
1461 2010 CONTINUE
1462* end of the jbc-loop
1463 2011 CONTINUE
1464*2011 bailed out of the jbc-loop
1465 DO 2012 p = igl, min( igl+kbl-1, n )
1466 sva( p ) = abs( sva( p ) )
1467 2012 CONTINUE
1468***
1469 2000 CONTINUE
1470*2000 :: end of the ibr-loop
1471*
1472* .. update SVA(N)
1473 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1474 $ THEN
1475 sva( n ) = snrm2( m, a( 1, n ), 1 )*work( n )
1476 ELSE
1477 t = zero
1478 aapp = one
1479 CALL slassq( m, a( 1, n ), 1, t, aapp )
1480 sva( n ) = t*sqrt( aapp )*work( n )
1481 END IF
1482*
1483* Additional steering devices
1484*
1485 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
1486 $ ( iswrot.LE.n ) ) )swband = i
1487*
1488 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( float( n ) )*
1489 $ tol ) .AND. ( float( n )*mxaapq*mxsinj.LT.tol ) ) THEN
1490 GO TO 1994
1491 END IF
1492*
1493 IF( notrot.GE.emptsw )GO TO 1994
1494*
1495 1993 CONTINUE
1496* end i=1:NSWEEP loop
1497*
1498* #:( Reaching this point means that the procedure has not converged.
1499 info = nsweep - 1
1500 GO TO 1995
1501*
1502 1994 CONTINUE
1503* #:) Reaching this point means numerical convergence after the i-th
1504* sweep.
1505*
1506 info = 0
1507* #:) INFO = 0 confirms successful iterations.
1508 1995 CONTINUE
1509*
1510* Sort the singular values and find how many are above
1511* the underflow threshold.
1512*
1513 n2 = 0
1514 n4 = 0
1515 DO 5991 p = 1, n - 1
1516 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
1517 IF( p.NE.q ) THEN
1518 temp1 = sva( p )
1519 sva( p ) = sva( q )
1520 sva( q ) = temp1
1521 temp1 = work( p )
1522 work( p ) = work( q )
1523 work( q ) = temp1
1524 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1525 IF( rsvec )CALL sswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
1526 END IF
1527 IF( sva( p ).NE.zero ) THEN
1528 n4 = n4 + 1
1529 IF( sva( p )*skl.GT.sfmin )n2 = n2 + 1
1530 END IF
1531 5991 CONTINUE
1532 IF( sva( n ).NE.zero ) THEN
1533 n4 = n4 + 1
1534 IF( sva( n )*skl.GT.sfmin )n2 = n2 + 1
1535 END IF
1536*
1537* Normalize the left singular vectors.
1538*
1539 IF( lsvec .OR. uctol ) THEN
1540 DO 1998 p = 1, n2
1541 CALL sscal( m, work( p ) / sva( p ), a( 1, p ), 1 )
1542 1998 CONTINUE
1543 END IF
1544*
1545* Scale the product of Jacobi rotations (assemble the fast rotations).
1546*
1547 IF( rsvec ) THEN
1548 IF( applv ) THEN
1549 DO 2398 p = 1, n
1550 CALL sscal( mvl, work( p ), v( 1, p ), 1 )
1551 2398 CONTINUE
1552 ELSE
1553 DO 2399 p = 1, n
1554 temp1 = one / snrm2( mvl, v( 1, p ), 1 )
1555 CALL sscal( mvl, temp1, v( 1, p ), 1 )
1556 2399 CONTINUE
1557 END IF
1558 END IF
1559*
1560* Undo scaling, if necessary (and possible).
1561 IF( ( ( skl.GT.one ) .AND. ( sva( 1 ).LT.( big / skl ) ) )
1562 $ .OR. ( ( skl.LT.one ) .AND. ( sva( max( n2, 1 ) ) .GT.
1563 $ ( sfmin / skl ) ) ) ) THEN
1564 DO 2400 p = 1, n
1565 sva( p ) = skl*sva( p )
1566 2400 CONTINUE
1567 skl = one
1568 END IF
1569*
1570 work( 1 ) = skl
1571* The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
1572* then some of the singular values may overflow or underflow and
1573* the spectrum is given in this factored representation.
1574*
1575 work( 2 ) = float( n4 )
1576* N4 is the number of computed nonzero singular values of A.
1577*
1578 work( 3 ) = float( n2 )
1579* N2 is the number of singular values of A greater than SFMIN.
1580* If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
1581* that may carry some information.
1582*
1583 work( 4 ) = float( i )
1584* i is the index of the last sweep before declaring convergence.
1585*
1586 work( 5 ) = mxaapq
1587* MXAAPQ is the largest absolute value of scaled pivots in the
1588* last sweep
1589*
1590 work( 6 ) = mxsinj
1591* MXSINJ is the largest absolute value of the sines of Jacobi angles
1592* in the last sweep
1593*
1594 RETURN
1595* ..
1596* .. END OF SGESVJ
1597* ..
1598 END
slassq
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f90:137
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
sgesvj
subroutine sgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
SGESVJ
Definition: sgesvj.f:323
sgsvj1
subroutine sgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
SGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivot...
Definition: sgsvj1.f:236
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:218
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
srotm
subroutine srotm(N, SX, INCX, SY, INCY, SPARAM)
SROTM
Definition: srotm.f:97
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89