LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cgsvj1.f
Go to the documentation of this file.
1 *> \brief \b CGSVJ1 pre-processor for the routine cgesvj, 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 CGSVJ1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgsvj1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgsvj1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgsvj1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
22 * EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * REAL EPS, SFMIN, TOL
26 * INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
27 * CHARACTER*1 JOBV
28 * ..
29 * .. Array Arguments ..
30 * COMPLEX A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
31 * REAL SVA( N )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CGSVJ1 is called from CGESVJ as a pre-processor and that is its main
41 *> purpose. It applies Jacobi rotations in the same way as CGESVJ 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 *> CGSVJ1 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 COMPLEX 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 COMPLEX 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 REAL 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 COMPLEX 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 REAL
176 *> EPS = SLAMCH('Epsilon')
177 *> \endverbatim
178 *>
179 *> \param[in] SFMIN
180 *> \verbatim
181 *> SFMIN is REAL
182 *> SFMIN = SLAMCH('Safe Minimum')
183 *> \endverbatim
184 *>
185 *> \param[in] TOL
186 *> \verbatim
187 *> TOL is REAL
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 ABS(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 COMPLEX 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 complexOTHERcomputational
227 *
228 *> \par Contributor:
229 * ==================
230 *>
231 *> Zlatko Drmac (Zagreb, Croatia)
232 *
233 * =====================================================================
234  SUBROUTINE cgsvj1( 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  REAL EPS, SFMIN, TOL
243  INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
244  CHARACTER*1 JOBV
245 * ..
246 * .. Array Arguments ..
247  COMPLEX A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
248  REAL SVA( N )
249 * ..
250 *
251 * =====================================================================
252 *
253 * .. Local Parameters ..
254  REAL ZERO, HALF, ONE
255  parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
256 * ..
257 * .. Local Scalars ..
258  COMPLEX AAPQ, OMPQ
259  REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
260  $ bigtheta, cs, mxaapq, mxsinj, rootbig,
261  $ rooteps, rootsfmin, roottol, small, sn, t,
262  $ temp1, theta, thsign
263  INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
264  $ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
265  $ p, pskipped, q, rowskip, swband
266  LOGICAL APPLV, ROTOK, RSVEC
267 * ..
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC abs, max, conjg, real, min, sign, sqrt
271 * ..
272 * .. External Functions ..
273  REAL SCNRM2
274  COMPLEX CDOTC
275  INTEGER ISAMAX
276  LOGICAL LSAME
277  EXTERNAL isamax, lsame, cdotc, scnrm2
278 * ..
279 * .. External Subroutines ..
280 * .. from BLAS
281  EXTERNAL ccopy, crot, cswap, caxpy
282 * .. from LAPACK
283  EXTERNAL clascl, classq, xerbla
284 * ..
285 * .. Executable Statements ..
286 *
287 * Test the input parameters.
288 *
289  applv = lsame( jobv, 'A' )
290  rsvec = lsame( jobv, 'V' )
291  IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
292  info = -1
293  ELSE IF( m.LT.0 ) THEN
294  info = -2
295  ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
296  info = -3
297  ELSE IF( n1.LT.0 ) THEN
298  info = -4
299  ELSE IF( lda.LT.m ) THEN
300  info = -6
301  ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
302  info = -9
303  ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
304  $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
305  info = -11
306  ELSE IF( tol.LE.eps ) THEN
307  info = -14
308  ELSE IF( nsweep.LT.0 ) THEN
309  info = -15
310  ELSE IF( lwork.LT.m ) THEN
311  info = -17
312  ELSE
313  info = 0
314  END IF
315 *
316 * #:(
317  IF( info.NE.0 ) THEN
318  CALL xerbla( 'CGSVJ1', -info )
319  RETURN
320  END IF
321 *
322  IF( rsvec ) THEN
323  mvl = n
324  ELSE IF( applv ) THEN
325  mvl = mv
326  END IF
327  rsvec = rsvec .OR. applv
328 
329  rooteps = sqrt( eps )
330  rootsfmin = sqrt( sfmin )
331  small = sfmin / eps
332  big = one / sfmin
333  rootbig = one / rootsfmin
334 * LARGE = BIG / SQRT( REAL( M*N ) )
335  bigtheta = one / rooteps
336  roottol = sqrt( tol )
337 *
338 * .. Initialize the right singular vector matrix ..
339 *
340 * RSVEC = LSAME( JOBV, 'Y' )
341 *
342  emptsw = n1*( n-n1 )
343  notrot = 0
344 *
345 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
346 *
347  kbl = min( 8, n )
348  nblr = n1 / kbl
349  IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
350 
351 * .. the tiling is nblr-by-nblc [tiles]
352 
353  nblc = ( n-n1 ) / kbl
354  IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
355  blskip = ( kbl**2 ) + 1
356 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
357 
358  rowskip = min( 5, kbl )
359 *[TP] ROWSKIP is a tuning parameter.
360  swband = 0
361 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
362 * if CGESVJ is used as a computational routine in the preconditioned
363 * Jacobi SVD algorithm CGEJSV.
364 *
365 *
366 * | * * * [x] [x] [x]|
367 * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
368 * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
369 * |[x] [x] [x] * * * |
370 * |[x] [x] [x] * * * |
371 * |[x] [x] [x] * * * |
372 *
373 *
374  DO 1993 i = 1, nsweep
375 *
376 * .. go go go ...
377 *
378  mxaapq = zero
379  mxsinj = zero
380  iswrot = 0
381 *
382  notrot = 0
383  pskipped = 0
384 *
385 * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
386 * 1 <= p < q <= N. This is the first step toward a blocked implementation
387 * of the rotations. New implementation, based on block transformations,
388 * is under development.
389 *
390  DO 2000 ibr = 1, nblr
391 *
392  igl = ( ibr-1 )*kbl + 1
393 *
394 
395 *
396 * ... go to the off diagonal blocks
397 *
398  igl = ( ibr-1 )*kbl + 1
399 *
400 * DO 2010 jbc = ibr + 1, NBL
401  DO 2010 jbc = 1, nblc
402 *
403  jgl = ( jbc-1 )*kbl + n1 + 1
404 *
405 * doing the block at ( ibr, jbc )
406 *
407  ijblsk = 0
408  DO 2100 p = igl, min( igl+kbl-1, n1 )
409 *
410  aapp = sva( p )
411  IF( aapp.GT.zero ) THEN
412 *
413  pskipped = 0
414 *
415  DO 2200 q = jgl, min( jgl+kbl-1, n )
416 *
417  aaqq = sva( q )
418  IF( aaqq.GT.zero ) THEN
419  aapp0 = aapp
420 *
421 * .. M x 2 Jacobi SVD ..
422 *
423 * Safe Gram matrix computation
424 *
425  IF( aaqq.GE.one ) THEN
426  IF( aapp.GE.aaqq ) THEN
427  rotok = ( small*aapp ).LE.aaqq
428  ELSE
429  rotok = ( small*aaqq ).LE.aapp
430  END IF
431  IF( aapp.LT.( big / aaqq ) ) THEN
432  aapq = ( cdotc( m, a( 1, p ), 1,
433  $ a( 1, q ), 1 ) / aaqq ) / aapp
434  ELSE
435  CALL ccopy( m, a( 1, p ), 1,
436  $ work, 1 )
437  CALL clascl( 'G', 0, 0, aapp,
438  $ one, m, 1,
439  $ work, lda, ierr )
440  aapq = cdotc( m, work, 1,
441  $ a( 1, q ), 1 ) / aaqq
442  END IF
443  ELSE
444  IF( aapp.GE.aaqq ) THEN
445  rotok = aapp.LE.( aaqq / small )
446  ELSE
447  rotok = aaqq.LE.( aapp / small )
448  END IF
449  IF( aapp.GT.( small / aaqq ) ) THEN
450  aapq = ( cdotc( m, a( 1, p ), 1,
451  $ a( 1, q ), 1 ) / max(aaqq,aapp) )
452  $ / min(aaqq,aapp)
453  ELSE
454  CALL ccopy( m, a( 1, q ), 1,
455  $ work, 1 )
456  CALL clascl( 'G', 0, 0, aaqq,
457  $ one, m, 1,
458  $ work, lda, ierr )
459  aapq = cdotc( m, a( 1, p ), 1,
460  $ work, 1 ) / aapp
461  END IF
462  END IF
463 *
464 * AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
465  aapq1 = -abs(aapq)
466  mxaapq = max( mxaapq, -aapq1 )
467 *
468 * TO rotate or NOT to rotate, THAT is the question ...
469 *
470  IF( abs( aapq1 ).GT.tol ) THEN
471  ompq = aapq / abs(aapq)
472  notrot = 0
473 *[RTD] ROTATED = ROTATED + 1
474  pskipped = 0
475  iswrot = iswrot + 1
476 *
477  IF( rotok ) THEN
478 *
479  aqoap = aaqq / aapp
480  apoaq = aapp / aaqq
481  theta = -half*abs( aqoap-apoaq )/ aapq1
482  IF( aaqq.GT.aapp0 )theta = -theta
483 *
484  IF( abs( theta ).GT.bigtheta ) THEN
485  t = half / theta
486  cs = one
487  CALL crot( m, a(1,p), 1, a(1,q), 1,
488  $ cs, conjg(ompq)*t )
489  IF( rsvec ) THEN
490  CALL crot( mvl, v(1,p), 1,
491  $ v(1,q), 1, cs, conjg(ompq)*t )
492  END IF
493  sva( q ) = aaqq*sqrt( max( zero,
494  $ one+t*apoaq*aapq1 ) )
495  aapp = aapp*sqrt( max( zero,
496  $ one-t*aqoap*aapq1 ) )
497  mxsinj = max( mxsinj, abs( t ) )
498  ELSE
499 *
500 * .. choose correct signum for THETA and rotate
501 *
502  thsign = -sign( one, aapq1 )
503  IF( aaqq.GT.aapp0 )thsign = -thsign
504  t = one / ( theta+thsign*
505  $ sqrt( one+theta*theta ) )
506  cs = sqrt( one / ( one+t*t ) )
507  sn = t*cs
508  mxsinj = max( mxsinj, abs( sn ) )
509  sva( q ) = aaqq*sqrt( max( zero,
510  $ one+t*apoaq*aapq1 ) )
511  aapp = aapp*sqrt( max( zero,
512  $ one-t*aqoap*aapq1 ) )
513 *
514  CALL crot( m, a(1,p), 1, a(1,q), 1,
515  $ cs, conjg(ompq)*sn )
516  IF( rsvec ) THEN
517  CALL crot( mvl, v(1,p), 1,
518  $ v(1,q), 1, cs, conjg(ompq)*sn )
519  END IF
520  END IF
521  d(p) = -d(q) * ompq
522 *
523  ELSE
524 * .. have to use modified Gram-Schmidt like transformation
525  IF( aapp.GT.aaqq ) THEN
526  CALL ccopy( m, a( 1, p ), 1,
527  $ work, 1 )
528  CALL clascl( 'G', 0, 0, aapp, one,
529  $ m, 1, work,lda,
530  $ ierr )
531  CALL clascl( 'G', 0, 0, aaqq, one,
532  $ m, 1, a( 1, q ), lda,
533  $ ierr )
534  CALL caxpy( m, -aapq, work,
535  $ 1, a( 1, q ), 1 )
536  CALL clascl( 'G', 0, 0, one, aaqq,
537  $ m, 1, a( 1, q ), lda,
538  $ ierr )
539  sva( q ) = aaqq*sqrt( max( zero,
540  $ one-aapq1*aapq1 ) )
541  mxsinj = max( mxsinj, sfmin )
542  ELSE
543  CALL ccopy( m, a( 1, q ), 1,
544  $ work, 1 )
545  CALL clascl( 'G', 0, 0, aaqq, one,
546  $ m, 1, work,lda,
547  $ ierr )
548  CALL clascl( 'G', 0, 0, aapp, one,
549  $ m, 1, a( 1, p ), lda,
550  $ ierr )
551  CALL caxpy( m, -conjg(aapq),
552  $ work, 1, a( 1, p ), 1 )
553  CALL clascl( 'G', 0, 0, one, aapp,
554  $ m, 1, a( 1, p ), lda,
555  $ ierr )
556  sva( p ) = aapp*sqrt( max( zero,
557  $ one-aapq1*aapq1 ) )
558  mxsinj = max( mxsinj, sfmin )
559  END IF
560  END IF
561 * END IF ROTOK THEN ... ELSE
562 *
563 * In the case of cancellation in updating SVA(q), SVA(p)
564 * .. recompute SVA(q), SVA(p)
565  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
566  $ THEN
567  IF( ( aaqq.LT.rootbig ) .AND.
568  $ ( aaqq.GT.rootsfmin ) ) THEN
569  sva( q ) = scnrm2( m, a( 1, q ), 1)
570  ELSE
571  t = zero
572  aaqq = one
573  CALL classq( m, a( 1, q ), 1, t,
574  $ aaqq )
575  sva( q ) = t*sqrt( aaqq )
576  END IF
577  END IF
578  IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
579  IF( ( aapp.LT.rootbig ) .AND.
580  $ ( aapp.GT.rootsfmin ) ) THEN
581  aapp = scnrm2( m, a( 1, p ), 1 )
582  ELSE
583  t = zero
584  aapp = one
585  CALL classq( m, a( 1, p ), 1, t,
586  $ aapp )
587  aapp = t*sqrt( aapp )
588  END IF
589  sva( p ) = aapp
590  END IF
591 * end of OK rotation
592  ELSE
593  notrot = notrot + 1
594 *[RTD] SKIPPED = SKIPPED + 1
595  pskipped = pskipped + 1
596  ijblsk = ijblsk + 1
597  END IF
598  ELSE
599  notrot = notrot + 1
600  pskipped = pskipped + 1
601  ijblsk = ijblsk + 1
602  END IF
603 *
604  IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
605  $ THEN
606  sva( p ) = aapp
607  notrot = 0
608  GO TO 2011
609  END IF
610  IF( ( i.LE.swband ) .AND.
611  $ ( pskipped.GT.rowskip ) ) THEN
612  aapp = -aapp
613  notrot = 0
614  GO TO 2203
615  END IF
616 *
617  2200 CONTINUE
618 * end of the q-loop
619  2203 CONTINUE
620 *
621  sva( p ) = aapp
622 *
623  ELSE
624 *
625  IF( aapp.EQ.zero )notrot = notrot +
626  $ min( jgl+kbl-1, n ) - jgl + 1
627  IF( aapp.LT.zero )notrot = 0
628 *
629  END IF
630 *
631  2100 CONTINUE
632 * end of the p-loop
633  2010 CONTINUE
634 * end of the jbc-loop
635  2011 CONTINUE
636 *2011 bailed out of the jbc-loop
637  DO 2012 p = igl, min( igl+kbl-1, n )
638  sva( p ) = abs( sva( p ) )
639  2012 CONTINUE
640 ***
641  2000 CONTINUE
642 *2000 :: end of the ibr-loop
643 *
644 * .. update SVA(N)
645  IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
646  $ THEN
647  sva( n ) = scnrm2( m, a( 1, n ), 1 )
648  ELSE
649  t = zero
650  aapp = one
651  CALL classq( m, a( 1, n ), 1, t, aapp )
652  sva( n ) = t*sqrt( aapp )
653  END IF
654 *
655 * Additional steering devices
656 *
657  IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
658  $ ( iswrot.LE.n ) ) )swband = i
659 *
660  IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( real( n ) )*
661  $ tol ) .AND. ( real( n )*mxaapq*mxsinj.LT.tol ) ) THEN
662  GO TO 1994
663  END IF
664 *
665  IF( notrot.GE.emptsw )GO TO 1994
666 *
667  1993 CONTINUE
668 * end i=1:NSWEEP loop
669 *
670 * #:( Reaching this point means that the procedure has not converged.
671  info = nsweep - 1
672  GO TO 1995
673 *
674  1994 CONTINUE
675 * #:) Reaching this point means numerical convergence after the i-th
676 * sweep.
677 *
678  info = 0
679 * #:) INFO = 0 confirms successful iterations.
680  1995 CONTINUE
681 *
682 * Sort the vector SVA() of column norms.
683  DO 5991 p = 1, n - 1
684  q = isamax( n-p+1, sva( p ), 1 ) + p - 1
685  IF( p.NE.q ) THEN
686  temp1 = sva( p )
687  sva( p ) = sva( q )
688  sva( q ) = temp1
689  aapq = d( p )
690  d( p ) = d( q )
691  d( q ) = aapq
692  CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
693  IF( rsvec )CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
694  END IF
695  5991 CONTINUE
696 *
697 *
698  RETURN
699 * ..
700 * .. END OF CGSVJ1
701 * ..
702  END
subroutine classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f90:126
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: crot.f:103
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:143
subroutine cgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
CGSVJ1 pre-processor for the routine cgesvj, applies Jacobi rotations targeting only particular pivot...
Definition: cgsvj1.f:236