Annotation of rpl/lapack/lapack/dgsvj1.f, revision 1.6
1.1 bertrand 1: SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
1.6 ! bertrand 2: $ EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
1.1 bertrand 3: *
1.6 ! bertrand 4: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 5: *
6: * -- Contributed by Zlatko Drmac of the University of Zagreb and --
7: * -- Kresimir Veselic of the Fernuniversitaet Hagen --
1.6 ! bertrand 8: * -- April 2011 --
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, * ),
1.6 ! bertrand 27: $ WORK( LWORK )
1.1 bertrand 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,
1.6 ! bertrand 165: $ TWO = 2.0D0 )
1.1 bertrand 166: * ..
167: * .. Local Scalars ..
168: DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
1.6 ! bertrand 169: $ BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
! 170: $ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
! 171: $ TEMP1, THETA, THSIGN
1.1 bertrand 172: INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
1.6 ! bertrand 173: $ ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
! 174: $ p, PSKIPPED, q, ROWSKIP, SWBAND
1.1 bertrand 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.
1.6 ! bertrand 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,
1.6 ! bertrand 336: $ q ), 1 )*D( p )*D( q ) / AAQQ )
! 337: $ / AAPP
1.1 bertrand 338: ELSE
339: CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
340: CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
1.6 ! bertrand 341: $ M, 1, WORK, LDA, IERR )
1.1 bertrand 342: AAPQ = DDOT( M, WORK, 1, A( 1, q ),
1.6 ! bertrand 343: $ 1 )*D( q ) / AAQQ
1.1 bertrand 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,
1.6 ! bertrand 353: $ q ), 1 )*D( p )*D( q ) / AAQQ )
! 354: $ / AAPP
1.1 bertrand 355: ELSE
356: CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
357: CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
1.6 ! bertrand 358: $ M, 1, WORK, LDA, IERR )
1.1 bertrand 359: AAPQ = DDOT( M, WORK, 1, A( 1, p ),
1.6 ! bertrand 360: $ 1 )*D( p ) / AAPP
1.1 bertrand 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
1.6 ! bertrand 378: THETA = -HALF*DABS(AQOAP-APOAQ) / AAPQ
1.1 bertrand 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,
1.6 ! bertrand 386: $ A( 1, q ), 1, FASTR )
1.1 bertrand 387: IF( RSVEC )CALL DROTM( MVL,
1.6 ! bertrand 388: $ V( 1, p ), 1,
! 389: $ V( 1, q ), 1,
! 390: $ FASTR )
1.1 bertrand 391: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6 ! bertrand 392: $ ONE+T*APOAQ*AAPQ ) )
1.1 bertrand 393: AAPP = AAPP*DSQRT( DMAX1( ZERO,
1.6 ! bertrand 394: $ ONE-T*AQOAP*AAPQ ) )
1.1 bertrand 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*
1.6 ! bertrand 403: $ DSQRT( ONE+THETA*THETA ) )
1.1 bertrand 404: CS = DSQRT( ONE / ( ONE+T*T ) )
405: SN = T*CS
406: MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
407: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6 ! bertrand 408: $ ONE+T*APOAQ*AAPQ ) )
1.4 bertrand 409: AAPP = AAPP*DSQRT( DMAX1( ZERO,
1.6 ! bertrand 410: $ ONE-T*AQOAP*AAPQ ) )
1.1 bertrand 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,
1.6 ! bertrand 422: $ A( 1, q ), 1,
! 423: $ FASTR )
1.1 bertrand 424: IF( RSVEC )CALL DROTM( MVL,
1.6 ! bertrand 425: $ V( 1, p ), 1, V( 1, q ),
! 426: $ 1, FASTR )
1.1 bertrand 427: ELSE
428: CALL DAXPY( M, -T*AQOAP,
1.6 ! bertrand 429: $ A( 1, q ), 1,
! 430: $ A( 1, p ), 1 )
1.1 bertrand 431: CALL DAXPY( M, CS*SN*APOAQ,
1.6 ! bertrand 432: $ A( 1, p ), 1,
! 433: $ A( 1, q ), 1 )
1.1 bertrand 434: IF( RSVEC ) THEN
435: CALL DAXPY( MVL, -T*AQOAP,
1.6 ! bertrand 436: $ V( 1, q ), 1,
! 437: $ V( 1, p ), 1 )
1.1 bertrand 438: CALL DAXPY( MVL,
1.6 ! bertrand 439: $ CS*SN*APOAQ,
! 440: $ V( 1, p ), 1,
! 441: $ V( 1, q ), 1 )
1.1 bertrand 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,
1.6 ! bertrand 449: $ A( 1, p ), 1,
! 450: $ A( 1, q ), 1 )
1.1 bertrand 451: CALL DAXPY( M, -CS*SN*AQOAP,
1.6 ! bertrand 452: $ A( 1, q ), 1,
! 453: $ A( 1, p ), 1 )
1.1 bertrand 454: IF( RSVEC ) THEN
455: CALL DAXPY( MVL, T*APOAQ,
1.6 ! bertrand 456: $ V( 1, p ), 1,
! 457: $ V( 1, q ), 1 )
1.1 bertrand 458: CALL DAXPY( MVL,
1.6 ! bertrand 459: $ -CS*SN*AQOAP,
! 460: $ V( 1, q ), 1,
! 461: $ V( 1, p ), 1 )
1.1 bertrand 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,
1.6 ! bertrand 468: $ A( 1, q ), 1,
! 469: $ A( 1, p ), 1 )
1.1 bertrand 470: CALL DAXPY( M, CS*SN*APOAQ,
1.6 ! bertrand 471: $ A( 1, p ), 1,
! 472: $ A( 1, q ), 1 )
1.1 bertrand 473: D( p ) = D( p )*CS
474: D( q ) = D( q ) / CS
475: IF( RSVEC ) THEN
476: CALL DAXPY( MVL,
1.6 ! bertrand 477: $ -T*AQOAP,
! 478: $ V( 1, q ), 1,
! 479: $ V( 1, p ), 1 )
1.1 bertrand 480: CALL DAXPY( MVL,
1.6 ! bertrand 481: $ CS*SN*APOAQ,
! 482: $ V( 1, p ), 1,
! 483: $ V( 1, q ), 1 )
1.1 bertrand 484: END IF
485: ELSE
486: CALL DAXPY( M, T*APOAQ,
1.6 ! bertrand 487: $ A( 1, p ), 1,
! 488: $ A( 1, q ), 1 )
1.1 bertrand 489: CALL DAXPY( M,
1.6 ! bertrand 490: $ -CS*SN*AQOAP,
! 491: $ A( 1, q ), 1,
! 492: $ A( 1, p ), 1 )
1.1 bertrand 493: D( p ) = D( p ) / CS
494: D( q ) = D( q )*CS
495: IF( RSVEC ) THEN
496: CALL DAXPY( MVL,
1.6 ! bertrand 497: $ T*APOAQ, V( 1, p ),
! 498: $ 1, V( 1, q ), 1 )
1.1 bertrand 499: CALL DAXPY( MVL,
1.6 ! bertrand 500: $ -CS*SN*AQOAP,
! 501: $ V( 1, q ), 1,
! 502: $ V( 1, p ), 1 )
1.1 bertrand 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,
1.6 ! bertrand 512: $ 1 )
1.1 bertrand 513: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
1.6 ! bertrand 514: $ M, 1, WORK, LDA, IERR )
1.1 bertrand 515: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
1.6 ! bertrand 516: $ M, 1, A( 1, q ), LDA,
! 517: $ IERR )
1.1 bertrand 518: TEMP1 = -AAPQ*D( p ) / D( q )
519: CALL DAXPY( M, TEMP1, WORK, 1,
1.6 ! bertrand 520: $ A( 1, q ), 1 )
1.1 bertrand 521: CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
1.6 ! bertrand 522: $ M, 1, A( 1, q ), LDA,
! 523: $ IERR )
1.1 bertrand 524: SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6 ! bertrand 525: $ ONE-AAPQ*AAPQ ) )
1.1 bertrand 526: MXSINJ = DMAX1( MXSINJ, SFMIN )
527: ELSE
528: CALL DCOPY( M, A( 1, q ), 1, WORK,
1.6 ! bertrand 529: $ 1 )
1.1 bertrand 530: CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
1.6 ! bertrand 531: $ M, 1, WORK, LDA, IERR )
1.1 bertrand 532: CALL DLASCL( 'G', 0, 0, AAPP, ONE,
1.6 ! bertrand 533: $ M, 1, A( 1, p ), LDA,
! 534: $ IERR )
1.1 bertrand 535: TEMP1 = -AAPQ*D( q ) / D( p )
536: CALL DAXPY( M, TEMP1, WORK, 1,
1.6 ! bertrand 537: $ A( 1, p ), 1 )
1.1 bertrand 538: CALL DLASCL( 'G', 0, 0, ONE, AAPP,
1.6 ! bertrand 539: $ M, 1, A( 1, p ), LDA,
! 540: $ IERR )
1.1 bertrand 541: SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
1.6 ! bertrand 542: $ ONE-AAPQ*AAPQ ) )
1.1 bertrand 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 )
1.6 ! bertrand 551: $ THEN
1.1 bertrand 552: IF( ( AAQQ.LT.ROOTBIG ) .AND.
1.6 ! bertrand 553: $ ( AAQQ.GT.ROOTSFMIN ) ) THEN
1.1 bertrand 554: SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
1.6 ! bertrand 555: $ D( q )
1.1 bertrand 556: ELSE
557: T = ZERO
1.4 bertrand 558: AAQQ = ONE
1.1 bertrand 559: CALL DLASSQ( M, A( 1, q ), 1, T,
1.6 ! bertrand 560: $ AAQQ )
1.1 bertrand 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.
1.6 ! bertrand 566: $ ( AAPP.GT.ROOTSFMIN ) ) THEN
1.1 bertrand 567: AAPP = DNRM2( M, A( 1, p ), 1 )*
1.6 ! bertrand 568: $ D( p )
1.1 bertrand 569: ELSE
570: T = ZERO
1.4 bertrand 571: AAPP = ONE
1.1 bertrand 572: CALL DLASSQ( M, A( 1, p ), 1, T,
1.6 ! bertrand 573: $ AAPP )
1.1 bertrand 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 ) )
1.6 ! bertrand 593: $ THEN
1.1 bertrand 594: SVA( p ) = AAPP
595: NOTROT = 0
596: GO TO 2011
597: END IF
598: IF( ( i.LE.SWBAND ) .AND.
1.6 ! bertrand 599: $ ( PSKIPPED.GT.ROWSKIP ) ) THEN
1.1 bertrand 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 +
1.6 ! bertrand 614: $ MIN0( jgl+KBL-1, N ) - jgl + 1
1.1 bertrand 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 ) )
1.6 ! bertrand 634: $ THEN
1.1 bertrand 635: SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
636: ELSE
637: T = ZERO
1.4 bertrand 638: AAPP = ONE
1.1 bertrand 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.
1.6 ! bertrand 646: $ ( ISWROT.LE.N ) ) )SWBAND = i
1.1 bertrand 647:
648: IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
1.6 ! bertrand 649: $ ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[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