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