LAPACK  3.7.0
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 Contributor:
207 * ==================
208 *>
209 *> Zlatko Drmac (Zagreb, Croatia)
210 *>
211 *> \par Bugs, Examples and Comments:
212 * =================================
213 *>
214 *> Please report all bugs and send interesting test examples and comments to
215 *> [email protected]. 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.7.0) --
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, max, conjg, REAL, MIN, 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, 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 = min( 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 = min( 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, min( lkahead, nbl-ibr )
388 *
389  igl = igl + ir1*kbl
390 *
391  DO 2001 p = igl, min( 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, min( 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 ) / aapp ) / aaqq
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 * AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
477  aapq1 = -abs(aapq)
478  mxaapq = max( mxaapq, -aapq1 )
479 *
480 * TO rotate or NOT to rotate, THAT is the question ...
481 *
482  IF( abs( aapq1 ).GT.tol ) THEN
483  ompq = aapq / abs(aapq)
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( max( zero,
513  $ one+t*apoaq*aapq1 ) )
514  aapp = aapp*sqrt( max( zero,
515  $ one-t*aqoap*aapq1 ) )
516  mxsinj = max( 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 = max( mxsinj, abs( sn ) )
529  sva( q ) = aaqq*sqrt( max( zero,
530  $ one+t*apoaq*aapq1 ) )
531  aapp = aapp*sqrt( max( 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( max( zero,
557  $ one-aapq1*aapq1 ) )
558  mxsinj = max( 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 + min( 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, min( 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, min( 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 ) / max(aaqq,aapp) )
686  $ / min(aaqq,aapp)
687  ELSE
688  CALL ccopy( m, a( 1, q ), 1,
689  $ work, 1 )
690  CALL clascl( 'G', 0, 0, aaqq,
691  $ one, m, 1,
692  $ work, lda, ierr )
693  aapq = cdotc( m, a( 1, p ), 1,
694  $ work, 1 ) / aapp
695  END IF
696  END IF
697 *
698 * AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
699  aapq1 = -abs(aapq)
700  mxaapq = max( mxaapq, -aapq1 )
701 *
702 * TO rotate or NOT to rotate, THAT is the question ...
703 *
704  IF( abs( aapq1 ).GT.tol ) THEN
705  ompq = aapq / abs(aapq)
706  notrot = 0
707 *[RTD] ROTATED = ROTATED + 1
708  pskipped = 0
709  iswrot = iswrot + 1
710 *
711  IF( rotok ) THEN
712 *
713  aqoap = aaqq / aapp
714  apoaq = aapp / aaqq
715  theta = -half*abs( aqoap-apoaq )/ aapq1
716  IF( aaqq.GT.aapp0 )theta = -theta
717 *
718  IF( abs( theta ).GT.bigtheta ) THEN
719  t = half / theta
720  cs = one
721  CALL crot( m, a(1,p), 1, a(1,q), 1,
722  $ cs, conjg(ompq)*t )
723  IF( rsvec ) THEN
724  CALL crot( mvl, v(1,p), 1,
725  $ v(1,q), 1, cs, conjg(ompq)*t )
726  END IF
727  sva( q ) = aaqq*sqrt( max( zero,
728  $ one+t*apoaq*aapq1 ) )
729  aapp = aapp*sqrt( max( zero,
730  $ one-t*aqoap*aapq1 ) )
731  mxsinj = max( mxsinj, abs( t ) )
732  ELSE
733 *
734 * .. choose correct signum for THETA and rotate
735 *
736  thsign = -sign( one, aapq1 )
737  IF( aaqq.GT.aapp0 )thsign = -thsign
738  t = one / ( theta+thsign*
739  $ sqrt( one+theta*theta ) )
740  cs = sqrt( one / ( one+t*t ) )
741  sn = t*cs
742  mxsinj = max( mxsinj, abs( sn ) )
743  sva( q ) = aaqq*sqrt( max( zero,
744  $ one+t*apoaq*aapq1 ) )
745  aapp = aapp*sqrt( max( zero,
746  $ one-t*aqoap*aapq1 ) )
747 *
748  CALL crot( m, a(1,p), 1, a(1,q), 1,
749  $ cs, conjg(ompq)*sn )
750  IF( rsvec ) THEN
751  CALL crot( mvl, v(1,p), 1,
752  $ v(1,q), 1, cs, conjg(ompq)*sn )
753  END IF
754  END IF
755  d(p) = -d(q) * ompq
756 *
757  ELSE
758 * .. have to use modified Gram-Schmidt like transformation
759  IF( aapp.GT.aaqq ) THEN
760  CALL ccopy( m, a( 1, p ), 1,
761  $ work, 1 )
762  CALL clascl( 'G', 0, 0, aapp, one,
763  $ m, 1, work,lda,
764  $ ierr )
765  CALL clascl( 'G', 0, 0, aaqq, one,
766  $ m, 1, a( 1, q ), lda,
767  $ ierr )
768  CALL caxpy( m, -aapq, work,
769  $ 1, a( 1, q ), 1 )
770  CALL clascl( 'G', 0, 0, one, aaqq,
771  $ m, 1, a( 1, q ), lda,
772  $ ierr )
773  sva( q ) = aaqq*sqrt( max( zero,
774  $ one-aapq1*aapq1 ) )
775  mxsinj = max( mxsinj, sfmin )
776  ELSE
777  CALL ccopy( m, a( 1, q ), 1,
778  $ work, 1 )
779  CALL clascl( 'G', 0, 0, aaqq, one,
780  $ m, 1, work,lda,
781  $ ierr )
782  CALL clascl( 'G', 0, 0, aapp, one,
783  $ m, 1, a( 1, p ), lda,
784  $ ierr )
785  CALL caxpy( m, -conjg(aapq),
786  $ work, 1, a( 1, p ), 1 )
787  CALL clascl( 'G', 0, 0, one, aapp,
788  $ m, 1, a( 1, p ), lda,
789  $ ierr )
790  sva( p ) = aapp*sqrt( max( zero,
791  $ one-aapq1*aapq1 ) )
792  mxsinj = max( mxsinj, sfmin )
793  END IF
794  END IF
795 * END IF ROTOK THEN ... ELSE
796 *
797 * In the case of cancellation in updating SVA(q), SVA(p)
798 * .. recompute SVA(q), SVA(p)
799  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
800  $ THEN
801  IF( ( aaqq.LT.rootbig ) .AND.
802  $ ( aaqq.GT.rootsfmin ) ) THEN
803  sva( q ) = scnrm2( m, a( 1, q ), 1)
804  ELSE
805  t = zero
806  aaqq = one
807  CALL classq( m, a( 1, q ), 1, t,
808  $ aaqq )
809  sva( q ) = t*sqrt( aaqq )
810  END IF
811  END IF
812  IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
813  IF( ( aapp.LT.rootbig ) .AND.
814  $ ( aapp.GT.rootsfmin ) ) THEN
815  aapp = scnrm2( m, a( 1, p ), 1 )
816  ELSE
817  t = zero
818  aapp = one
819  CALL classq( m, a( 1, p ), 1, t,
820  $ aapp )
821  aapp = t*sqrt( aapp )
822  END IF
823  sva( p ) = aapp
824  END IF
825 * end of OK rotation
826  ELSE
827  notrot = notrot + 1
828 *[RTD] SKIPPED = SKIPPED + 1
829  pskipped = pskipped + 1
830  ijblsk = ijblsk + 1
831  END IF
832  ELSE
833  notrot = notrot + 1
834  pskipped = pskipped + 1
835  ijblsk = ijblsk + 1
836  END IF
837 *
838  IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
839  $ THEN
840  sva( p ) = aapp
841  notrot = 0
842  GO TO 2011
843  END IF
844  IF( ( i.LE.swband ) .AND.
845  $ ( pskipped.GT.rowskip ) ) THEN
846  aapp = -aapp
847  notrot = 0
848  GO TO 2203
849  END IF
850 *
851  2200 CONTINUE
852 * end of the q-loop
853  2203 CONTINUE
854 *
855  sva( p ) = aapp
856 *
857  ELSE
858 *
859  IF( aapp.EQ.zero )notrot = notrot +
860  $ min( jgl+kbl-1, n ) - jgl + 1
861  IF( aapp.LT.zero )notrot = 0
862 *
863  END IF
864 *
865  2100 CONTINUE
866 * end of the p-loop
867  2010 CONTINUE
868 * end of the jbc-loop
869  2011 CONTINUE
870 *2011 bailed out of the jbc-loop
871  DO 2012 p = igl, min( igl+kbl-1, n )
872  sva( p ) = abs( sva( p ) )
873  2012 CONTINUE
874 ***
875  2000 CONTINUE
876 *2000 :: end of the ibr-loop
877 *
878 * .. update SVA(N)
879  IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
880  $ THEN
881  sva( n ) = scnrm2( m, a( 1, n ), 1 )
882  ELSE
883  t = zero
884  aapp = one
885  CALL classq( m, a( 1, n ), 1, t, aapp )
886  sva( n ) = t*sqrt( aapp )
887  END IF
888 *
889 * Additional steering devices
890 *
891  IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
892  $ ( iswrot.LE.n ) ) )swband = i
893 *
894  IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( REAL( N ) )*
895  $ tol ) .AND. ( REAL( n )*MXAAPQ*MXSINJ.LT.TOL ) ) then
896  GO TO 1994
897  END IF
898 *
899  IF( notrot.GE.emptsw )GO TO 1994
900 *
901  1993 CONTINUE
902 * end i=1:NSWEEP loop
903 *
904 * #:( Reaching this point means that the procedure has not converged.
905  info = nsweep - 1
906  GO TO 1995
907 *
908  1994 CONTINUE
909 * #:) Reaching this point means numerical convergence after the i-th
910 * sweep.
911 *
912  info = 0
913 * #:) INFO = 0 confirms successful iterations.
914  1995 CONTINUE
915 *
916 * Sort the vector SVA() of column norms.
917  DO 5991 p = 1, n - 1
918  q = isamax( n-p+1, sva( p ), 1 ) + p - 1
919  IF( p.NE.q ) THEN
920  temp1 = sva( p )
921  sva( p ) = sva( q )
922  sva( q ) = temp1
923  aapq = d( p )
924  d( p ) = d( q )
925  d( q ) = aapq
926  CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
927  IF( rsvec )CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
928  END IF
929  5991 CONTINUE
930 *
931  RETURN
932 * ..
933 * .. END OF CGSVJ0
934 * ..
935  END
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
classq
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f:108
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
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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