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