001:       SUBROUTINE DGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
002:      +                   EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2.1)                                    --
005: *
006: *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
007: *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
008: *  -- April 2009                                                      --
009: *
010: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
011: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
012: *
013: * This routine is also part of SIGMA (version 1.23, October 23. 2008.)
014: * SIGMA is a library of algorithms for highly accurate algorithms for
015: * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the
016: * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.
017: *
018:       IMPLICIT           NONE
019: *     ..
020: *     .. Scalar Arguments ..
021:       DOUBLE PRECISION   EPS, SFMIN, TOL
022:       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
023:       CHARACTER*1        JOBV
024: *     ..
025: *     .. Array Arguments ..
026:       DOUBLE PRECISION   A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
027:      +                   WORK( LWORK )
028: *     ..
029: *
030: *  Purpose
031: *  =======
032: *
033: *  DGSVJ1 is called from SGESVJ as a pre-processor and that is its main
034: *  purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
035: *  it targets only particular pivots and it does not check convergence
036: *  (stopping criterion). Few tunning parameters (marked by [TP]) are
037: *  available for the implementer.
038: *
039: *  Further Details
040: *  ~~~~~~~~~~~~~~~
041: *  DGSVJ1 applies few sweeps of Jacobi rotations in the column space of
042: *  the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
043: *  off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
044: *  block-entries (tiles) of the (1,2) off-diagonal block are marked by the
045: *  [x]'s in the following scheme:
046: *
047: *     | *   *   * [x] [x] [x]|
048: *     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
049: *     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
050: *     |[x] [x] [x] *   *   * |
051: *     |[x] [x] [x] *   *   * |
052: *     |[x] [x] [x] *   *   * |
053: *
054: *  In terms of the columns of A, the first N1 columns are rotated 'against'
055: *  the remaining N-N1 columns, trying to increase the angle between the
056: *  corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
057: *  tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
058: *  The number of sweeps is given in NSWEEP and the orthogonality threshold
059: *  is given in TOL.
060: *
061: *  Contributors
062: *  ~~~~~~~~~~~~
063: *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
064: *
065: *  Arguments
066: *  =========
067: *
068: *  JOBV    (input) CHARACTER*1
069: *          Specifies whether the output from this procedure is used
070: *          to compute the matrix V:
071: *          = 'V': the product of the Jacobi rotations is accumulated
072: *                 by postmulyiplying the N-by-N array V.
073: *                (See the description of V.)
074: *          = 'A': the product of the Jacobi rotations is accumulated
075: *                 by postmulyiplying the MV-by-N array V.
076: *                (See the descriptions of MV and V.)
077: *          = 'N': the Jacobi rotations are not accumulated.
078: *
079: *  M       (input) INTEGER
080: *          The number of rows of the input matrix A.  M >= 0.
081: *
082: *  N       (input) INTEGER
083: *          The number of columns of the input matrix A.
084: *          M >= N >= 0.
085: *
086: *  N1      (input) INTEGER
087: *          N1 specifies the 2 x 2 block partition, the first N1 columns are
088: *          rotated 'against' the remaining N-N1 columns of A.
089: *
090: *  A       (input/output) REAL array, dimension (LDA,N)
091: *          On entry, M-by-N matrix A, such that A*diag(D) represents
092: *          the input matrix.
093: *          On exit,
094: *          A_onexit * D_onexit represents the input matrix A*diag(D)
095: *          post-multiplied by a sequence of Jacobi rotations, where the
096: *          rotation threshold and the total number of sweeps are given in
097: *          TOL and NSWEEP, respectively.
098: *          (See the descriptions of N1, D, TOL and NSWEEP.)
099: *
100: *  LDA     (input) INTEGER
101: *          The leading dimension of the array A.  LDA >= max(1,M).
102: *
103: *  D       (input/workspace/output) REAL 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) REAL 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) REAL 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) INTEGER
137: *          EPS = SLAMCH('Epsilon')
138: *
139: *  SFMIN   (input) INTEGER
140: *          SFMIN = SLAMCH('Safe Minimum')
141: *
142: *  TOL     (input) REAL
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) REAL 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
691: