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