LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
zgsvj1.f
Go to the documentation of this file.
1 *> \brief \b ZGSVJ1 pre-processor for the routine zgesvj, applies Jacobi rotations targeting only particular pivots.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
22 * EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * DOUBLE PRECISION EPS, SFMIN, TOL
26 * INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
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 *> ZGSVJ1 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 targets only particular pivots and it does not check convergence
43 *> (stopping criterion). Few tunning parameters (marked by [TP]) are
44 *> available for the implementer.
45 *>
46 *> Further Details
47 *> ~~~~~~~~~~~~~~~
48 *> ZGSVJ1 applies few sweeps of Jacobi rotations in the column space of
49 *> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
50 *> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
51 *> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
52 *> [x]'s in the following scheme:
53 *>
54 *> | * * * [x] [x] [x]|
55 *> | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
56 *> | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
57 *> |[x] [x] [x] * * * |
58 *> |[x] [x] [x] * * * |
59 *> |[x] [x] [x] * * * |
60 *>
61 *> In terms of the columns of A, the first N1 columns are rotated 'against'
62 *> the remaining N-N1 columns, trying to increase the angle between the
63 *> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
64 *> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
65 *> The number of sweeps is given in NSWEEP and the orthogonality threshold
66 *> is given in TOL.
67 *> \endverbatim
68 *
69 * Arguments:
70 * ==========
71 *
72 *> \param[in] JOBV
73 *> \verbatim
74 *> JOBV is CHARACTER*1
75 *> Specifies whether the output from this procedure is used
76 *> to compute the matrix V:
77 *> = 'V': the product of the Jacobi rotations is accumulated
78 *> by postmulyiplying the N-by-N array V.
79 *> (See the description of V.)
80 *> = 'A': the product of the Jacobi rotations is accumulated
81 *> by postmulyiplying the MV-by-N array V.
82 *> (See the descriptions of MV and V.)
83 *> = 'N': the Jacobi rotations are not accumulated.
84 *> \endverbatim
85 *>
86 *> \param[in] M
87 *> \verbatim
88 *> M is INTEGER
89 *> The number of rows of the input matrix A. M >= 0.
90 *> \endverbatim
91 *>
92 *> \param[in] N
93 *> \verbatim
94 *> N is INTEGER
95 *> The number of columns of the input matrix A.
96 *> M >= N >= 0.
97 *> \endverbatim
98 *>
99 *> \param[in] N1
100 *> \verbatim
101 *> N1 is INTEGER
102 *> N1 specifies the 2 x 2 block partition, the first N1 columns are
103 *> rotated 'against' the remaining N-N1 columns of A.
104 *> \endverbatim
105 *>
106 *> \param[in,out] A
107 *> \verbatim
108 *> A is COMPLEX*16 array, dimension (LDA,N)
109 *> On entry, M-by-N matrix A, such that A*diag(D) represents
110 *> the input matrix.
111 *> On exit,
112 *> A_onexit * D_onexit represents the input matrix A*diag(D)
113 *> post-multiplied by a sequence of Jacobi rotations, where the
114 *> rotation threshold and the total number of sweeps are given in
115 *> TOL and NSWEEP, respectively.
116 *> (See the descriptions of N1, D, TOL and NSWEEP.)
117 *> \endverbatim
118 *>
119 *> \param[in] LDA
120 *> \verbatim
121 *> LDA is INTEGER
122 *> The leading dimension of the array A. LDA >= max(1,M).
123 *> \endverbatim
124 *>
125 *> \param[in,out] D
126 *> \verbatim
127 *> D is COMPLEX*16 array, dimension (N)
128 *> The array D accumulates the scaling factors from the fast scaled
129 *> Jacobi rotations.
130 *> On entry, A*diag(D) represents the input matrix.
131 *> On exit, A_onexit*diag(D_onexit) represents the input matrix
132 *> post-multiplied by a sequence of Jacobi rotations, where the
133 *> rotation threshold and the total number of sweeps are given in
134 *> TOL and NSWEEP, respectively.
135 *> (See the descriptions of N1, A, TOL and NSWEEP.)
136 *> \endverbatim
137 *>
138 *> \param[in,out] SVA
139 *> \verbatim
140 *> SVA is DOUBLE PRECISION array, dimension (N)
141 *> On entry, SVA contains the Euclidean norms of the columns of
142 *> the matrix A*diag(D).
143 *> On exit, SVA contains the Euclidean norms of the columns of
144 *> the matrix onexit*diag(D_onexit).
145 *> \endverbatim
146 *>
147 *> \param[in] MV
148 *> \verbatim
149 *> MV is INTEGER
150 *> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
151 *> sequence of Jacobi rotations.
152 *> If JOBV = 'N', then MV is not referenced.
153 *> \endverbatim
154 *>
155 *> \param[in,out] V
156 *> \verbatim
157 *> V is COMPLEX*16 array, dimension (LDV,N)
158 *> If JOBV .EQ. 'V' then N rows of V are post-multipled by a
159 *> sequence of Jacobi rotations.
160 *> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
161 *> sequence of Jacobi rotations.
162 *> If JOBV = 'N', then V is not referenced.
163 *> \endverbatim
164 *>
165 *> \param[in] LDV
166 *> \verbatim
167 *> LDV is INTEGER
168 *> The leading dimension of the array V, LDV >= 1.
169 *> If JOBV = 'V', LDV .GE. N.
170 *> If JOBV = 'A', LDV .GE. MV.
171 *> \endverbatim
172 *>
173 *> \param[in] EPS
174 *> \verbatim
175 *> EPS is DOUBLE PRECISION
176 *> EPS = DLAMCH('Epsilon')
177 *> \endverbatim
178 *>
179 *> \param[in] SFMIN
180 *> \verbatim
181 *> SFMIN is DOUBLE PRECISION
182 *> SFMIN = DLAMCH('Safe Minimum')
183 *> \endverbatim
184 *>
185 *> \param[in] TOL
186 *> \verbatim
187 *> TOL is DOUBLE PRECISION
188 *> TOL is the threshold for Jacobi rotations. For a pair
189 *> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
190 *> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
191 *> \endverbatim
192 *>
193 *> \param[in] NSWEEP
194 *> \verbatim
195 *> NSWEEP is INTEGER
196 *> NSWEEP is the number of sweeps of Jacobi rotations to be
197 *> performed.
198 *> \endverbatim
199 *>
200 *> \param[out] WORK
201 *> \verbatim
202 *> WORK is COMPLEX*16 array, dimension (LWORK)
203 *> \endverbatim
204 *>
205 *> \param[in] LWORK
206 *> \verbatim
207 *> LWORK is INTEGER
208 *> LWORK is the dimension of WORK. LWORK .GE. M.
209 *> \endverbatim
210 *>
211 *> \param[out] INFO
212 *> \verbatim
213 *> INFO is INTEGER
214 *> = 0 : successful exit.
215 *> < 0 : if INFO = -i, then the i-th argument had an illegal value
216 *> \endverbatim
217 *
218 * Authors:
219 * ========
220 *
221 *> \author Univ. of Tennessee
222 *> \author Univ. of California Berkeley
223 *> \author Univ. of Colorado Denver
224 *> \author NAG Ltd.
225 *
226 *> \date June 2016
227 *
228 *> \ingroup complex16OTHERcomputational
229 *
230 *> \par Contributors:
231 * ==================
232 *>
233 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
234 *
235 * =====================================================================
236  SUBROUTINE zgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
237  \$ eps, sfmin, tol, nsweep, work, lwork, info )
238 *
239 * -- LAPACK computational routine (version 3.6.1) --
240 * -- LAPACK is a software package provided by Univ. of Tennessee, --
241 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242 * June 2016
243 *
244  IMPLICIT NONE
245 * .. Scalar Arguments ..
246  DOUBLE PRECISION EPS, SFMIN, TOL
247  INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
248  CHARACTER*1 JOBV
249 * ..
250 * .. Array Arguments ..
251  COMPLEX*16 A( lda, * ), D( n ), V( ldv, * ), WORK( lwork )
252  DOUBLE PRECISION SVA( n )
253 * ..
254 *
255 * =====================================================================
256 *
257 * .. Local Parameters ..
258  DOUBLE PRECISION ZERO, HALF, ONE
259  parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
260 * ..
261 * .. Local Scalars ..
262  COMPLEX*16 AAPQ, OMPQ
263  DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
264  \$ bigtheta, cs, large, mxaapq, mxsinj, rootbig,
265  \$ rooteps, rootsfmin, roottol, small, sn, t,
266  \$ temp1, theta, thsign
267  INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
268  \$ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
269  \$ p, pskipped, q, rowskip, swband
270  LOGICAL APPLV, ROTOK, RSVEC
271 * ..
272 * ..
273 * .. Intrinsic Functions ..
274  INTRINSIC abs, dconjg, dmax1, dble, min0, dsign, dsqrt
275 * ..
276 * .. External Functions ..
277  DOUBLE PRECISION DZNRM2
278  COMPLEX*16 ZDOTC
279  INTEGER IDAMAX
280  LOGICAL LSAME
281  EXTERNAL idamax, lsame, zdotc, dznrm2
282 * ..
283 * .. External Subroutines ..
284 * .. from BLAS
285  EXTERNAL zcopy, zrot, zswap
286 * .. from LAPACK
287  EXTERNAL zlascl, zlassq, xerbla
288 * ..
289 * .. Executable Statements ..
290 *
291 * Test the input parameters.
292 *
293  applv = lsame( jobv, 'A' )
294  rsvec = lsame( jobv, 'V' )
295  IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
296  info = -1
297  ELSE IF( m.LT.0 ) THEN
298  info = -2
299  ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
300  info = -3
301  ELSE IF( n1.LT.0 ) THEN
302  info = -4
303  ELSE IF( lda.LT.m ) THEN
304  info = -6
305  ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
306  info = -9
307  ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
308  \$ ( applv.AND.( ldv.LT.mv ) ) ) THEN
309  info = -11
310  ELSE IF( tol.LE.eps ) THEN
311  info = -14
312  ELSE IF( nsweep.LT.0 ) THEN
313  info = -15
314  ELSE IF( lwork.LT.m ) THEN
315  info = -17
316  ELSE
317  info = 0
318  END IF
319 *
320 * #:(
321  IF( info.NE.0 ) THEN
322  CALL xerbla( 'ZGSVJ1', -info )
323  RETURN
324  END IF
325 *
326  IF( rsvec ) THEN
327  mvl = n
328  ELSE IF( applv ) THEN
329  mvl = mv
330  END IF
331  rsvec = rsvec .OR. applv
332
333  rooteps = dsqrt( eps )
334  rootsfmin = dsqrt( sfmin )
335  small = sfmin / eps
336  big = one / sfmin
337  rootbig = one / rootsfmin
338  large = big / dsqrt( dble( m*n ) )
339  bigtheta = one / rooteps
340  roottol = dsqrt( tol )
341 *
342 * .. Initialize the right singular vector matrix ..
343 *
344 * RSVEC = LSAME( JOBV, 'Y' )
345 *
346  emptsw = n1*( n-n1 )
347  notrot = 0
348 *
349 * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
350 *
351  kbl = min0( 8, n )
352  nblr = n1 / kbl
353  IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
354
355 * .. the tiling is nblr-by-nblc [tiles]
356
357  nblc = ( n-n1 ) / kbl
358  IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
359  blskip = ( kbl**2 ) + 1
360 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
361
362  rowskip = min0( 5, kbl )
363 *[TP] ROWSKIP is a tuning parameter.
364  swband = 0
365 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
366 * if ZGESVJ is used as a computational routine in the preconditioned
367 * Jacobi SVD algorithm ZGEJSV.
368 *
369 *
370 * | * * * [x] [x] [x]|
371 * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
372 * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
373 * |[x] [x] [x] * * * |
374 * |[x] [x] [x] * * * |
375 * |[x] [x] [x] * * * |
376 *
377 *
378  DO 1993 i = 1, nsweep
379 *
380 * .. go go go ...
381 *
382  mxaapq = zero
383  mxsinj = zero
384  iswrot = 0
385 *
386  notrot = 0
387  pskipped = 0
388 *
389 * Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
390 * 1 <= p < q <= N. This is the first step toward a blocked implementation
391 * of the rotations. New implementation, based on block transformations,
392 * is under development.
393 *
394  DO 2000 ibr = 1, nblr
395 *
396  igl = ( ibr-1 )*kbl + 1
397 *
398
399 *
400 * ... go to the off diagonal blocks
401 *
402  igl = ( ibr-1 )*kbl + 1
403 *
404 * DO 2010 jbc = ibr + 1, NBL
405  DO 2010 jbc = 1, nblc
406 *
407  jgl = ( jbc-1 )*kbl + n1 + 1
408 *
409 * doing the block at ( ibr, jbc )
410 *
411  ijblsk = 0
412  DO 2100 p = igl, min0( igl+kbl-1, n1 )
413 *
414  aapp = sva( p )
415  IF( aapp.GT.zero ) THEN
416 *
417  pskipped = 0
418 *
419  DO 2200 q = jgl, min0( jgl+kbl-1, n )
420 *
421  aaqq = sva( q )
422  IF( aaqq.GT.zero ) THEN
423  aapp0 = aapp
424 *
425 * .. M x 2 Jacobi SVD ..
426 *
427 * Safe Gram matrix computation
428 *
429  IF( aaqq.GE.one ) THEN
430  IF( aapp.GE.aaqq ) THEN
431  rotok = ( small*aapp ).LE.aaqq
432  ELSE
433  rotok = ( small*aaqq ).LE.aapp
434  END IF
435  IF( aapp.LT.( big / aaqq ) ) THEN
436  aapq = ( zdotc( m, a( 1, p ), 1,
437  \$ a( 1, q ), 1 ) / aaqq ) / aapp
438  ELSE
439  CALL zcopy( m, a( 1, p ), 1,
440  \$ work, 1 )
441  CALL zlascl( 'G', 0, 0, aapp,
442  \$ one, m, 1,
443  \$ work, lda, ierr )
444  aapq = zdotc( m, work, 1,
445  \$ a( 1, q ), 1 ) / aaqq
446  END IF
447  ELSE
448  IF( aapp.GE.aaqq ) THEN
449  rotok = aapp.LE.( aaqq / small )
450  ELSE
451  rotok = aaqq.LE.( aapp / small )
452  END IF
453  IF( aapp.GT.( small / aaqq ) ) THEN
454  aapq = ( zdotc( m, a( 1, p ), 1,
455  \$ a( 1, q ), 1 ) / aaqq ) / aapp
456  ELSE
457  CALL zcopy( m, a( 1, q ), 1,
458  \$ work, 1 )
459  CALL zlascl( 'G', 0, 0, aaqq,
460  \$ one, m, 1,
461  \$ work, lda, ierr )
462  aapq = zdotc( m, a( 1, p ), 1,
463  \$ work, 1 ) / aapp
464  END IF
465  END IF
466 *
467  ompq = aapq / abs(aapq)
468 * AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)
469  aapq1 = -abs(aapq)
470  mxaapq = dmax1( mxaapq, -aapq1 )
471 *
472 * TO rotate or NOT to rotate, THAT is the question ...
473 *
474  IF( abs( aapq1 ).GT.tol ) THEN
475  notrot = 0
476 *[RTD] ROTATED = ROTATED + 1
477  pskipped = 0
478  iswrot = iswrot + 1
479 *
480  IF( rotok ) THEN
481 *
482  aqoap = aaqq / aapp
483  apoaq = aapp / aaqq
484  theta = -half*abs( aqoap-apoaq )/ aapq1
485  IF( aaqq.GT.aapp0 )theta = -theta
486 *
487  IF( abs( theta ).GT.bigtheta ) THEN
488  t = half / theta
489  cs = one
490  CALL zrot( m, a(1,p), 1, a(1,q), 1,
491  \$ cs, dconjg(ompq)*t )
492  IF( rsvec ) THEN
493  CALL zrot( mvl, v(1,p), 1,
494  \$ v(1,q), 1, cs, dconjg(ompq)*t )
495  END IF
496  sva( q ) = aaqq*dsqrt( dmax1( zero,
497  \$ one+t*apoaq*aapq1 ) )
498  aapp = aapp*dsqrt( dmax1( zero,
499  \$ one-t*aqoap*aapq1 ) )
500  mxsinj = dmax1( mxsinj, abs( t ) )
501  ELSE
502 *
503 * .. choose correct signum for THETA and rotate
504 *
505  thsign = -dsign( one, aapq1 )
506  IF( aaqq.GT.aapp0 )thsign = -thsign
507  t = one / ( theta+thsign*
508  \$ dsqrt( one+theta*theta ) )
509  cs = dsqrt( one / ( one+t*t ) )
510  sn = t*cs
511  mxsinj = dmax1( mxsinj, abs( sn ) )
512  sva( q ) = aaqq*dsqrt( dmax1( zero,
513  \$ one+t*apoaq*aapq1 ) )
514  aapp = aapp*dsqrt( dmax1( zero,
515  \$ one-t*aqoap*aapq1 ) )
516 *
517  CALL zrot( m, a(1,p), 1, a(1,q), 1,
518  \$ cs, dconjg(ompq)*sn )
519  IF( rsvec ) THEN
520  CALL zrot( mvl, v(1,p), 1,
521  \$ v(1,q), 1, cs, dconjg(ompq)*sn )
522  END IF
523  END IF
524  d(p) = -d(q) * ompq
525 *
526  ELSE
527 * .. have to use modified Gram-Schmidt like transformation
528  IF( aapp.GT.aaqq ) THEN
529  CALL zcopy( m, a( 1, p ), 1,
530  \$ work, 1 )
531  CALL zlascl( 'G', 0, 0, aapp, one,
532  \$ m, 1, work,lda,
533  \$ ierr )
534  CALL zlascl( 'G', 0, 0, aaqq, one,
535  \$ m, 1, a( 1, q ), lda,
536  \$ ierr )
537  CALL zaxpy( m, -aapq, work,
538  \$ 1, a( 1, q ), 1 )
539  CALL zlascl( 'G', 0, 0, one, aaqq,
540  \$ m, 1, a( 1, q ), lda,
541  \$ ierr )
542  sva( q ) = aaqq*dsqrt( dmax1( zero,
543  \$ one-aapq1*aapq1 ) )
544  mxsinj = dmax1( mxsinj, sfmin )
545  ELSE
546  CALL zcopy( m, a( 1, q ), 1,
547  \$ work, 1 )
548  CALL zlascl( 'G', 0, 0, aaqq, one,
549  \$ m, 1, work,lda,
550  \$ ierr )
551  CALL zlascl( 'G', 0, 0, aapp, one,
552  \$ m, 1, a( 1, p ), lda,
553  \$ ierr )
554  CALL zaxpy( m, -dconjg(aapq),
555  \$ work, 1, a( 1, p ), 1 )
556  CALL zlascl( 'G', 0, 0, one, aapp,
557  \$ m, 1, a( 1, p ), lda,
558  \$ ierr )
559  sva( p ) = aapp*dsqrt( dmax1( zero,
560  \$ one-aapq1*aapq1 ) )
561  mxsinj = dmax1( mxsinj, sfmin )
562  END IF
563  END IF
564 * END IF ROTOK THEN ... ELSE
565 *
566 * In the case of cancellation in updating SVA(q), SVA(p)
567 * .. recompute SVA(q), SVA(p)
568  IF( ( sva( q ) / aaqq )**2.LE.rooteps )
569  \$ THEN
570  IF( ( aaqq.LT.rootbig ) .AND.
571  \$ ( aaqq.GT.rootsfmin ) ) THEN
572  sva( q ) = dznrm2( m, a( 1, q ), 1)
573  ELSE
574  t = zero
575  aaqq = one
576  CALL zlassq( m, a( 1, q ), 1, t,
577  \$ aaqq )
578  sva( q ) = t*dsqrt( aaqq )
579  END IF
580  END IF
581  IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
582  IF( ( aapp.LT.rootbig ) .AND.
583  \$ ( aapp.GT.rootsfmin ) ) THEN
584  aapp = dznrm2( m, a( 1, p ), 1 )
585  ELSE
586  t = zero
587  aapp = one
588  CALL zlassq( m, a( 1, p ), 1, t,
589  \$ aapp )
590  aapp = t*dsqrt( aapp )
591  END IF
592  sva( p ) = aapp
593  END IF
594 * end of OK rotation
595  ELSE
596  notrot = notrot + 1
597 *[RTD] SKIPPED = SKIPPED + 1
598  pskipped = pskipped + 1
599  ijblsk = ijblsk + 1
600  END IF
601  ELSE
602  notrot = notrot + 1
603  pskipped = pskipped + 1
604  ijblsk = ijblsk + 1
605  END IF
606 *
607  IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
608  \$ THEN
609  sva( p ) = aapp
610  notrot = 0
611  GO TO 2011
612  END IF
613  IF( ( i.LE.swband ) .AND.
614  \$ ( pskipped.GT.rowskip ) ) THEN
615  aapp = -aapp
616  notrot = 0
617  GO TO 2203
618  END IF
619 *
620  2200 CONTINUE
621 * end of the q-loop
622  2203 CONTINUE
623 *
624  sva( p ) = aapp
625 *
626  ELSE
627 *
628  IF( aapp.EQ.zero )notrot = notrot +
629  \$ min0( jgl+kbl-1, n ) - jgl + 1
630  IF( aapp.LT.zero )notrot = 0
631 *
632  END IF
633 *
634  2100 CONTINUE
635 * end of the p-loop
636  2010 CONTINUE
637 * end of the jbc-loop
638  2011 CONTINUE
639 *2011 bailed out of the jbc-loop
640  DO 2012 p = igl, min0( igl+kbl-1, n )
641  sva( p ) = abs( sva( p ) )
642  2012 CONTINUE
643 ***
644  2000 CONTINUE
645 *2000 :: end of the ibr-loop
646 *
647 * .. update SVA(N)
648  IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
649  \$ THEN
650  sva( n ) = dznrm2( m, a( 1, n ), 1 )
651  ELSE
652  t = zero
653  aapp = one
654  CALL zlassq( m, a( 1, n ), 1, t, aapp )
655  sva( n ) = t*dsqrt( aapp )
656  END IF
657 *
659 *
660  IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
661  \$ ( iswrot.LE.n ) ) )swband = i
662 *
663  IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dsqrt( dble( n ) )*
664  \$ tol ) .AND. ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
665  GO TO 1994
666  END IF
667 *
668  IF( notrot.GE.emptsw )GO TO 1994
669 *
670  1993 CONTINUE
671 * end i=1:NSWEEP loop
672 *
673 * #:( Reaching this point means that the procedure has not converged.
674  info = nsweep - 1
675  GO TO 1995
676 *
677  1994 CONTINUE
678 * #:) Reaching this point means numerical convergence after the i-th
679 * sweep.
680 *
681  info = 0
682 * #:) INFO = 0 confirms successful iterations.
683  1995 CONTINUE
684 *
685 * Sort the vector SVA() of column norms.
686  DO 5991 p = 1, n - 1
687  q = idamax( n-p+1, sva( p ), 1 ) + p - 1
688  IF( p.NE.q ) THEN
689  temp1 = sva( p )
690  sva( p ) = sva( q )
691  sva( q ) = temp1
692  aapq = d( p )
693  d( p ) = d( q )
694  d( q ) = aapq
695  CALL zswap( m, a( 1, p ), 1, a( 1, q ), 1 )
696  IF( rsvec )CALL zswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
697  END IF
698  5991 CONTINUE
699 *
700 *
701  RETURN
702 * ..
703 * .. END OF ZGSVJ1
704 * ..
705  END
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
ZGSVJ1 pre-processor for the routine zgesvj, applies Jacobi rotations targeting only particular pivot...
Definition: zgsvj1.f:238
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
Definition: zlassq.f:108
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
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
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53