![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgsvj1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dgsvj1.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
2: + EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2.2) --
5: *
6: * -- Contributed by Zlatko Drmac of the University of Zagreb and --
7: * -- Kresimir Veselic of the Fernuniversitaet Hagen --
8: * -- June 2010 --
9: *
10: * -- LAPACK is a software package provided by Univ. of Tennessee, --
11: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
12: *
13: * This routine is also part of SIGMA (version 1.23, October 23. 2008.)
14: * SIGMA is a library of algorithms for highly accurate algorithms for
15: * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the
16: * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
17: *
18: IMPLICIT NONE
19: * ..
20: * .. Scalar Arguments ..
21: DOUBLE PRECISION EPS, SFMIN, TOL
22: INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
23: CHARACTER*1 JOBV
24: * ..
25: * .. Array Arguments ..
26: DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
27: + WORK( LWORK )
28: * ..
29: *
30: * Purpose
31: * =======
32: *
33: * DGSVJ1 is called from SGESVJ as a pre-processor and that is its main
34: * purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
35: * it targets only particular pivots and it does not check convergence
36: * (stopping criterion). Few tunning parameters (marked by [TP]) are
37: * available for the implementer.
38: *
39: * Further Details
40: * ~~~~~~~~~~~~~~~
41: * DGSVJ1 applies few sweeps of Jacobi rotations in the column space of
42: * the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
43: * off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
44: * block-entries (tiles) of the (1,2) off-diagonal block are marked by the
45: * [x]'s in the following scheme:
46: *
47: * | * * * [x] [x] [x]|
48: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
49: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
50: * |[x] [x] [x] * * * |
51: * |[x] [x] [x] * * * |
52: * |[x] [x] [x] * * * |
53: *
54: * In terms of the columns of A, the first N1 columns are rotated 'against'
55: * the remaining N-N1 columns, trying to increase the angle between the
56: * corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
57: * tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
58: * The number of sweeps is given in NSWEEP and the orthogonality threshold
59: * is given in TOL.
60: *
61: * Contributors
62: * ~~~~~~~~~~~~
63: * Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
64: *
65: * Arguments
66: * =========
67: *
68: * JOBV (input) CHARACTER*1
69: * Specifies whether the output from this procedure is used
70: * to compute the matrix V:
71: * = 'V': the product of the Jacobi rotations is accumulated
72: * by postmulyiplying the N-by-N array V.
73: * (See the description of V.)
74: * = 'A': the product of the Jacobi rotations is accumulated
75: * by postmulyiplying the MV-by-N array V.
76: * (See the descriptions of MV and V.)
77: * = 'N': the Jacobi rotations are not accumulated.
78: *
79: * M (input) INTEGER
80: * The number of rows of the input matrix A. M >= 0.
81: *
82: * N (input) INTEGER
83: * The number of columns of the input matrix A.
84: * M >= N >= 0.
85: *
86: * N1 (input) INTEGER
87: * N1 specifies the 2 x 2 block partition, the first N1 columns are
88: * rotated 'against' the remaining N-N1 columns of A.
89: *
90: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
91: * On entry, M-by-N matrix A, such that A*diag(D) represents
92: * the input matrix.
93: * On exit,
94: * A_onexit * D_onexit represents the input matrix A*diag(D)
95: * post-multiplied by a sequence of Jacobi rotations, where the
96: * rotation threshold and the total number of sweeps are given in
97: * TOL and NSWEEP, respectively.
98: * (See the descriptions of N1, D, TOL and NSWEEP.)
99: *
100: * LDA (input) INTEGER
101: * The leading dimension of the array A. LDA >= max(1,M).
102: *
103: * D (input/workspace/output) DOUBLE PRECISION array, dimension (N)
104: * The array D accumulates the scaling factors from the fast scaled
105: * Jacobi rotations.
106: * On entry, A*diag(D) represents the input matrix.
107: * On exit, A_onexit*diag(D_onexit) represents the input matrix
108: * post-multiplied by a sequence of Jacobi rotations, where the
109: * rotation threshold and the total number of sweeps are given in
110: * TOL and NSWEEP, respectively.
111: * (See the descriptions of N1, A, TOL and NSWEEP.)
112: *
113: * SVA (input/workspace/output) DOUBLE PRECISION array, dimension (N)
114: * On entry, SVA contains the Euclidean norms of the columns of
115: * the matrix A*diag(D).
116: * On exit, SVA contains the Euclidean norms of the columns of
117: * the matrix onexit*diag(D_onexit).
118: *
119: * MV (input) 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: *
124: * V (input/output) DOUBLE PRECISION array, dimension (LDV,N)
125: * If JOBV .EQ. 'V' then N rows of V are post-multipled by a
126: * sequence of Jacobi rotations.
127: * If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
128: * sequence of Jacobi rotations.
129: * If JOBV = 'N', then V is not referenced.
130: *
131: * LDV (input) INTEGER
132: * The leading dimension of the array V, LDV >= 1.
133: * If JOBV = 'V', LDV .GE. N.
134: * If JOBV = 'A', LDV .GE. MV.
135: *
136: * EPS (input) DOUBLE PRECISION
137: * EPS = DLAMCH('Epsilon')
138: *
139: * SFMIN (input) DOUBLE PRECISION
140: * SFMIN = DLAMCH('Safe Minimum')
141: *
142: * TOL (input) DOUBLE PRECISION
143: * TOL is the threshold for Jacobi rotations. For a pair
144: * A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
145: * applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
146: *
147: * NSWEEP (input) INTEGER
148: * NSWEEP is the number of sweeps of Jacobi rotations to be
149: * performed.
150: *
151: * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
152: *
153: * LWORK (input) INTEGER
154: * LWORK is the dimension of WORK. LWORK .GE. M.
155: *
156: * INFO (output) INTEGER
157: * = 0 : successful exit.
158: * < 0 : if INFO = -i, then the i-th argument had an illegal value
159: *
160: * =====================================================================
161: *
162: * .. Local Parameters ..
163: DOUBLE PRECISION ZERO, HALF, ONE, TWO
164: PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
165: + TWO = 2.0D0 )
166: * ..
167: * .. Local Scalars ..
168: DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
169: + BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
170: + ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
171: + TEMP1, THETA, THSIGN
172: INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
173: + ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
174: + p, PSKIPPED, q, ROWSKIP, SWBAND
175: LOGICAL APPLV, ROTOK, RSVEC
176: * ..
177: * .. Local Arrays ..
178: DOUBLE PRECISION FASTR( 5 )
179: * ..
180: * .. Intrinsic Functions ..
181: INTRINSIC DABS, DMAX1, DBLE, MIN0, DSIGN, DSQRT
182: * ..
183: * .. External Functions ..
184: DOUBLE PRECISION DDOT, DNRM2
185: INTEGER IDAMAX
186: LOGICAL LSAME
187: EXTERNAL IDAMAX, LSAME, DDOT, DNRM2
188: * ..
189: * .. External Subroutines ..
190: EXTERNAL DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP
191: * ..
192: * .. Executable Statements ..
193: *
194: * Test the input parameters.
195: *
196: APPLV = LSAME( JOBV, 'A' )
197: RSVEC = LSAME( JOBV, 'V' )
198: IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
199: INFO = -1
200: ELSE IF( M.LT.0 ) THEN
201: INFO = -2
202: ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
203: INFO = -3
204: ELSE IF( N1.LT.0 ) THEN
205: INFO = -4
206: ELSE IF( LDA.LT.M ) THEN
207: INFO = -6
208: ELSE IF( MV.LT.0 ) THEN
209: INFO = -9
210: ELSE IF( LDV.LT.M ) THEN
211: INFO = -11
212: ELSE IF( TOL.LE.EPS ) THEN
213: INFO = -14
214: ELSE IF( NSWEEP.LT.0 ) THEN
215: INFO = -15
216: ELSE IF( LWORK.LT.M ) THEN
217: INFO = -17
218: ELSE
219: INFO = 0
220: END IF
221: *
222: * #:(
223: IF( INFO.NE.0 ) THEN
224: CALL XERBLA( 'DGSVJ1', -INFO )
225: RETURN
226: END IF
227: *
228: IF( RSVEC ) THEN
229: MVL = N
230: ELSE IF( APPLV ) THEN
231: MVL = MV
232: END IF
233: RSVEC = RSVEC .OR. APPLV
234:
235: ROOTEPS = DSQRT( EPS )
236: ROOTSFMIN = DSQRT( SFMIN )
237: SMALL = SFMIN / EPS
238: BIG = ONE / SFMIN
239: ROOTBIG = ONE / ROOTSFMIN
240: LARGE = BIG / DSQRT( DBLE( M*N ) )
241: BIGTHETA = ONE / ROOTEPS
242: ROOTTOL = DSQRT( TOL )
243: *
244: * .. Initialize the right singular vector matrix ..
245: *
246: * RSVEC = LSAME( JOBV, 'Y' )
247: *
248: EMPTSW = N1*( N-N1 )
249: NOTROT = 0
250: FASTR( 1 ) = ZERO
251: *
252: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
253: *
254: KBL = MIN0( 8, N )
255: NBLR = N1 / KBL
256: IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
257:
258: * .. the tiling is nblr-by-nblc [tiles]
259:
260: NBLC = ( N-N1 ) / KBL
261: IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
262: BLSKIP = ( KBL**2 ) + 1
263: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
264:
265: ROWSKIP = MIN0( 5, KBL )
266: *[TP] ROWSKIP is a tuning parameter.
267: SWBAND = 0
268: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
269: * if SGESVJ is used as a computational routine in the preconditioned
270: * Jacobi SVD algorithm SGESVJ.
271: *
272: *
273: * | * * * [x] [x] [x]|
274: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
275: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
276: * |[x] [x] [x] * * * |
277: * |[x] [x] [x] * * * |
278: * |[x] [x] [x] * * * |
279: *
280: *
281: DO 1993 i = 1, NSWEEP
282: * .. go go go ...
283: *
284: MXAAPQ = ZERO
285: MXSINJ = ZERO
286: ISWROT = 0
287: *
288: NOTROT = 0
289: PSKIPPED = 0
290: *
291: DO 2000 ibr = 1, NBLR
292:
293: igl = ( ibr-1 )*KBL + 1
294: *
295: *
296: *........................................................
297: * ... go to the off diagonal blocks
298:
299: igl = ( ibr-1 )*KBL + 1
300:
301: DO 2010 jbc = 1, NBLC
302:
303: jgl = N1 + ( jbc-1 )*KBL + 1
304:
305: * doing the block at ( ibr, jbc )
306:
307: IJBLSK = 0
308: DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
309:
310: AAPP = SVA( p )
311:
312: IF( AAPP.GT.ZERO ) THEN
313:
314: PSKIPPED = 0
315:
316: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
317: *
318: AAQQ = SVA( q )
319:
320: IF( AAQQ.GT.ZERO ) THEN
321: AAPP0 = AAPP
322: *
323: * .. M x 2 Jacobi SVD ..
324: *
325: * .. Safe Gram matrix computation ..
326: *
327: IF( AAQQ.GE.ONE ) THEN
328: IF( AAPP.GE.AAQQ ) THEN
329: ROTOK = ( SMALL*AAPP ).LE.AAQQ
330: ELSE
331: ROTOK = ( SMALL*AAQQ ).LE.AAPP
332: END IF
333: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
334: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
335: + q ), 1 )*D( p )*D( q ) / AAQQ )
336: + / AAPP
337: ELSE
338: CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
339: CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
340: + M, 1, WORK, LDA, IERR )
341: AAPQ = DDOT( M, WORK, 1, A( 1, q ),
342: + 1 )*D( q ) / AAQQ
343: END IF
344: ELSE
345: IF( AAPP.GE.AAQQ ) THEN
346: ROTOK = AAPP.LE.( AAQQ / SMALL )
347: ELSE
348: ROTOK = AAQQ.LE.( AAPP / SMALL )
349: END IF
350: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
351: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
352: + q ), 1 )*D( p )*D( q ) / AAQQ )
353: + / AAPP
354: ELSE
355: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
356: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
357: + M, 1, WORK, LDA, IERR )
358: AAPQ = DDOT( M, WORK, 1, A( 1, p ),
359: + 1 )*D( p ) / AAPP
360: END IF
361: END IF
362:
363: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
364:
365: * TO rotate or NOT to rotate, THAT is the question ...
366: *
367: IF( DABS( AAPQ ).GT.TOL ) THEN
368: NOTROT = 0
369: * ROTATED = ROTATED + 1
370: PSKIPPED = 0
371: ISWROT = ISWROT + 1
372: *
373: IF( ROTOK ) THEN
374: *
375: AQOAP = AAQQ / AAPP
376: APOAQ = AAPP / AAQQ
377: THETA = -HALF*DABS( AQOAP-APOAQ ) /
378: + AAPQ
379: IF( AAQQ.GT.AAPP0 )THETA = -THETA
380:
381: IF( DABS( THETA ).GT.BIGTHETA ) THEN
382: T = HALF / THETA
383: FASTR( 3 ) = T*D( p ) / D( q )
384: FASTR( 4 ) = -T*D( q ) / D( p )
385: CALL DROTM( M, A( 1, p ), 1,
386: + A( 1, q ), 1, FASTR )
387: IF( RSVEC )CALL DROTM( MVL,
388: + V( 1, p ), 1,
389: + V( 1, q ), 1,
390: + FASTR )
391: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
392: + ONE+T*APOAQ*AAPQ ) )
393: AAPP = AAPP*DSQRT( DMAX1( ZERO,
394: + ONE-T*AQOAP*AAPQ ) )
395: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
396: ELSE
397: *
398: * .. choose correct signum for THETA and rotate
399: *
400: THSIGN = -DSIGN( ONE, AAPQ )
401: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
402: T = ONE / ( THETA+THSIGN*
403: + DSQRT( ONE+THETA*THETA ) )
404: CS = DSQRT( ONE / ( ONE+T*T ) )
405: SN = T*CS
406: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
407: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
408: + ONE+T*APOAQ*AAPQ ) )
409: AAPP = AAPP*DSQRT( ONE-T*AQOAP*
410: + AAPQ )
411:
412: APOAQ = D( p ) / D( q )
413: AQOAP = D( q ) / D( p )
414: IF( D( p ).GE.ONE ) THEN
415: *
416: IF( D( q ).GE.ONE ) THEN
417: FASTR( 3 ) = T*APOAQ
418: FASTR( 4 ) = -T*AQOAP
419: D( p ) = D( p )*CS
420: D( q ) = D( q )*CS
421: CALL DROTM( M, A( 1, p ), 1,
422: + A( 1, q ), 1,
423: + FASTR )
424: IF( RSVEC )CALL DROTM( MVL,
425: + V( 1, p ), 1, V( 1, q ),
426: + 1, FASTR )
427: ELSE
428: CALL DAXPY( M, -T*AQOAP,
429: + A( 1, q ), 1,
430: + A( 1, p ), 1 )
431: CALL DAXPY( M, CS*SN*APOAQ,
432: + A( 1, p ), 1,
433: + A( 1, q ), 1 )
434: IF( RSVEC ) THEN
435: CALL DAXPY( MVL, -T*AQOAP,
436: + V( 1, q ), 1,
437: + V( 1, p ), 1 )
438: CALL DAXPY( MVL,
439: + CS*SN*APOAQ,
440: + V( 1, p ), 1,
441: + V( 1, q ), 1 )
442: END IF
443: D( p ) = D( p )*CS
444: D( q ) = D( q ) / CS
445: END IF
446: ELSE
447: IF( D( q ).GE.ONE ) THEN
448: CALL DAXPY( M, T*APOAQ,
449: + A( 1, p ), 1,
450: + A( 1, q ), 1 )
451: CALL DAXPY( M, -CS*SN*AQOAP,
452: + A( 1, q ), 1,
453: + A( 1, p ), 1 )
454: IF( RSVEC ) THEN
455: CALL DAXPY( MVL, T*APOAQ,
456: + V( 1, p ), 1,
457: + V( 1, q ), 1 )
458: CALL DAXPY( MVL,
459: + -CS*SN*AQOAP,
460: + V( 1, q ), 1,
461: + V( 1, p ), 1 )
462: END IF
463: D( p ) = D( p ) / CS
464: D( q ) = D( q )*CS
465: ELSE
466: IF( D( p ).GE.D( q ) ) THEN
467: CALL DAXPY( M, -T*AQOAP,
468: + A( 1, q ), 1,
469: + A( 1, p ), 1 )
470: CALL DAXPY( M, CS*SN*APOAQ,
471: + A( 1, p ), 1,
472: + A( 1, q ), 1 )
473: D( p ) = D( p )*CS
474: D( q ) = D( q ) / CS
475: IF( RSVEC ) THEN
476: CALL DAXPY( MVL,
477: + -T*AQOAP,
478: + V( 1, q ), 1,
479: + V( 1, p ), 1 )
480: CALL DAXPY( MVL,
481: + CS*SN*APOAQ,
482: + V( 1, p ), 1,
483: + V( 1, q ), 1 )
484: END IF
485: ELSE
486: CALL DAXPY( M, T*APOAQ,
487: + A( 1, p ), 1,
488: + A( 1, q ), 1 )
489: CALL DAXPY( M,
490: + -CS*SN*AQOAP,
491: + A( 1, q ), 1,
492: + A( 1, p ), 1 )
493: D( p ) = D( p ) / CS
494: D( q ) = D( q )*CS
495: IF( RSVEC ) THEN
496: CALL DAXPY( MVL,
497: + T*APOAQ, V( 1, p ),
498: + 1, V( 1, q ), 1 )
499: CALL DAXPY( MVL,
500: + -CS*SN*AQOAP,
501: + V( 1, q ), 1,
502: + V( 1, p ), 1 )
503: END IF
504: END IF
505: END IF
506: END IF
507: END IF
508:
509: ELSE
510: IF( AAPP.GT.AAQQ ) THEN
511: CALL DCOPY( M, A( 1, p ), 1, WORK,
512: + 1 )
513: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
514: + M, 1, WORK, LDA, IERR )
515: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
516: + M, 1, A( 1, q ), LDA,
517: + IERR )
518: TEMP1 = -AAPQ*D( p ) / D( q )
519: CALL DAXPY( M, TEMP1, WORK, 1,
520: + A( 1, q ), 1 )
521: CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
522: + M, 1, A( 1, q ), LDA,
523: + IERR )
524: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
525: + ONE-AAPQ*AAPQ ) )
526: MXSINJ = DMAX1( MXSINJ, SFMIN )
527: ELSE
528: CALL DCOPY( M, A( 1, q ), 1, WORK,
529: + 1 )
530: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
531: + M, 1, WORK, LDA, IERR )
532: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
533: + M, 1, A( 1, p ), LDA,
534: + IERR )
535: TEMP1 = -AAPQ*D( q ) / D( p )
536: CALL DAXPY( M, TEMP1, WORK, 1,
537: + A( 1, p ), 1 )
538: CALL DLASCL( 'G', 0, 0, ONE, AAPP,
539: + M, 1, A( 1, p ), LDA,
540: + IERR )
541: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
542: + ONE-AAPQ*AAPQ ) )
543: MXSINJ = DMAX1( MXSINJ, SFMIN )
544: END IF
545: END IF
546: * END IF ROTOK THEN ... ELSE
547: *
548: * In the case of cancellation in updating SVA(q)
549: * .. recompute SVA(q)
550: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
551: + THEN
552: IF( ( AAQQ.LT.ROOTBIG ) .AND.
553: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
554: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
555: + D( q )
556: ELSE
557: T = ZERO
558: AAQQ = ZERO
559: CALL DLASSQ( M, A( 1, q ), 1, T,
560: + AAQQ )
561: SVA( q ) = T*DSQRT( AAQQ )*D( q )
562: END IF
563: END IF
564: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
565: IF( ( AAPP.LT.ROOTBIG ) .AND.
566: + ( AAPP.GT.ROOTSFMIN ) ) THEN
567: AAPP = DNRM2( M, A( 1, p ), 1 )*
568: + D( p )
569: ELSE
570: T = ZERO
571: AAPP = ZERO
572: CALL DLASSQ( M, A( 1, p ), 1, T,
573: + AAPP )
574: AAPP = T*DSQRT( AAPP )*D( p )
575: END IF
576: SVA( p ) = AAPP
577: END IF
578: * end of OK rotation
579: ELSE
580: NOTROT = NOTROT + 1
581: * SKIPPED = SKIPPED + 1
582: PSKIPPED = PSKIPPED + 1
583: IJBLSK = IJBLSK + 1
584: END IF
585: ELSE
586: NOTROT = NOTROT + 1
587: PSKIPPED = PSKIPPED + 1
588: IJBLSK = IJBLSK + 1
589: END IF
590:
591: * IF ( NOTROT .GE. EMPTSW ) GO TO 2011
592: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
593: + THEN
594: SVA( p ) = AAPP
595: NOTROT = 0
596: GO TO 2011
597: END IF
598: IF( ( i.LE.SWBAND ) .AND.
599: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
600: AAPP = -AAPP
601: NOTROT = 0
602: GO TO 2203
603: END IF
604:
605: *
606: 2200 CONTINUE
607: * end of the q-loop
608: 2203 CONTINUE
609:
610: SVA( p ) = AAPP
611: *
612: ELSE
613: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
614: + MIN0( jgl+KBL-1, N ) - jgl + 1
615: IF( AAPP.LT.ZERO )NOTROT = 0
616: *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
617: END IF
618:
619: 2100 CONTINUE
620: * end of the p-loop
621: 2010 CONTINUE
622: * end of the jbc-loop
623: 2011 CONTINUE
624: *2011 bailed out of the jbc-loop
625: DO 2012 p = igl, MIN0( igl+KBL-1, N )
626: SVA( p ) = DABS( SVA( p ) )
627: 2012 CONTINUE
628: *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
629: 2000 CONTINUE
630: *2000 :: end of the ibr-loop
631: *
632: * .. update SVA(N)
633: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
634: + THEN
635: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
636: ELSE
637: T = ZERO
638: AAPP = ZERO
639: CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
640: SVA( N ) = T*DSQRT( AAPP )*D( N )
641: END IF
642: *
643: * Additional steering devices
644: *
645: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
646: + ( ISWROT.LE.N ) ) )SWBAND = i
647:
648: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
649: + ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
650: GO TO 1994
651: END IF
652:
653: *
654: IF( NOTROT.GE.EMPTSW )GO TO 1994
655:
656: 1993 CONTINUE
657: * end i=1:NSWEEP loop
658: * #:) Reaching this point means that the procedure has completed the given
659: * number of sweeps.
660: INFO = NSWEEP - 1
661: GO TO 1995
662: 1994 CONTINUE
663: * #:) Reaching this point means that during the i-th sweep all pivots were
664: * below the given threshold, causing early exit.
665:
666: INFO = 0
667: * #:) INFO = 0 confirms successful iterations.
668: 1995 CONTINUE
669: *
670: * Sort the vector D
671: *
672: DO 5991 p = 1, N - 1
673: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
674: IF( p.NE.q ) THEN
675: TEMP1 = SVA( p )
676: SVA( p ) = SVA( q )
677: SVA( q ) = TEMP1
678: TEMP1 = D( p )
679: D( p ) = D( q )
680: D( q ) = TEMP1
681: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
682: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
683: END IF
684: 5991 CONTINUE
685: *
686: RETURN
687: * ..
688: * .. END OF DGSVJ1
689: * ..
690: END