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