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