![[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: *
1.4 ! bertrand 4: * -- LAPACK routine (version 3.3.0) --
1.1 bertrand 5: *
6: * -- Contributed by Zlatko Drmac of the University of Zagreb and --
7: * -- Kresimir Veselic of the Fernuniversitaet Hagen --
1.4 ! bertrand 8: * November 2010
1.1 bertrand 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
1.4 ! bertrand 208: ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
1.1 bertrand 209: INFO = -9
1.4 ! bertrand 210: ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
! 211: & ( APPLV.AND.( LDV.LT.MV ) ) ) THEN
1.1 bertrand 212: INFO = -11
213: ELSE IF( TOL.LE.EPS ) THEN
214: INFO = -14
215: ELSE IF( NSWEEP.LT.0 ) THEN
216: INFO = -15
217: ELSE IF( LWORK.LT.M ) THEN
218: INFO = -17
219: ELSE
220: INFO = 0
221: END IF
222: *
223: * #:(
224: IF( INFO.NE.0 ) THEN
225: CALL XERBLA( 'DGSVJ1', -INFO )
226: RETURN
227: END IF
228: *
229: IF( RSVEC ) THEN
230: MVL = N
231: ELSE IF( APPLV ) THEN
232: MVL = MV
233: END IF
234: RSVEC = RSVEC .OR. APPLV
235:
236: ROOTEPS = DSQRT( EPS )
237: ROOTSFMIN = DSQRT( SFMIN )
238: SMALL = SFMIN / EPS
239: BIG = ONE / SFMIN
240: ROOTBIG = ONE / ROOTSFMIN
241: LARGE = BIG / DSQRT( DBLE( M*N ) )
242: BIGTHETA = ONE / ROOTEPS
243: ROOTTOL = DSQRT( TOL )
244: *
245: * .. Initialize the right singular vector matrix ..
246: *
247: * RSVEC = LSAME( JOBV, 'Y' )
248: *
249: EMPTSW = N1*( N-N1 )
250: NOTROT = 0
251: FASTR( 1 ) = ZERO
252: *
253: * .. Row-cyclic pivot strategy with de Rijk's pivoting ..
254: *
255: KBL = MIN0( 8, N )
256: NBLR = N1 / KBL
257: IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
258:
259: * .. the tiling is nblr-by-nblc [tiles]
260:
261: NBLC = ( N-N1 ) / KBL
262: IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
263: BLSKIP = ( KBL**2 ) + 1
264: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
265:
266: ROWSKIP = MIN0( 5, KBL )
267: *[TP] ROWSKIP is a tuning parameter.
268: SWBAND = 0
269: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
270: * if SGESVJ is used as a computational routine in the preconditioned
271: * Jacobi SVD algorithm SGESVJ.
272: *
273: *
274: * | * * * [x] [x] [x]|
275: * | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
276: * | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
277: * |[x] [x] [x] * * * |
278: * |[x] [x] [x] * * * |
279: * |[x] [x] [x] * * * |
280: *
281: *
282: DO 1993 i = 1, NSWEEP
283: * .. go go go ...
284: *
285: MXAAPQ = ZERO
286: MXSINJ = ZERO
287: ISWROT = 0
288: *
289: NOTROT = 0
290: PSKIPPED = 0
291: *
292: DO 2000 ibr = 1, NBLR
293:
294: igl = ( ibr-1 )*KBL + 1
295: *
296: *
297: *........................................................
298: * ... go to the off diagonal blocks
299:
300: igl = ( ibr-1 )*KBL + 1
301:
302: DO 2010 jbc = 1, NBLC
303:
304: jgl = N1 + ( jbc-1 )*KBL + 1
305:
306: * doing the block at ( ibr, jbc )
307:
308: IJBLSK = 0
309: DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
310:
311: AAPP = SVA( p )
312:
313: IF( AAPP.GT.ZERO ) THEN
314:
315: PSKIPPED = 0
316:
317: DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
318: *
319: AAQQ = SVA( q )
320:
321: IF( AAQQ.GT.ZERO ) THEN
322: AAPP0 = AAPP
323: *
324: * .. M x 2 Jacobi SVD ..
325: *
326: * .. Safe Gram matrix computation ..
327: *
328: IF( AAQQ.GE.ONE ) THEN
329: IF( AAPP.GE.AAQQ ) THEN
330: ROTOK = ( SMALL*AAPP ).LE.AAQQ
331: ELSE
332: ROTOK = ( SMALL*AAQQ ).LE.AAPP
333: END IF
334: IF( AAPP.LT.( BIG / AAQQ ) ) THEN
335: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
336: + q ), 1 )*D( p )*D( q ) / AAQQ )
337: + / AAPP
338: ELSE
339: CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
340: CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
341: + M, 1, WORK, LDA, IERR )
342: AAPQ = DDOT( M, WORK, 1, A( 1, q ),
343: + 1 )*D( q ) / AAQQ
344: END IF
345: ELSE
346: IF( AAPP.GE.AAQQ ) THEN
347: ROTOK = AAPP.LE.( AAQQ / SMALL )
348: ELSE
349: ROTOK = AAQQ.LE.( AAPP / SMALL )
350: END IF
351: IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
352: AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
353: + q ), 1 )*D( p )*D( q ) / AAQQ )
354: + / AAPP
355: ELSE
356: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
357: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
358: + M, 1, WORK, LDA, IERR )
359: AAPQ = DDOT( M, WORK, 1, A( 1, p ),
360: + 1 )*D( p ) / AAPP
361: END IF
362: END IF
363:
364: MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
365:
366: * TO rotate or NOT to rotate, THAT is the question ...
367: *
368: IF( DABS( AAPQ ).GT.TOL ) THEN
369: NOTROT = 0
370: * ROTATED = ROTATED + 1
371: PSKIPPED = 0
372: ISWROT = ISWROT + 1
373: *
374: IF( ROTOK ) THEN
375: *
376: AQOAP = AAQQ / AAPP
377: APOAQ = AAPP / AAQQ
378: THETA = -HALF*DABS( AQOAP-APOAQ ) /
379: + AAPQ
380: IF( AAQQ.GT.AAPP0 )THETA = -THETA
381:
382: IF( DABS( THETA ).GT.BIGTHETA ) THEN
383: T = HALF / THETA
384: FASTR( 3 ) = T*D( p ) / D( q )
385: FASTR( 4 ) = -T*D( q ) / D( p )
386: CALL DROTM( M, A( 1, p ), 1,
387: + A( 1, q ), 1, FASTR )
388: IF( RSVEC )CALL DROTM( MVL,
389: + V( 1, p ), 1,
390: + V( 1, q ), 1,
391: + FASTR )
392: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
393: + ONE+T*APOAQ*AAPQ ) )
394: AAPP = AAPP*DSQRT( DMAX1( ZERO,
395: + ONE-T*AQOAP*AAPQ ) )
396: MXSINJ = DMAX1( MXSINJ, DABS( T ) )
397: ELSE
398: *
399: * .. choose correct signum for THETA and rotate
400: *
401: THSIGN = -DSIGN( ONE, AAPQ )
402: IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
403: T = ONE / ( THETA+THSIGN*
404: + DSQRT( ONE+THETA*THETA ) )
405: CS = DSQRT( ONE / ( ONE+T*T ) )
406: SN = T*CS
407: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
408: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
409: + ONE+T*APOAQ*AAPQ ) )
1.4 ! bertrand 410: AAPP = AAPP*DSQRT( DMAX1( ZERO,
! 411: + ONE-T*AQOAP*AAPQ ) )
1.1 bertrand 412:
413: APOAQ = D( p ) / D( q )
414: AQOAP = D( q ) / D( p )
415: IF( D( p ).GE.ONE ) THEN
416: *
417: IF( D( q ).GE.ONE ) THEN
418: FASTR( 3 ) = T*APOAQ
419: FASTR( 4 ) = -T*AQOAP
420: D( p ) = D( p )*CS
421: D( q ) = D( q )*CS
422: CALL DROTM( M, A( 1, p ), 1,
423: + A( 1, q ), 1,
424: + FASTR )
425: IF( RSVEC )CALL DROTM( MVL,
426: + V( 1, p ), 1, V( 1, q ),
427: + 1, FASTR )
428: ELSE
429: CALL DAXPY( M, -T*AQOAP,
430: + A( 1, q ), 1,
431: + A( 1, p ), 1 )
432: CALL DAXPY( M, CS*SN*APOAQ,
433: + A( 1, p ), 1,
434: + A( 1, q ), 1 )
435: IF( RSVEC ) THEN
436: CALL DAXPY( MVL, -T*AQOAP,
437: + V( 1, q ), 1,
438: + V( 1, p ), 1 )
439: CALL DAXPY( MVL,
440: + CS*SN*APOAQ,
441: + V( 1, p ), 1,
442: + V( 1, q ), 1 )
443: END IF
444: D( p ) = D( p )*CS
445: D( q ) = D( q ) / CS
446: END IF
447: ELSE
448: IF( D( q ).GE.ONE ) THEN
449: CALL DAXPY( M, T*APOAQ,
450: + A( 1, p ), 1,
451: + A( 1, q ), 1 )
452: CALL DAXPY( M, -CS*SN*AQOAP,
453: + A( 1, q ), 1,
454: + A( 1, p ), 1 )
455: IF( RSVEC ) THEN
456: CALL DAXPY( MVL, T*APOAQ,
457: + V( 1, p ), 1,
458: + V( 1, q ), 1 )
459: CALL DAXPY( MVL,
460: + -CS*SN*AQOAP,
461: + V( 1, q ), 1,
462: + V( 1, p ), 1 )
463: END IF
464: D( p ) = D( p ) / CS
465: D( q ) = D( q )*CS
466: ELSE
467: IF( D( p ).GE.D( q ) ) THEN
468: CALL DAXPY( M, -T*AQOAP,
469: + A( 1, q ), 1,
470: + A( 1, p ), 1 )
471: CALL DAXPY( M, CS*SN*APOAQ,
472: + A( 1, p ), 1,
473: + A( 1, q ), 1 )
474: D( p ) = D( p )*CS
475: D( q ) = D( q ) / CS
476: IF( RSVEC ) THEN
477: CALL DAXPY( MVL,
478: + -T*AQOAP,
479: + V( 1, q ), 1,
480: + V( 1, p ), 1 )
481: CALL DAXPY( MVL,
482: + CS*SN*APOAQ,
483: + V( 1, p ), 1,
484: + V( 1, q ), 1 )
485: END IF
486: ELSE
487: CALL DAXPY( M, T*APOAQ,
488: + A( 1, p ), 1,
489: + A( 1, q ), 1 )
490: CALL DAXPY( M,
491: + -CS*SN*AQOAP,
492: + A( 1, q ), 1,
493: + A( 1, p ), 1 )
494: D( p ) = D( p ) / CS
495: D( q ) = D( q )*CS
496: IF( RSVEC ) THEN
497: CALL DAXPY( MVL,
498: + T*APOAQ, V( 1, p ),
499: + 1, V( 1, q ), 1 )
500: CALL DAXPY( MVL,
501: + -CS*SN*AQOAP,
502: + V( 1, q ), 1,
503: + V( 1, p ), 1 )
504: END IF
505: END IF
506: END IF
507: END IF
508: END IF
509:
510: ELSE
511: IF( AAPP.GT.AAQQ ) THEN
512: CALL DCOPY( M, A( 1, p ), 1, WORK,
513: + 1 )
514: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
515: + M, 1, WORK, LDA, IERR )
516: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
517: + M, 1, A( 1, q ), LDA,
518: + IERR )
519: TEMP1 = -AAPQ*D( p ) / D( q )
520: CALL DAXPY( M, TEMP1, WORK, 1,
521: + A( 1, q ), 1 )
522: CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
523: + M, 1, A( 1, q ), LDA,
524: + IERR )
525: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
526: + ONE-AAPQ*AAPQ ) )
527: MXSINJ = DMAX1( MXSINJ, SFMIN )
528: ELSE
529: CALL DCOPY( M, A( 1, q ), 1, WORK,
530: + 1 )
531: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
532: + M, 1, WORK, LDA, IERR )
533: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
534: + M, 1, A( 1, p ), LDA,
535: + IERR )
536: TEMP1 = -AAPQ*D( q ) / D( p )
537: CALL DAXPY( M, TEMP1, WORK, 1,
538: + A( 1, p ), 1 )
539: CALL DLASCL( 'G', 0, 0, ONE, AAPP,
540: + M, 1, A( 1, p ), LDA,
541: + IERR )
542: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
543: + ONE-AAPQ*AAPQ ) )
544: MXSINJ = DMAX1( MXSINJ, SFMIN )
545: END IF
546: END IF
547: * END IF ROTOK THEN ... ELSE
548: *
549: * In the case of cancellation in updating SVA(q)
550: * .. recompute SVA(q)
551: IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
552: + THEN
553: IF( ( AAQQ.LT.ROOTBIG ) .AND.
554: + ( AAQQ.GT.ROOTSFMIN ) ) THEN
555: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
556: + D( q )
557: ELSE
558: T = ZERO
1.4 ! bertrand 559: AAQQ = ONE
1.1 bertrand 560: CALL DLASSQ( M, A( 1, q ), 1, T,
561: + AAQQ )
562: SVA( q ) = T*DSQRT( AAQQ )*D( q )
563: END IF
564: END IF
565: IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
566: IF( ( AAPP.LT.ROOTBIG ) .AND.
567: + ( AAPP.GT.ROOTSFMIN ) ) THEN
568: AAPP = DNRM2( M, A( 1, p ), 1 )*
569: + D( p )
570: ELSE
571: T = ZERO
1.4 ! bertrand 572: AAPP = ONE
1.1 bertrand 573: CALL DLASSQ( M, A( 1, p ), 1, T,
574: + AAPP )
575: AAPP = T*DSQRT( AAPP )*D( p )
576: END IF
577: SVA( p ) = AAPP
578: END IF
579: * end of OK rotation
580: ELSE
581: NOTROT = NOTROT + 1
582: * SKIPPED = SKIPPED + 1
583: PSKIPPED = PSKIPPED + 1
584: IJBLSK = IJBLSK + 1
585: END IF
586: ELSE
587: NOTROT = NOTROT + 1
588: PSKIPPED = PSKIPPED + 1
589: IJBLSK = IJBLSK + 1
590: END IF
591:
592: * IF ( NOTROT .GE. EMPTSW ) GO TO 2011
593: IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
594: + THEN
595: SVA( p ) = AAPP
596: NOTROT = 0
597: GO TO 2011
598: END IF
599: IF( ( i.LE.SWBAND ) .AND.
600: + ( PSKIPPED.GT.ROWSKIP ) ) THEN
601: AAPP = -AAPP
602: NOTROT = 0
603: GO TO 2203
604: END IF
605:
606: *
607: 2200 CONTINUE
608: * end of the q-loop
609: 2203 CONTINUE
610:
611: SVA( p ) = AAPP
612: *
613: ELSE
614: IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
615: + MIN0( jgl+KBL-1, N ) - jgl + 1
616: IF( AAPP.LT.ZERO )NOTROT = 0
617: *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
618: END IF
619:
620: 2100 CONTINUE
621: * end of the p-loop
622: 2010 CONTINUE
623: * end of the jbc-loop
624: 2011 CONTINUE
625: *2011 bailed out of the jbc-loop
626: DO 2012 p = igl, MIN0( igl+KBL-1, N )
627: SVA( p ) = DABS( SVA( p ) )
628: 2012 CONTINUE
629: *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
630: 2000 CONTINUE
631: *2000 :: end of the ibr-loop
632: *
633: * .. update SVA(N)
634: IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
635: + THEN
636: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
637: ELSE
638: T = ZERO
1.4 ! bertrand 639: AAPP = ONE
1.1 bertrand 640: CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
641: SVA( N ) = T*DSQRT( AAPP )*D( N )
642: END IF
643: *
644: * Additional steering devices
645: *
646: IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
647: + ( ISWROT.LE.N ) ) )SWBAND = i
648:
649: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
650: + ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
651: GO TO 1994
652: END IF
653:
654: *
655: IF( NOTROT.GE.EMPTSW )GO TO 1994
656:
657: 1993 CONTINUE
658: * end i=1:NSWEEP loop
659: * #:) Reaching this point means that the procedure has completed the given
660: * number of sweeps.
661: INFO = NSWEEP - 1
662: GO TO 1995
663: 1994 CONTINUE
664: * #:) Reaching this point means that during the i-th sweep all pivots were
665: * below the given threshold, causing early exit.
666:
667: INFO = 0
668: * #:) INFO = 0 confirms successful iterations.
669: 1995 CONTINUE
670: *
671: * Sort the vector D
672: *
673: DO 5991 p = 1, N - 1
674: q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
675: IF( p.NE.q ) THEN
676: TEMP1 = SVA( p )
677: SVA( p ) = SVA( q )
678: SVA( q ) = TEMP1
679: TEMP1 = D( p )
680: D( p ) = D( q )
681: D( q ) = TEMP1
682: CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
683: IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
684: END IF
685: 5991 CONTINUE
686: *
687: RETURN
688: * ..
689: * .. END OF DGSVJ1
690: * ..
691: END