LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
cgsvj1.f
Go to the documentation of this file.
1*> \brief \b CGSVJ1 pre-processor for the routine cgesvj, 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/cgsvj1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgsvj1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgsvj1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
22* EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* REAL EPS, SFMIN, TOL
26* INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
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*> CGSVJ1 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 targets only particular pivots and it does not check convergence
43*> (stopping criterion). Few tuning parameters (marked by [TP]) are
44*> available for the implementer.
45*>
46*> Further Details
47*> ~~~~~~~~~~~~~~~
48*> CGSVJ1 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 tuning parameter.
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 postmultiplying the N-by-N array V.
79*> (See the description of V.)
80*> = 'A': the product of the Jacobi rotations is accumulated
81*> by postmultiplying 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 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 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 REAL 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 = 'A', then MV rows of V are post-multiplied 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 array, dimension (LDV,N)
158*> If JOBV = 'V' then N rows of V are post-multiplied by a
159*> sequence of Jacobi rotations.
160*> If JOBV = 'A' then MV rows of V are post-multiplied 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 >= N.
170*> If JOBV = 'A', LDV >= MV.
171*> \endverbatim
172*>
173*> \param[in] EPS
174*> \verbatim
175*> EPS is REAL
176*> EPS = SLAMCH('Epsilon')
177*> \endverbatim
178*>
179*> \param[in] SFMIN
180*> \verbatim
181*> SFMIN is REAL
182*> SFMIN = SLAMCH('Safe Minimum')
183*> \endverbatim
184*>
185*> \param[in] TOL
186*> \verbatim
187*> TOL is REAL
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)))) > 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 array, dimension (LWORK)
203*> \endverbatim
204*>
205*> \param[in] LWORK
206*> \verbatim
207*> LWORK is INTEGER
208*> LWORK is the dimension of WORK. LWORK >= 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*> \ingroup gsvj1
227*
228*> \par Contributor:
229* ==================
230*>
231*> Zlatko Drmac (Zagreb, Croatia)
232*
233* =====================================================================
234 SUBROUTINE cgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
235 \$ EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
236*
237* -- LAPACK computational routine --
238* -- LAPACK is a software package provided by Univ. of Tennessee, --
239* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240*
241* .. Scalar Arguments ..
242 REAL EPS, SFMIN, TOL
243 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
244 CHARACTER*1 JOBV
245* ..
246* .. Array Arguments ..
247 COMPLEX A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
248 REAL SVA( N )
249* ..
250*
251* =====================================================================
252*
253* .. Local Parameters ..
254 REAL ZERO, HALF, ONE
255 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
256* ..
257* .. Local Scalars ..
258 COMPLEX AAPQ, OMPQ
259 REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
260 \$ bigtheta, cs, mxaapq, mxsinj, rootbig,
261 \$ rooteps, rootsfmin, roottol, small, sn, t,
262 \$ temp1, theta, thsign
263 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
264 \$ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
265 \$ p, pskipped, q, rowskip, swband
266 LOGICAL APPLV, ROTOK, RSVEC
267* ..
268* ..
269* .. Intrinsic Functions ..
270 INTRINSIC abs, max, conjg, real, min, sign, sqrt
271* ..
272* .. External Functions ..
273 REAL SCNRM2
274 COMPLEX CDOTC
275 INTEGER ISAMAX
276 LOGICAL LSAME
277 EXTERNAL isamax, lsame, cdotc, scnrm2
278* ..
279* .. External Subroutines ..
280* .. from BLAS
281 EXTERNAL ccopy, crot, cswap, caxpy
282* .. from LAPACK
283 EXTERNAL clascl, classq, xerbla
284* ..
285* .. Executable Statements ..
286*
287* Test the input parameters.
288*
289 applv = lsame( jobv, 'A' )
290 rsvec = lsame( jobv, 'V' )
291 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
292 info = -1
293 ELSE IF( m.LT.0 ) THEN
294 info = -2
295 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
296 info = -3
297 ELSE IF( n1.LT.0 ) THEN
298 info = -4
299 ELSE IF( lda.LT.m ) THEN
300 info = -6
301 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
302 info = -9
303 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
304 \$ ( applv.AND.( ldv.LT.mv ) ) ) THEN
305 info = -11
306 ELSE IF( tol.LE.eps ) THEN
307 info = -14
308 ELSE IF( nsweep.LT.0 ) THEN
309 info = -15
310 ELSE IF( lwork.LT.m ) THEN
311 info = -17
312 ELSE
313 info = 0
314 END IF
315*
316* #:(
317 IF( info.NE.0 ) THEN
318 CALL xerbla( 'CGSVJ1', -info )
319 RETURN
320 END IF
321*
322 IF( rsvec ) THEN
323 mvl = n
324 ELSE IF( applv ) THEN
325 mvl = mv
326 END IF
327 rsvec = rsvec .OR. applv
328
329 rooteps = sqrt( eps )
330 rootsfmin = sqrt( sfmin )
331 small = sfmin / eps
332 big = one / sfmin
333 rootbig = one / rootsfmin
334* LARGE = BIG / SQRT( REAL( M*N ) )
335 bigtheta = one / rooteps
336 roottol = sqrt( tol )
337*
338* .. Initialize the right singular vector matrix ..
339*
340* RSVEC = LSAME( JOBV, 'Y' )
341*
342 emptsw = n1*( n-n1 )
343 notrot = 0
344*
345* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
346*
347 kbl = min( 8, n )
348 nblr = n1 / kbl
349 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
350
351* .. the tiling is nblr-by-nblc [tiles]
352
353 nblc = ( n-n1 ) / kbl
354 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
355 blskip = ( kbl**2 ) + 1
356*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
357
358 rowskip = min( 5, kbl )
359*[TP] ROWSKIP is a tuning parameter.
360 swband = 0
361*[TP] SWBAND is a tuning parameter. It is meaningful and effective
362* if CGESVJ is used as a computational routine in the preconditioned
363* Jacobi SVD algorithm CGEJSV.
364*
365*
366* | * * * [x] [x] [x]|
367* | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
368* | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
369* |[x] [x] [x] * * * |
370* |[x] [x] [x] * * * |
371* |[x] [x] [x] * * * |
372*
373*
374 DO 1993 i = 1, nsweep
375*
376* .. go go go ...
377*
378 mxaapq = zero
379 mxsinj = zero
380 iswrot = 0
381*
382 notrot = 0
383 pskipped = 0
384*
385* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
386* 1 <= p < q <= N. This is the first step toward a blocked implementation
387* of the rotations. New implementation, based on block transformations,
388* is under development.
389*
390 DO 2000 ibr = 1, nblr
391*
392 igl = ( ibr-1 )*kbl + 1
393*
394
395*
396* ... go to the off diagonal blocks
397*
398 igl = ( ibr-1 )*kbl + 1
399*
400* DO 2010 jbc = ibr + 1, NBL
401 DO 2010 jbc = 1, nblc
402*
403 jgl = ( jbc-1 )*kbl + n1 + 1
404*
405* doing the block at ( ibr, jbc )
406*
407 ijblsk = 0
408 DO 2100 p = igl, min( igl+kbl-1, n1 )
409*
410 aapp = sva( p )
411 IF( aapp.GT.zero ) THEN
412*
413 pskipped = 0
414*
415 DO 2200 q = jgl, min( jgl+kbl-1, n )
416*
417 aaqq = sva( q )
418 IF( aaqq.GT.zero ) THEN
419 aapp0 = aapp
420*
421* .. M x 2 Jacobi SVD ..
422*
423* Safe Gram matrix computation
424*
425 IF( aaqq.GE.one ) THEN
426 IF( aapp.GE.aaqq ) THEN
427 rotok = ( small*aapp ).LE.aaqq
428 ELSE
429 rotok = ( small*aaqq ).LE.aapp
430 END IF
431 IF( aapp.LT.( big / aaqq ) ) THEN
432 aapq = ( cdotc( m, a( 1, p ), 1,
433 \$ a( 1, q ), 1 ) / aaqq ) / aapp
434 ELSE
435 CALL ccopy( m, a( 1, p ), 1,
436 \$ work, 1 )
437 CALL clascl( 'G', 0, 0, aapp,
438 \$ one, m, 1,
439 \$ work, lda, ierr )
440 aapq = cdotc( m, work, 1,
441 \$ a( 1, q ), 1 ) / aaqq
442 END IF
443 ELSE
444 IF( aapp.GE.aaqq ) THEN
445 rotok = aapp.LE.( aaqq / small )
446 ELSE
447 rotok = aaqq.LE.( aapp / small )
448 END IF
449 IF( aapp.GT.( small / aaqq ) ) THEN
450 aapq = ( cdotc( m, a( 1, p ), 1,
451 \$ a( 1, q ), 1 ) / max(aaqq,aapp) )
452 \$ / min(aaqq,aapp)
453 ELSE
454 CALL ccopy( m, a( 1, q ), 1,
455 \$ work, 1 )
456 CALL clascl( 'G', 0, 0, aaqq,
457 \$ one, m, 1,
458 \$ work, lda, ierr )
459 aapq = cdotc( m, a( 1, p ), 1,
460 \$ work, 1 ) / aapp
461 END IF
462 END IF
463*
464* AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
465 aapq1 = -abs(aapq)
466 mxaapq = max( mxaapq, -aapq1 )
467*
468* TO rotate or NOT to rotate, THAT is the question ...
469*
470 IF( abs( aapq1 ).GT.tol ) THEN
471 ompq = aapq / abs(aapq)
472 notrot = 0
473*[RTD] ROTATED = ROTATED + 1
474 pskipped = 0
475 iswrot = iswrot + 1
476*
477 IF( rotok ) THEN
478*
479 aqoap = aaqq / aapp
480 apoaq = aapp / aaqq
481 theta = -half*abs( aqoap-apoaq )/ aapq1
482 IF( aaqq.GT.aapp0 )theta = -theta
483*
484 IF( abs( theta ).GT.bigtheta ) THEN
485 t = half / theta
486 cs = one
487 CALL crot( m, a(1,p), 1, a(1,q), 1,
488 \$ cs, conjg(ompq)*t )
489 IF( rsvec ) THEN
490 CALL crot( mvl, v(1,p), 1,
491 \$ v(1,q), 1, cs, conjg(ompq)*t )
492 END IF
493 sva( q ) = aaqq*sqrt( max( zero,
494 \$ one+t*apoaq*aapq1 ) )
495 aapp = aapp*sqrt( max( zero,
496 \$ one-t*aqoap*aapq1 ) )
497 mxsinj = max( mxsinj, abs( t ) )
498 ELSE
499*
500* .. choose correct signum for THETA and rotate
501*
502 thsign = -sign( one, aapq1 )
503 IF( aaqq.GT.aapp0 )thsign = -thsign
504 t = one / ( theta+thsign*
505 \$ sqrt( one+theta*theta ) )
506 cs = sqrt( one / ( one+t*t ) )
507 sn = t*cs
508 mxsinj = max( mxsinj, abs( sn ) )
509 sva( q ) = aaqq*sqrt( max( zero,
510 \$ one+t*apoaq*aapq1 ) )
511 aapp = aapp*sqrt( max( zero,
512 \$ one-t*aqoap*aapq1 ) )
513*
514 CALL crot( m, a(1,p), 1, a(1,q), 1,
515 \$ cs, conjg(ompq)*sn )
516 IF( rsvec ) THEN
517 CALL crot( mvl, v(1,p), 1,
518 \$ v(1,q), 1, cs, conjg(ompq)*sn )
519 END IF
520 END IF
521 d(p) = -d(q) * ompq
522*
523 ELSE
524* .. have to use modified Gram-Schmidt like transformation
525 IF( aapp.GT.aaqq ) THEN
526 CALL ccopy( m, a( 1, p ), 1,
527 \$ work, 1 )
528 CALL clascl( 'G', 0, 0, aapp, one,
529 \$ m, 1, work,lda,
530 \$ ierr )
531 CALL clascl( 'G', 0, 0, aaqq, one,
532 \$ m, 1, a( 1, q ), lda,
533 \$ ierr )
534 CALL caxpy( m, -aapq, work,
535 \$ 1, a( 1, q ), 1 )
536 CALL clascl( 'G', 0, 0, one, aaqq,
537 \$ m, 1, a( 1, q ), lda,
538 \$ ierr )
539 sva( q ) = aaqq*sqrt( max( zero,
540 \$ one-aapq1*aapq1 ) )
541 mxsinj = max( mxsinj, sfmin )
542 ELSE
543 CALL ccopy( m, a( 1, q ), 1,
544 \$ work, 1 )
545 CALL clascl( 'G', 0, 0, aaqq, one,
546 \$ m, 1, work,lda,
547 \$ ierr )
548 CALL clascl( 'G', 0, 0, aapp, one,
549 \$ m, 1, a( 1, p ), lda,
550 \$ ierr )
551 CALL caxpy( m, -conjg(aapq),
552 \$ work, 1, a( 1, p ), 1 )
553 CALL clascl( 'G', 0, 0, one, aapp,
554 \$ m, 1, a( 1, p ), lda,
555 \$ ierr )
556 sva( p ) = aapp*sqrt( max( zero,
557 \$ one-aapq1*aapq1 ) )
558 mxsinj = max( mxsinj, sfmin )
559 END IF
560 END IF
561* END IF ROTOK THEN ... ELSE
562*
563* In the case of cancellation in updating SVA(q), SVA(p)
564* .. recompute SVA(q), SVA(p)
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 )**2.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* end of OK rotation
592 ELSE
593 notrot = notrot + 1
594*[RTD] SKIPPED = SKIPPED + 1
595 pskipped = pskipped + 1
596 ijblsk = ijblsk + 1
597 END IF
598 ELSE
599 notrot = notrot + 1
600 pskipped = pskipped + 1
601 ijblsk = ijblsk + 1
602 END IF
603*
604 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
605 \$ THEN
606 sva( p ) = aapp
607 notrot = 0
608 GO TO 2011
609 END IF
610 IF( ( i.LE.swband ) .AND.
611 \$ ( pskipped.GT.rowskip ) ) THEN
612 aapp = -aapp
613 notrot = 0
614 GO TO 2203
615 END IF
616*
617 2200 CONTINUE
618* end of the q-loop
619 2203 CONTINUE
620*
621 sva( p ) = aapp
622*
623 ELSE
624*
625 IF( aapp.EQ.zero )notrot = notrot +
626 \$ min( jgl+kbl-1, n ) - jgl + 1
627 IF( aapp.LT.zero )notrot = 0
628*
629 END IF
630*
631 2100 CONTINUE
632* end of the p-loop
633 2010 CONTINUE
634* end of the jbc-loop
635 2011 CONTINUE
636*2011 bailed out of the jbc-loop
637 DO 2012 p = igl, min( igl+kbl-1, n )
638 sva( p ) = abs( sva( p ) )
639 2012 CONTINUE
640***
641 2000 CONTINUE
642*2000 :: end of the ibr-loop
643*
644* .. update SVA(N)
645 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
646 \$ THEN
647 sva( n ) = scnrm2( m, a( 1, n ), 1 )
648 ELSE
649 t = zero
650 aapp = one
651 CALL classq( m, a( 1, n ), 1, t, aapp )
652 sva( n ) = t*sqrt( aapp )
653 END IF
654*
656*
657 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
658 \$ ( iswrot.LE.n ) ) )swband = i
659*
660 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( real( n ) )*
661 \$ tol ) .AND. ( real( n )*mxaapq*mxsinj.LT.tol ) ) THEN
662 GO TO 1994
663 END IF
664*
665 IF( notrot.GE.emptsw )GO TO 1994
666*
667 1993 CONTINUE
668* end i=1:NSWEEP loop
669*
670* #:( Reaching this point means that the procedure has not converged.
671 info = nsweep - 1
672 GO TO 1995
673*
674 1994 CONTINUE
675* #:) Reaching this point means numerical convergence after the i-th
676* sweep.
677*
678 info = 0
679* #:) INFO = 0 confirms successful iterations.
680 1995 CONTINUE
681*
682* Sort the vector SVA() of column norms.
683 DO 5991 p = 1, n - 1
684 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
685 IF( p.NE.q ) THEN
686 temp1 = sva( p )
687 sva( p ) = sva( q )
688 sva( q ) = temp1
689 aapq = d( p )
690 d( p ) = d( q )
691 d( q ) = aapq
692 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
693 IF( rsvec )CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
694 END IF
695 5991 CONTINUE
696*
697*
698 RETURN
699* ..
700* .. END OF CGSVJ1
701* ..
702 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgsvj1(jobv, m, n, n1, a, lda, d, sva, mv, v, ldv, eps, sfmin, tol, nsweep, work, lwork, info)
CGSVJ1 pre-processor for the routine cgesvj, applies Jacobi rotations targeting only particular pivot...
Definition cgsvj1.f:236
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
subroutine classq(n, x, incx, scale, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:124
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition crot.f:103
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81