LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
cgsvj0.f
1 *> \brief \b CGSVJ0 pre-processor for the routine cgesvj.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGSVJ0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgsvj0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgsvj0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgsvj0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGSVJ0( 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 * REAL EPS, SFMIN, TOL
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 *> CGSVJ0 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 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 COMPLEX 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 * diag(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 COMPLEX array, dimension (N)
98 *> The array D accumulates the scaling factors from the complex 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 REAL 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 A_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 COMPLEX 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 REAL
146 *> EPS = SLAMCH('Epsilon')
147 *> \endverbatim
148 *>
149 *> \param[in] SFMIN
150 *> \verbatim
151 *> SFMIN is REAL
152 *> SFMIN = SLAMCH('Safe Minimum')
153 *> \endverbatim
154 *>
155 *> \param[in] TOL
156 *> \verbatim
157 *> TOL is REAL
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 ABS(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 COMPLEX 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 *> CGSVJ0 is used just to enable CGESVJ to call a simplified version of
202 *> itself to work on a submatrix of the original matrix.
203 *>
204 *> \par Contributor:
205 * ==================
206 *>
207 *> Zlatko Drmac (Zagreb, Croatia)
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 cgsvj0( 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  IMPLICIT NONE
224 * .. Scalar Arguments ..
225  INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
226  REAL EPS, SFMIN, TOL
227  CHARACTER*1 JOBV
228 * ..
229 * .. Array Arguments ..
230  COMPLEX A( lda, * ), D( n ), V( ldv, * ), WORK( lwork )
231  REAL SVA( n )
232 * ..
233 *
234 * =====================================================================
235 *
236 * .. Local Parameters ..
237  REAL ZERO, HALF, ONE
238  parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
239  COMPLEX CZERO, CONE
240  parameter( czero = (0.0e0, 0.0e0), cone = (1.0e0, 0.0e0) )
241 * ..
242 * .. Local Scalars ..
243  COMPLEX AAPQ, OMPQ
244  REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
245  $ bigtheta, cs, mxaapq, mxsinj, rootbig, rooteps,
246  $ rootsfmin, roottol, small, sn, t, temp1, theta,
247  $ thsign
248  INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
249  $ iswrot, jbc, jgl, kbl, lkahead, mvl, nbl,
250  $ notrot, p, pskipped, q, rowskip, swband
251  LOGICAL APPLV, ROTOK, RSVEC
252 * ..
253 * ..
254 * .. Intrinsic Functions ..
255  INTRINSIC abs, max, conjg, REAL, MIN, SIGN, SQRT
256 * ..
257 * .. External Functions ..
258  REAL SCNRM2
259  COMPLEX CDOTC
260  INTEGER ISAMAX
261  LOGICAL LSAME
262  EXTERNAL isamax, lsame, cdotc, scnrm2
263 * ..
264 * ..
265 * .. External Subroutines ..
266 * ..
267 * from BLAS
268  EXTERNAL ccopy, crot, cswap, caxpy
269 * from LAPACK
270  EXTERNAL clascl, classq, xerbla
271 * ..
272 * .. Executable Statements ..
273 *
274 * Test the input parameters.
275 *
276  applv = lsame( jobv, 'A' )
277  rsvec = lsame( jobv, 'V' )
278  IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
279  info = -1
280  ELSE IF( m.LT.0 ) THEN
281  info = -2
282  ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
283  info = -3
284  ELSE IF( lda.LT.m ) THEN
285  info = -5
286  ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
287  info = -8
288  ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
289  $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
290  info = -10
291  ELSE IF( tol.LE.eps ) THEN
292  info = -13
293  ELSE IF( nsweep.LT.0 ) THEN
294  info = -14
295  ELSE IF( lwork.LT.m ) THEN
296  info = -16
297  ELSE
298  info = 0
299  END IF
300 *
301 * #:(
302  IF( info.NE.0 ) THEN
303  CALL xerbla( 'CGSVJ0', -info )
304  RETURN
305  END IF
306 *
307  IF( rsvec ) THEN
308  mvl = n
309  ELSE IF( applv ) THEN
310  mvl = mv
311  END IF
312  rsvec = rsvec .OR. applv
313 
314  rooteps = sqrt( eps )
315  rootsfmin = sqrt( sfmin )
316  small = sfmin / eps
317  big = one / sfmin
318  rootbig = one / rootsfmin
319  bigtheta = one / rooteps
320  roottol = sqrt( tol )
321 *
322 * .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
323 *
324  emptsw = ( n*( n-1 ) ) / 2
325  notrot = 0
326 *
327 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
328 *
329 
330  swband = 0
331 *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
332 * if CGESVJ is used as a computational routine in the preconditioned
333 * Jacobi SVD algorithm CGEJSV. For sweeps i=1:SWBAND the procedure
334 * works on pivots inside a band-like region around the diagonal.
335 * The boundaries are determined dynamically, based on the number of
336 * pivots above a threshold.
337 *
338  kbl = min( 8, n )
339 *[TP] KBL is a tuning parameter that defines the tile size in the
340 * tiling of the p-q loops of pivot pairs. In general, an optimal
341 * value of KBL depends on the matrix dimensions and on the
342 * parameters of the computer's memory.
343 *
344  nbl = n / kbl
345  IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
346 *
347  blskip = kbl**2
348 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
349 *
350  rowskip = min( 5, kbl )
351 *[TP] ROWSKIP is a tuning parameter.
352 *
353  lkahead = 1
354 *[TP] LKAHEAD is a tuning parameter.
355 *
356 * Quasi block transformations, using the lower (upper) triangular
357 * structure of the input matrix. The quasi-block-cycling usually
358 * invokes cubic convergence. Big part of this cycle is done inside
359 * canonical subspaces of dimensions less than M.
360 *
361 *
362 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
363 *
364  DO 1993 i = 1, nsweep
365 *
366 * .. go go go ...
367 *
368  mxaapq = zero
369  mxsinj = zero
370  iswrot = 0
371 *
372  notrot = 0
373  pskipped = 0
374 *
375 * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
376 * 1 <= p < q <= N. This is the first step toward a blocked implementation
377 * of the rotations. New implementation, based on block transformations,
378 * is under development.
379 *
380  DO 2000 ibr = 1, nbl
381 *
382  igl = ( ibr-1 )*kbl + 1
383 *
384  DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
385 *
386  igl = igl + ir1*kbl
387 *
388  DO 2001 p = igl, min( igl+kbl-1, n-1 )
389 *
390 * .. de Rijk's pivoting
391 *
392  q = isamax( n-p+1, sva( p ), 1 ) + p - 1
393  IF( p.NE.q ) THEN
394  CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
395  IF( rsvec )CALL cswap( mvl, v( 1, p ), 1,
396  $ v( 1, q ), 1 )
397  temp1 = sva( p )
398  sva( p ) = sva( q )
399  sva( q ) = temp1
400  aapq = d(p)
401  d(p) = d(q)
402  d(q) = aapq
403  END IF
404 *
405  IF( ir1.EQ.0 ) THEN
406 *
407 * Column norms are periodically updated by explicit
408 * norm computation.
409 * Caveat:
410 * Unfortunately, some BLAS implementations compute SNCRM2(M,A(1,p),1)
411 * as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
412 * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
413 * underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
414 * Hence, SCNRM2 cannot be trusted, not even in the case when
415 * the true norm is far from the under(over)flow boundaries.
416 * If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF
417 * below should be replaced with "AAPP = SCNRM2( M, A(1,p), 1 )".
418 *
419  IF( ( sva( p ).LT.rootbig ) .AND.
420  $ ( sva( p ).GT.rootsfmin ) ) THEN
421  sva( p ) = scnrm2( m, a( 1, p ), 1 )
422  ELSE
423  temp1 = zero
424  aapp = one
425  CALL classq( m, a( 1, p ), 1, temp1, aapp )
426  sva( p ) = temp1*sqrt( aapp )
427  END IF
428  aapp = sva( p )
429  ELSE
430  aapp = sva( p )
431  END IF
432 *
433  IF( aapp.GT.zero ) THEN
434 *
435  pskipped = 0
436 *
437  DO 2002 q = p + 1, min( igl+kbl-1, n )
438 *
439  aaqq = sva( q )
440 *
441  IF( aaqq.GT.zero ) THEN
442 *
443  aapp0 = aapp
444  IF( aaqq.GE.one ) THEN
445  rotok = ( small*aapp ).LE.aaqq
446  IF( aapp.LT.( big / aaqq ) ) THEN
447  aapq = ( cdotc( m, a( 1, p ), 1,
448  $ a( 1, q ), 1 ) / aaqq ) / aapp
449  ELSE
450  CALL ccopy( m, a( 1, p ), 1,
451  $ work, 1 )
452  CALL clascl( 'G', 0, 0, aapp, one,
453  $ m, 1, work, lda, ierr )
454  aapq = cdotc( m, work, 1,
455  $ a( 1, q ), 1 ) / aaqq
456  END IF
457  ELSE
458  rotok = aapp.LE.( aaqq / small )
459  IF( aapp.GT.( small / aaqq ) ) THEN
460  aapq = ( cdotc( m, a( 1, p ), 1,
461  $ a( 1, q ), 1 ) / aapp ) / aaqq
462  ELSE
463  CALL ccopy( m, a( 1, q ), 1,
464  $ work, 1 )
465  CALL clascl( 'G', 0, 0, aaqq,
466  $ one, m, 1,
467  $ work, lda, ierr )
468  aapq = cdotc( m, a( 1, p ), 1,
469  $ work, 1 ) / aapp
470  END IF
471  END IF
472 *
473 * AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
474  aapq1 = -abs(aapq)
475  mxaapq = max( mxaapq, -aapq1 )
476 *
477 * TO rotate or NOT to rotate, THAT is the question ...
478 *
479  IF( abs( aapq1 ).GT.tol ) THEN
480  ompq = aapq / abs(aapq)
481 *
482 * .. rotate
483 *[RTD] ROTATED = ROTATED + ONE
484 *
485  IF( ir1.EQ.0 ) THEN
486  notrot = 0
487  pskipped = 0
488  iswrot = iswrot + 1
489  END IF
490 *
491  IF( rotok ) THEN
492 *
493  aqoap = aaqq / aapp
494  apoaq = aapp / aaqq
495  theta = -half*abs( aqoap-apoaq )/aapq1
496 *
497  IF( abs( theta ).GT.bigtheta ) THEN
498 *
499  t = half / theta
500  cs = one
501 
502  CALL crot( m, a(1,p), 1, a(1,q), 1,
503  $ cs, conjg(ompq)*t )
504  IF ( rsvec ) THEN
505  CALL crot( mvl, v(1,p), 1,
506  $ v(1,q), 1, cs, conjg(ompq)*t )
507  END IF
508 
509  sva( q ) = aaqq*sqrt( max( zero,
510  $ one+t*apoaq*aapq1 ) )
511  aapp = aapp*sqrt( max( zero,
512  $ one-t*aqoap*aapq1 ) )
513  mxsinj = max( mxsinj, abs( t ) )
514 *
515  ELSE
516 *
517 * .. choose correct signum for THETA and rotate
518 *
519  thsign = -sign( one, aapq1 )
520  t = one / ( theta+thsign*
521  $ sqrt( one+theta*theta ) )
522  cs = sqrt( one / ( one+t*t ) )
523  sn = t*cs
524 *
525  mxsinj = max( mxsinj, abs( sn ) )
526  sva( q ) = aaqq*sqrt( max( zero,
527  $ one+t*apoaq*aapq1 ) )
528  aapp = aapp*sqrt( max( zero,
529  $ one-t*aqoap*aapq1 ) )
530 *
531  CALL crot( m, a(1,p), 1, a(1,q), 1,
532  $ cs, conjg(ompq)*sn )
533  IF ( rsvec ) THEN
534  CALL crot( mvl, v(1,p), 1,
535  $ v(1,q), 1, cs, conjg(ompq)*sn )
536  END IF
537  END IF
538  d(p) = -d(q) * ompq
539 *
540  ELSE
541 * .. have to use modified Gram-Schmidt like transformation
542  CALL ccopy( m, a( 1, p ), 1,
543  $ work, 1 )
544  CALL clascl( 'G', 0, 0, aapp, one, m,
545  $ 1, work, lda,
546  $ ierr )
547  CALL clascl( 'G', 0, 0, aaqq, one, m,
548  $ 1, a( 1, q ), lda, ierr )
549  CALL caxpy( m, -aapq, work, 1,
550  $ a( 1, q ), 1 )
551  CALL clascl( 'G', 0, 0, one, aaqq, m,
552  $ 1, a( 1, q ), lda, ierr )
553  sva( q ) = aaqq*sqrt( max( zero,
554  $ one-aapq1*aapq1 ) )
555  mxsinj = max( mxsinj, sfmin )
556  END IF
557 * END IF ROTOK THEN ... ELSE
558 *
559 * In the case of cancellation in updating SVA(q), SVA(p)
560 * recompute SVA(q), SVA(p).
561 *
562  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
563  $ THEN
564  IF( ( aaqq.LT.rootbig ) .AND.
565  $ ( aaqq.GT.rootsfmin ) ) THEN
566  sva( q ) = scnrm2( m, a( 1, q ), 1 )
567  ELSE
568  t = zero
569  aaqq = one
570  CALL classq( m, a( 1, q ), 1, t,
571  $ aaqq )
572  sva( q ) = t*sqrt( aaqq )
573  END IF
574  END IF
575  IF( ( aapp / aapp0 ).LE.rooteps ) THEN
576  IF( ( aapp.LT.rootbig ) .AND.
577  $ ( aapp.GT.rootsfmin ) ) THEN
578  aapp = scnrm2( m, a( 1, p ), 1 )
579  ELSE
580  t = zero
581  aapp = one
582  CALL classq( m, a( 1, p ), 1, t,
583  $ aapp )
584  aapp = t*sqrt( aapp )
585  END IF
586  sva( p ) = aapp
587  END IF
588 *
589  ELSE
590 * A(:,p) and A(:,q) already numerically orthogonal
591  IF( ir1.EQ.0 )notrot = notrot + 1
592 *[RTD] SKIPPED = SKIPPED + 1
593  pskipped = pskipped + 1
594  END IF
595  ELSE
596 * A(:,q) is zero column
597  IF( ir1.EQ.0 )notrot = notrot + 1
598  pskipped = pskipped + 1
599  END IF
600 *
601  IF( ( i.LE.swband ) .AND.
602  $ ( pskipped.GT.rowskip ) ) THEN
603  IF( ir1.EQ.0 )aapp = -aapp
604  notrot = 0
605  GO TO 2103
606  END IF
607 *
608  2002 CONTINUE
609 * END q-LOOP
610 *
611  2103 CONTINUE
612 * bailed out of q-loop
613 *
614  sva( p ) = aapp
615 *
616  ELSE
617  sva( p ) = aapp
618  IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
619  $ notrot = notrot + min( igl+kbl-1, n ) - p
620  END IF
621 *
622  2001 CONTINUE
623 * end of the p-loop
624 * end of doing the block ( ibr, ibr )
625  1002 CONTINUE
626 * end of ir1-loop
627 *
628 * ... go to the off diagonal blocks
629 *
630  igl = ( ibr-1 )*kbl + 1
631 *
632  DO 2010 jbc = ibr + 1, nbl
633 *
634  jgl = ( jbc-1 )*kbl + 1
635 *
636 * doing the block at ( ibr, jbc )
637 *
638  ijblsk = 0
639  DO 2100 p = igl, min( igl+kbl-1, n )
640 *
641  aapp = sva( p )
642  IF( aapp.GT.zero ) THEN
643 *
644  pskipped = 0
645 *
646  DO 2200 q = jgl, min( jgl+kbl-1, n )
647 *
648  aaqq = sva( q )
649  IF( aaqq.GT.zero ) THEN
650  aapp0 = aapp
651 *
652 * .. M x 2 Jacobi SVD ..
653 *
654 * Safe Gram matrix computation
655 *
656  IF( aaqq.GE.one ) THEN
657  IF( aapp.GE.aaqq ) THEN
658  rotok = ( small*aapp ).LE.aaqq
659  ELSE
660  rotok = ( small*aaqq ).LE.aapp
661  END IF
662  IF( aapp.LT.( big / aaqq ) ) THEN
663  aapq = ( cdotc( m, a( 1, p ), 1,
664  $ a( 1, q ), 1 ) / aaqq ) / aapp
665  ELSE
666  CALL ccopy( m, a( 1, p ), 1,
667  $ work, 1 )
668  CALL clascl( 'G', 0, 0, aapp,
669  $ one, m, 1,
670  $ work, lda, ierr )
671  aapq = cdotc( m, work, 1,
672  $ a( 1, q ), 1 ) / aaqq
673  END IF
674  ELSE
675  IF( aapp.GE.aaqq ) THEN
676  rotok = aapp.LE.( aaqq / small )
677  ELSE
678  rotok = aaqq.LE.( aapp / small )
679  END IF
680  IF( aapp.GT.( small / aaqq ) ) THEN
681  aapq = ( cdotc( m, a( 1, p ), 1,
682  $ a( 1, q ), 1 ) / max(aaqq,aapp) )
683  $ / min(aaqq,aapp)
684  ELSE
685  CALL ccopy( m, a( 1, q ), 1,
686  $ work, 1 )
687  CALL clascl( 'G', 0, 0, aaqq,
688  $ one, m, 1,
689  $ work, lda, ierr )
690  aapq = cdotc( m, a( 1, p ), 1,
691  $ work, 1 ) / aapp
692  END IF
693  END IF
694 *
695 * AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
696  aapq1 = -abs(aapq)
697  mxaapq = max( mxaapq, -aapq1 )
698 *
699 * TO rotate or NOT to rotate, THAT is the question ...
700 *
701  IF( abs( aapq1 ).GT.tol ) THEN
702  ompq = aapq / abs(aapq)
703  notrot = 0
704 *[RTD] ROTATED = ROTATED + 1
705  pskipped = 0
706  iswrot = iswrot + 1
707 *
708  IF( rotok ) THEN
709 *
710  aqoap = aaqq / aapp
711  apoaq = aapp / aaqq
712  theta = -half*abs( aqoap-apoaq )/ aapq1
713  IF( aaqq.GT.aapp0 )theta = -theta
714 *
715  IF( abs( theta ).GT.bigtheta ) THEN
716  t = half / theta
717  cs = one
718  CALL crot( m, a(1,p), 1, a(1,q), 1,
719  $ cs, conjg(ompq)*t )
720  IF( rsvec ) THEN
721  CALL crot( mvl, v(1,p), 1,
722  $ v(1,q), 1, cs, conjg(ompq)*t )
723  END IF
724  sva( q ) = aaqq*sqrt( max( zero,
725  $ one+t*apoaq*aapq1 ) )
726  aapp = aapp*sqrt( max( zero,
727  $ one-t*aqoap*aapq1 ) )
728  mxsinj = max( mxsinj, abs( t ) )
729  ELSE
730 *
731 * .. choose correct signum for THETA and rotate
732 *
733  thsign = -sign( one, aapq1 )
734  IF( aaqq.GT.aapp0 )thsign = -thsign
735  t = one / ( theta+thsign*
736  $ sqrt( one+theta*theta ) )
737  cs = sqrt( one / ( one+t*t ) )
738  sn = t*cs
739  mxsinj = max( mxsinj, abs( sn ) )
740  sva( q ) = aaqq*sqrt( max( zero,
741  $ one+t*apoaq*aapq1 ) )
742  aapp = aapp*sqrt( max( zero,
743  $ one-t*aqoap*aapq1 ) )
744 *
745  CALL crot( m, a(1,p), 1, a(1,q), 1,
746  $ cs, conjg(ompq)*sn )
747  IF( rsvec ) THEN
748  CALL crot( mvl, v(1,p), 1,
749  $ v(1,q), 1, cs, conjg(ompq)*sn )
750  END IF
751  END IF
752  d(p) = -d(q) * ompq
753 *
754  ELSE
755 * .. have to use modified Gram-Schmidt like transformation
756  IF( aapp.GT.aaqq ) THEN
757  CALL ccopy( m, a( 1, p ), 1,
758  $ work, 1 )
759  CALL clascl( 'G', 0, 0, aapp, one,
760  $ m, 1, work,lda,
761  $ ierr )
762  CALL clascl( 'G', 0, 0, aaqq, one,
763  $ m, 1, a( 1, q ), lda,
764  $ ierr )
765  CALL caxpy( m, -aapq, work,
766  $ 1, a( 1, q ), 1 )
767  CALL clascl( 'G', 0, 0, one, aaqq,
768  $ m, 1, a( 1, q ), lda,
769  $ ierr )
770  sva( q ) = aaqq*sqrt( max( zero,
771  $ one-aapq1*aapq1 ) )
772  mxsinj = max( mxsinj, sfmin )
773  ELSE
774  CALL ccopy( m, a( 1, q ), 1,
775  $ work, 1 )
776  CALL clascl( 'G', 0, 0, aaqq, one,
777  $ m, 1, work,lda,
778  $ ierr )
779  CALL clascl( 'G', 0, 0, aapp, one,
780  $ m, 1, a( 1, p ), lda,
781  $ ierr )
782  CALL caxpy( m, -conjg(aapq),
783  $ work, 1, a( 1, p ), 1 )
784  CALL clascl( 'G', 0, 0, one, aapp,
785  $ m, 1, a( 1, p ), lda,
786  $ ierr )
787  sva( p ) = aapp*sqrt( max( zero,
788  $ one-aapq1*aapq1 ) )
789  mxsinj = max( mxsinj, sfmin )
790  END IF
791  END IF
792 * END IF ROTOK THEN ... ELSE
793 *
794 * In the case of cancellation in updating SVA(q), SVA(p)
795 * .. recompute SVA(q), SVA(p)
796  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
797  $ THEN
798  IF( ( aaqq.LT.rootbig ) .AND.
799  $ ( aaqq.GT.rootsfmin ) ) THEN
800  sva( q ) = scnrm2( m, a( 1, q ), 1)
801  ELSE
802  t = zero
803  aaqq = one
804  CALL classq( m, a( 1, q ), 1, t,
805  $ aaqq )
806  sva( q ) = t*sqrt( aaqq )
807  END IF
808  END IF
809  IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
810  IF( ( aapp.LT.rootbig ) .AND.
811  $ ( aapp.GT.rootsfmin ) ) THEN
812  aapp = scnrm2( m, a( 1, p ), 1 )
813  ELSE
814  t = zero
815  aapp = one
816  CALL classq( m, a( 1, p ), 1, t,
817  $ aapp )
818  aapp = t*sqrt( aapp )
819  END IF
820  sva( p ) = aapp
821  END IF
822 * end of OK rotation
823  ELSE
824  notrot = notrot + 1
825 *[RTD] SKIPPED = SKIPPED + 1
826  pskipped = pskipped + 1
827  ijblsk = ijblsk + 1
828  END IF
829  ELSE
830  notrot = notrot + 1
831  pskipped = pskipped + 1
832  ijblsk = ijblsk + 1
833  END IF
834 *
835  IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
836  $ THEN
837  sva( p ) = aapp
838  notrot = 0
839  GO TO 2011
840  END IF
841  IF( ( i.LE.swband ) .AND.
842  $ ( pskipped.GT.rowskip ) ) THEN
843  aapp = -aapp
844  notrot = 0
845  GO TO 2203
846  END IF
847 *
848  2200 CONTINUE
849 * end of the q-loop
850  2203 CONTINUE
851 *
852  sva( p ) = aapp
853 *
854  ELSE
855 *
856  IF( aapp.EQ.zero )notrot = notrot +
857  $ min( jgl+kbl-1, n ) - jgl + 1
858  IF( aapp.LT.zero )notrot = 0
859 *
860  END IF
861 *
862  2100 CONTINUE
863 * end of the p-loop
864  2010 CONTINUE
865 * end of the jbc-loop
866  2011 CONTINUE
867 *2011 bailed out of the jbc-loop
868  DO 2012 p = igl, min( igl+kbl-1, n )
869  sva( p ) = abs( sva( p ) )
870  2012 CONTINUE
871 ***
872  2000 CONTINUE
873 *2000 :: end of the ibr-loop
874 *
875 * .. update SVA(N)
876  IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
877  $ THEN
878  sva( n ) = scnrm2( m, a( 1, n ), 1 )
879  ELSE
880  t = zero
881  aapp = one
882  CALL classq( m, a( 1, n ), 1, t, aapp )
883  sva( n ) = t*sqrt( aapp )
884  END IF
885 *
886 * Additional steering devices
887 *
888  IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
889  $ ( iswrot.LE.n ) ) )swband = i
890 *
891  IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( REAL( N ) )*
892  $ tol ) .AND. ( REAL( n )*MXAAPQ*MXSINJ.LT.TOL ) ) then
893  GO TO 1994
894  END IF
895 *
896  IF( notrot.GE.emptsw )GO TO 1994
897 *
898  1993 CONTINUE
899 * end i=1:NSWEEP loop
900 *
901 * #:( Reaching this point means that the procedure has not converged.
902  info = nsweep - 1
903  GO TO 1995
904 *
905  1994 CONTINUE
906 * #:) Reaching this point means numerical convergence after the i-th
907 * sweep.
908 *
909  info = 0
910 * #:) INFO = 0 confirms successful iterations.
911  1995 CONTINUE
912 *
913 * Sort the vector SVA() of column norms.
914  DO 5991 p = 1, n - 1
915  q = isamax( n-p+1, sva( p ), 1 ) + p - 1
916  IF( p.NE.q ) THEN
917  temp1 = sva( p )
918  sva( p ) = sva( q )
919  sva( q ) = temp1
920  aapq = d( p )
921  d( p ) = d( q )
922  d( q ) = aapq
923  CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
924  IF( rsvec )CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
925  END IF
926  5991 CONTINUE
927 *
928  RETURN
929 * ..
930 * .. END OF CGSVJ0
931 * ..
932  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
crot
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
clascl
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
cgsvj0
subroutine cgsvj0(JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
CGSVJ0 pre-processor for the routine cgesvj.
Definition: cgsvj0.f:218
cswap
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
caxpy
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88