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