001:       SUBROUTINE DGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
002:      +                   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: *     .. Scalar Arguments ..
020:       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
021:       DOUBLE PRECISION   EPS, SFMIN, TOL
022:       CHARACTER*1        JOBV
023: *     ..
024: *     .. Array Arguments ..
025:       DOUBLE PRECISION   A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
026:      +                   WORK( LWORK )
027: *     ..
028: *
029: *  Purpose
030: *  =======
031: *
032: *  DGSVJ0 is called from DGESVJ as a pre-processor and that is its main
033: *  purpose. It applies Jacobi rotations in the same way as DGESVJ does, but
034: *  it does not check convergence (stopping criterion). Few tuning
035: *  parameters (marked by [TP]) are available for the implementer.
036: *
037: *  Further Details
038: *  ~~~~~~~~~~~~~~~
039: *  DGSVJ0 is used just to enable SGESVJ to call a simplified version of
040: *  itself to work on a submatrix of the original matrix.
041: *
042: *  Contributors
043: *  ~~~~~~~~~~~~
044: *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
045: *
046: *  Bugs, Examples and Comments
047: *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
048: *  Please report all bugs and send interesting test examples and comments to
049: *  drmac@math.hr. Thank you.
050: *
051: *  Arguments
052: *  =========
053: *
054: *  JOBV    (input) CHARACTER*1
055: *          Specifies whether the output from this procedure is used
056: *          to compute the matrix V:
057: *          = 'V': the product of the Jacobi rotations is accumulated
058: *                 by postmulyiplying the N-by-N array V.
059: *                (See the description of V.)
060: *          = 'A': the product of the Jacobi rotations is accumulated
061: *                 by postmulyiplying the MV-by-N array V.
062: *                (See the descriptions of MV and V.)
063: *          = 'N': the Jacobi rotations are not accumulated.
064: *
065: *  M       (input) INTEGER
066: *          The number of rows of the input matrix A.  M >= 0.
067: *
068: *  N       (input) INTEGER
069: *          The number of columns of the input matrix A.
070: *          M >= N >= 0.
071: *
072: *  A       (input/output) REAL array, dimension (LDA,N)
073: *          On entry, M-by-N matrix A, such that A*diag(D) represents
074: *          the input matrix.
075: *          On exit,
076: *          A_onexit * D_onexit represents the input matrix A*diag(D)
077: *          post-multiplied by a sequence of Jacobi rotations, where the
078: *          rotation threshold and the total number of sweeps are given in
079: *          TOL and NSWEEP, respectively.
080: *          (See the descriptions of D, TOL and NSWEEP.)
081: *
082: *  LDA     (input) INTEGER
083: *          The leading dimension of the array A.  LDA >= max(1,M).
084: *
085: *  D       (input/workspace/output) REAL array, dimension (N)
086: *          The array D accumulates the scaling factors from the fast scaled
087: *          Jacobi rotations.
088: *          On entry, A*diag(D) represents the input matrix.
089: *          On exit, A_onexit*diag(D_onexit) represents the input matrix
090: *          post-multiplied by a sequence of Jacobi rotations, where the
091: *          rotation threshold and the total number of sweeps are given in
092: *          TOL and NSWEEP, respectively.
093: *          (See the descriptions of A, TOL and NSWEEP.)
094: *
095: *  SVA     (input/workspace/output) REAL array, dimension (N)
096: *          On entry, SVA contains the Euclidean norms of the columns of
097: *          the matrix A*diag(D).
098: *          On exit, SVA contains the Euclidean norms of the columns of
099: *          the matrix onexit*diag(D_onexit).
100: *
101: *  MV      (input) INTEGER
102: *          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
103: *                           sequence of Jacobi rotations.
104: *          If JOBV = 'N',   then MV is not referenced.
105: *
106: *  V       (input/output) REAL array, dimension (LDV,N)
107: *          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
108: *                           sequence of Jacobi rotations.
109: *          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
110: *                           sequence of Jacobi rotations.
111: *          If JOBV = 'N',   then V is not referenced.
112: *
113: *  LDV     (input) INTEGER
114: *          The leading dimension of the array V,  LDV >= 1.
115: *          If JOBV = 'V', LDV .GE. N.
116: *          If JOBV = 'A', LDV .GE. MV.
117: *
118: *  EPS     (input) INTEGER
119: *          EPS = SLAMCH('Epsilon')
120: *
121: *  SFMIN   (input) INTEGER
122: *          SFMIN = SLAMCH('Safe Minimum')
123: *
124: *  TOL     (input) REAL
125: *          TOL is the threshold for Jacobi rotations. For a pair
126: *          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
127: *          applied only if DABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
128: *
129: *  NSWEEP  (input) INTEGER
130: *          NSWEEP is the number of sweeps of Jacobi rotations to be
131: *          performed.
132: *
133: *  WORK    (workspace) REAL array, dimension LWORK.
134: *
135: *  LWORK   (input) INTEGER
136: *          LWORK is the dimension of WORK. LWORK .GE. M.
137: *
138: *  INFO    (output) INTEGER
139: *          = 0 : successful exit.
140: *          < 0 : if INFO = -i, then the i-th argument had an illegal value
141: *
142: *  =====================================================================
143: *
144: *     .. Local Parameters ..
145:       DOUBLE PRECISION   ZERO, HALF, ONE, TWO
146:       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
147:      +                   TWO = 2.0D0 )
148: *     ..
149: *     .. Local Scalars ..
150:       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
151:      +                   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
152:      +                   ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
153:      +                   THSIGN
154:       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
155:      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
156:      +                   NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
157:       LOGICAL            APPLV, ROTOK, RSVEC
158: *     ..
159: *     .. Local Arrays ..
160:       DOUBLE PRECISION   FASTR( 5 )
161: *     ..
162: *     .. Intrinsic Functions ..
163:       INTRINSIC          DABS, DMAX1, DBLE, MIN0, DSIGN, DSQRT
164: *     ..
165: *     .. External Functions ..
166:       DOUBLE PRECISION   DDOT, DNRM2
167:       INTEGER            IDAMAX
168:       LOGICAL            LSAME
169:       EXTERNAL           IDAMAX, LSAME, DDOT, DNRM2
170: *     ..
171: *     .. External Subroutines ..
172:       EXTERNAL           DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP
173: *     ..
174: *     .. Executable Statements ..
175: *
176:       APPLV = LSAME( JOBV, 'A' )
177:       RSVEC = LSAME( JOBV, 'V' )
178:       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
179:          INFO = -1
180:       ELSE IF( M.LT.0 ) THEN
181:          INFO = -2
182:       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
183:          INFO = -3
184:       ELSE IF( LDA.LT.M ) THEN
185:          INFO = -5
186:       ELSE IF( MV.LT.0 ) THEN
187:          INFO = -8
188:       ELSE IF( LDV.LT.M ) THEN
189:          INFO = -10
190:       ELSE IF( TOL.LE.EPS ) THEN
191:          INFO = -13
192:       ELSE IF( NSWEEP.LT.0 ) THEN
193:          INFO = -14
194:       ELSE IF( LWORK.LT.M ) THEN
195:          INFO = -16
196:       ELSE
197:          INFO = 0
198:       END IF
199: *
200: *     #:(
201:       IF( INFO.NE.0 ) THEN
202:          CALL XERBLA( 'DGSVJ0', -INFO )
203:          RETURN
204:       END IF
205: *
206:       IF( RSVEC ) THEN
207:          MVL = N
208:       ELSE IF( APPLV ) THEN
209:          MVL = MV
210:       END IF
211:       RSVEC = RSVEC .OR. APPLV
212: 
213:       ROOTEPS = DSQRT( EPS )
214:       ROOTSFMIN = DSQRT( SFMIN )
215:       SMALL = SFMIN / EPS
216:       BIG = ONE / SFMIN
217:       ROOTBIG = ONE / ROOTSFMIN
218:       BIGTHETA = ONE / ROOTEPS
219:       ROOTTOL = DSQRT( TOL )
220: *
221: *
222: *     -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#-
223: *
224:       EMPTSW = ( N*( N-1 ) ) / 2
225:       NOTROT = 0
226:       FASTR( 1 ) = ZERO
227: *
228: *     -#- Row-cyclic pivot strategy with de Rijk's pivoting -#-
229: *
230: 
231:       SWBAND = 0
232: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
233: *     if SGESVJ is used as a computational routine in the preconditioned
234: *     Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
235: *     ......
236: 
237:       KBL = MIN0( 8, N )
238: *[TP] KBL is a tuning parameter that defines the tile size in the
239: *     tiling of the p-q loops of pivot pairs. In general, an optimal
240: *     value of KBL depends on the matrix dimensions and on the
241: *     parameters of the computer's memory.
242: *
243:       NBL = N / KBL
244:       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
245: 
246:       BLSKIP = ( KBL**2 ) + 1
247: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
248: 
249:       ROWSKIP = MIN0( 5, KBL )
250: *[TP] ROWSKIP is a tuning parameter.
251: 
252:       LKAHEAD = 1
253: *[TP] LKAHEAD is a tuning parameter.
254:       SWBAND = 0
255:       PSKIPPED = 0
256: *
257:       DO 1993 i = 1, NSWEEP
258: *     .. go go go ...
259: *
260:          MXAAPQ = ZERO
261:          MXSINJ = ZERO
262:          ISWROT = 0
263: *
264:          NOTROT = 0
265:          PSKIPPED = 0
266: *
267:          DO 2000 ibr = 1, NBL
268: 
269:             igl = ( ibr-1 )*KBL + 1
270: *
271:             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
272: *
273:                igl = igl + ir1*KBL
274: *
275:                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
276: 
277: *     .. de Rijk's pivoting
278:                   q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
279:                   IF( p.NE.q ) THEN
280:                      CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
281:                      IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
282:      +                                      V( 1, q ), 1 )
283:                      TEMP1 = SVA( p )
284:                      SVA( p ) = SVA( q )
285:                      SVA( q ) = TEMP1
286:                      TEMP1 = D( p )
287:                      D( p ) = D( q )
288:                      D( q ) = TEMP1
289:                   END IF
290: *
291:                   IF( ir1.EQ.0 ) THEN
292: *
293: *        Column norms are periodically updated by explicit
294: *        norm computation.
295: *        Caveat:
296: *        Some BLAS implementations compute DNRM2(M,A(1,p),1)
297: *        as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may result in
298: *        overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and
299: *        undeflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
300: *        Hence, DNRM2 cannot be trusted, not even in the case when
301: *        the true norm is far from the under(over)flow boundaries.
302: *        If properly implemented DNRM2 is available, the IF-THEN-ELSE
303: *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * D(p)".
304: *
305:                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
306:      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
307:                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*D( p )
308:                      ELSE
309:                         TEMP1 = ZERO
310:                         AAPP = ZERO
311:                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
312:                         SVA( p ) = TEMP1*DSQRT( AAPP )*D( p )
313:                      END IF
314:                      AAPP = SVA( p )
315:                   ELSE
316:                      AAPP = SVA( p )
317:                   END IF
318: 
319: *
320:                   IF( AAPP.GT.ZERO ) THEN
321: *
322:                      PSKIPPED = 0
323: *
324:                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
325: *
326:                         AAQQ = SVA( q )
327: 
328:                         IF( AAQQ.GT.ZERO ) THEN
329: *
330:                            AAPP0 = AAPP
331:                            IF( AAQQ.GE.ONE ) THEN
332:                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
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:                               ROTOK = AAPP.LE.( AAQQ / SMALL )
346:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
347:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
348:      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
349:      +                                  / AAPP
350:                               ELSE
351:                                  CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
352:                                  CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
353:      +                                        M, 1, WORK, LDA, IERR )
354:                                  AAPQ = DDOT( M, WORK, 1, A( 1, p ),
355:      +                                  1 )*D( p ) / AAPP
356:                               END IF
357:                            END IF
358: *
359:                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
360: *
361: *        TO rotate or NOT to rotate, THAT is the question ...
362: *
363:                            IF( DABS( AAPQ ).GT.TOL ) THEN
364: *
365: *           .. rotate
366: *           ROTATED = ROTATED + ONE
367: *
368:                               IF( ir1.EQ.0 ) THEN
369:                                  NOTROT = 0
370:                                  PSKIPPED = 0
371:                                  ISWROT = ISWROT + 1
372:                               END IF
373: *
374:                               IF( ROTOK ) THEN
375: *
376:                                  AQOAP = AAQQ / AAPP
377:                                  APOAQ = AAPP / AAQQ
378:                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /
379:      +                                   AAPQ
380: *
381:                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
382: *
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( ONE-T*AQOAP*
395:      +                                     AAPQ )
396:                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
397: *
398:                                  ELSE
399: *
400: *                 .. choose correct signum for THETA and rotate
401: *
402:                                     THSIGN = -DSIGN( ONE, AAPQ )
403:                                     T = ONE / ( THETA+THSIGN*
404:      +                                  DSQRT( ONE+THETA*THETA ) )
405:                                     CS = DSQRT( ONE / ( ONE+T*T ) )
406:                                     SN = T*CS
407: *
408:                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
409:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
410:      +                                         ONE+T*APOAQ*AAPQ ) )
411:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
412:      +                                     ONE-T*AQOAP*AAPQ ) )
413: *
414:                                     APOAQ = D( p ) / D( q )
415:                                     AQOAP = D( q ) / D( p )
416:                                     IF( D( p ).GE.ONE ) THEN
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:                                           D( p ) = D( p )*CS
436:                                           D( q ) = D( q ) / CS
437:                                           IF( RSVEC ) THEN
438:                                              CALL DAXPY( MVL, -T*AQOAP,
439:      +                                                   V( 1, q ), 1,
440:      +                                                   V( 1, p ), 1 )
441:                                              CALL DAXPY( MVL,
442:      +                                                   CS*SN*APOAQ,
443:      +                                                   V( 1, p ), 1,
444:      +                                                   V( 1, q ), 1 )
445:                                           END IF
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:                                           D( p ) = D( p ) / CS
456:                                           D( q ) = D( q )*CS
457:                                           IF( RSVEC ) THEN
458:                                              CALL DAXPY( MVL, T*APOAQ,
459:      +                                                   V( 1, p ), 1,
460:      +                                                   V( 1, q ), 1 )
461:                                              CALL DAXPY( MVL,
462:      +                                                   -CS*SN*AQOAP,
463:      +                                                   V( 1, q ), 1,
464:      +                                                   V( 1, p ), 1 )
465:                                           END IF
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: *              .. have to use modified Gram-Schmidt like transformation
512:                                  CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
513:                                  CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
514:      +                                        1, WORK, LDA, IERR )
515:                                  CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
516:      +                                        1, A( 1, q ), LDA, IERR )
517:                                  TEMP1 = -AAPQ*D( p ) / D( q )
518:                                  CALL DAXPY( M, TEMP1, WORK, 1,
519:      +                                       A( 1, q ), 1 )
520:                                  CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
521:      +                                        1, A( 1, q ), LDA, IERR )
522:                                  SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
523:      +                                      ONE-AAPQ*AAPQ ) )
524:                                  MXSINJ = DMAX1( MXSINJ, SFMIN )
525:                               END IF
526: *           END IF ROTOK THEN ... ELSE
527: *
528: *           In the case of cancellation in updating SVA(q), SVA(p)
529: *           recompute SVA(q), SVA(p).
530:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
531:      +                            THEN
532:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
533:      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
534:                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
535:      +                                         D( q )
536:                                  ELSE
537:                                     T = ZERO
538:                                     AAQQ = ZERO
539:                                     CALL DLASSQ( M, A( 1, q ), 1, T,
540:      +                                           AAQQ )
541:                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
542:                                  END IF
543:                               END IF
544:                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
545:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
546:      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
547:                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
548:      +                                     D( p )
549:                                  ELSE
550:                                     T = ZERO
551:                                     AAPP = ZERO
552:                                     CALL DLASSQ( M, A( 1, p ), 1, T,
553:      +                                           AAPP )
554:                                     AAPP = T*DSQRT( AAPP )*D( p )
555:                                  END IF
556:                                  SVA( p ) = AAPP
557:                               END IF
558: *
559:                            ELSE
560: *        A(:,p) and A(:,q) already numerically orthogonal
561:                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
562:                               PSKIPPED = PSKIPPED + 1
563:                            END IF
564:                         ELSE
565: *        A(:,q) is zero column
566:                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
567:                            PSKIPPED = PSKIPPED + 1
568:                         END IF
569: *
570:                         IF( ( i.LE.SWBAND ) .AND.
571:      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
572:                            IF( ir1.EQ.0 )AAPP = -AAPP
573:                            NOTROT = 0
574:                            GO TO 2103
575:                         END IF
576: *
577:  2002                CONTINUE
578: *     END q-LOOP
579: *
580:  2103                CONTINUE
581: *     bailed out of q-loop
582: 
583:                      SVA( p ) = AAPP
584: 
585:                   ELSE
586:                      SVA( p ) = AAPP
587:                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
588:      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
589:                   END IF
590: *
591:  2001          CONTINUE
592: *     end of the p-loop
593: *     end of doing the block ( ibr, ibr )
594:  1002       CONTINUE
595: *     end of ir1-loop
596: *
597: *........................................................
598: * ... go to the off diagonal blocks
599: *
600:             igl = ( ibr-1 )*KBL + 1
601: *
602:             DO 2010 jbc = ibr + 1, NBL
603: *
604:                jgl = ( jbc-1 )*KBL + 1
605: *
606: *        doing the block at ( ibr, jbc )
607: *
608:                IJBLSK = 0
609:                DO 2100 p = igl, MIN0( igl+KBL-1, N )
610: *
611:                   AAPP = SVA( p )
612: *
613:                   IF( AAPP.GT.ZERO ) THEN
614: *
615:                      PSKIPPED = 0
616: *
617:                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
618: *
619:                         AAQQ = SVA( q )
620: *
621:                         IF( AAQQ.GT.ZERO ) THEN
622:                            AAPP0 = AAPP
623: *
624: *     -#- M x 2 Jacobi SVD -#-
625: *
626: *        -#- Safe Gram matrix computation -#-
627: *
628:                            IF( AAQQ.GE.ONE ) THEN
629:                               IF( AAPP.GE.AAQQ ) THEN
630:                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
631:                               ELSE
632:                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
633:                               END IF
634:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
635:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
636:      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
637:      +                                  / AAPP
638:                               ELSE
639:                                  CALL DCOPY( M, A( 1, p ), 1, WORK, 1 )
640:                                  CALL DLASCL( 'G', 0, 0, AAPP, D( p ),
641:      +                                        M, 1, WORK, LDA, IERR )
642:                                  AAPQ = DDOT( M, WORK, 1, A( 1, q ),
643:      +                                  1 )*D( q ) / AAQQ
644:                               END IF
645:                            ELSE
646:                               IF( AAPP.GE.AAQQ ) THEN
647:                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
648:                               ELSE
649:                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
650:                               END IF
651:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
652:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
653:      +                                  q ), 1 )*D( p )*D( q ) / AAQQ )
654:      +                                  / AAPP
655:                               ELSE
656:                                  CALL DCOPY( M, A( 1, q ), 1, WORK, 1 )
657:                                  CALL DLASCL( 'G', 0, 0, AAQQ, D( q ),
658:      +                                        M, 1, WORK, LDA, IERR )
659:                                  AAPQ = DDOT( M, WORK, 1, A( 1, p ),
660:      +                                  1 )*D( p ) / AAPP
661:                               END IF
662:                            END IF
663: *
664:                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
665: *
666: *        TO rotate or NOT to rotate, THAT is the question ...
667: *
668:                            IF( DABS( AAPQ ).GT.TOL ) THEN
669:                               NOTROT = 0
670: *           ROTATED  = ROTATED + 1
671:                               PSKIPPED = 0
672:                               ISWROT = ISWROT + 1
673: *
674:                               IF( ROTOK ) THEN
675: *
676:                                  AQOAP = AAQQ / AAPP
677:                                  APOAQ = AAPP / AAQQ
678:                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /
679:      +                                   AAPQ
680:                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
681: *
682:                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
683:                                     T = HALF / THETA
684:                                     FASTR( 3 ) = T*D( p ) / D( q )
685:                                     FASTR( 4 ) = -T*D( q ) / D( p )
686:                                     CALL DROTM( M, A( 1, p ), 1,
687:      +                                          A( 1, q ), 1, FASTR )
688:                                     IF( RSVEC )CALL DROTM( MVL,
689:      +                                              V( 1, p ), 1,
690:      +                                              V( 1, q ), 1,
691:      +                                              FASTR )
692:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
693:      +                                         ONE+T*APOAQ*AAPQ ) )
694:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
695:      +                                     ONE-T*AQOAP*AAPQ ) )
696:                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
697:                                  ELSE
698: *
699: *                 .. choose correct signum for THETA and rotate
700: *
701:                                     THSIGN = -DSIGN( ONE, AAPQ )
702:                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
703:                                     T = ONE / ( THETA+THSIGN*
704:      +                                  DSQRT( ONE+THETA*THETA ) )
705:                                     CS = DSQRT( ONE / ( ONE+T*T ) )
706:                                     SN = T*CS
707:                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
708:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
709:      +                                         ONE+T*APOAQ*AAPQ ) )
710:                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*
711:      +                                     AAPQ )
712: *
713:                                     APOAQ = D( p ) / D( q )
714:                                     AQOAP = D( q ) / D( p )
715:                                     IF( D( p ).GE.ONE ) THEN
716: *
717:                                        IF( D( q ).GE.ONE ) THEN
718:                                           FASTR( 3 ) = T*APOAQ
719:                                           FASTR( 4 ) = -T*AQOAP
720:                                           D( p ) = D( p )*CS
721:                                           D( q ) = D( q )*CS
722:                                           CALL DROTM( M, A( 1, p ), 1,
723:      +                                                A( 1, q ), 1,
724:      +                                                FASTR )
725:                                           IF( RSVEC )CALL DROTM( MVL,
726:      +                                        V( 1, p ), 1, V( 1, q ),
727:      +                                        1, FASTR )
728:                                        ELSE
729:                                           CALL DAXPY( M, -T*AQOAP,
730:      +                                                A( 1, q ), 1,
731:      +                                                A( 1, p ), 1 )
732:                                           CALL DAXPY( M, CS*SN*APOAQ,
733:      +                                                A( 1, p ), 1,
734:      +                                                A( 1, q ), 1 )
735:                                           IF( RSVEC ) THEN
736:                                              CALL DAXPY( MVL, -T*AQOAP,
737:      +                                                   V( 1, q ), 1,
738:      +                                                   V( 1, p ), 1 )
739:                                              CALL DAXPY( MVL,
740:      +                                                   CS*SN*APOAQ,
741:      +                                                   V( 1, p ), 1,
742:      +                                                   V( 1, q ), 1 )
743:                                           END IF
744:                                           D( p ) = D( p )*CS
745:                                           D( q ) = D( q ) / CS
746:                                        END IF
747:                                     ELSE
748:                                        IF( D( q ).GE.ONE ) THEN
749:                                           CALL DAXPY( M, T*APOAQ,
750:      +                                                A( 1, p ), 1,
751:      +                                                A( 1, q ), 1 )
752:                                           CALL DAXPY( M, -CS*SN*AQOAP,
753:      +                                                A( 1, q ), 1,
754:      +                                                A( 1, p ), 1 )
755:                                           IF( RSVEC ) THEN
756:                                              CALL DAXPY( MVL, T*APOAQ,
757:      +                                                   V( 1, p ), 1,
758:      +                                                   V( 1, q ), 1 )
759:                                              CALL DAXPY( MVL,
760:      +                                                   -CS*SN*AQOAP,
761:      +                                                   V( 1, q ), 1,
762:      +                                                   V( 1, p ), 1 )
763:                                           END IF
764:                                           D( p ) = D( p ) / CS
765:                                           D( q ) = D( q )*CS
766:                                        ELSE
767:                                           IF( D( p ).GE.D( q ) ) THEN
768:                                              CALL DAXPY( M, -T*AQOAP,
769:      +                                                   A( 1, q ), 1,
770:      +                                                   A( 1, p ), 1 )
771:                                              CALL DAXPY( M, CS*SN*APOAQ,
772:      +                                                   A( 1, p ), 1,
773:      +                                                   A( 1, q ), 1 )
774:                                              D( p ) = D( p )*CS
775:                                              D( q ) = D( q ) / CS
776:                                              IF( RSVEC ) THEN
777:                                                 CALL DAXPY( MVL,
778:      +                                               -T*AQOAP,
779:      +                                               V( 1, q ), 1,
780:      +                                               V( 1, p ), 1 )
781:                                                 CALL DAXPY( MVL,
782:      +                                               CS*SN*APOAQ,
783:      +                                               V( 1, p ), 1,
784:      +                                               V( 1, q ), 1 )
785:                                              END IF
786:                                           ELSE
787:                                              CALL DAXPY( M, T*APOAQ,
788:      +                                                   A( 1, p ), 1,
789:      +                                                   A( 1, q ), 1 )
790:                                              CALL DAXPY( M,
791:      +                                                   -CS*SN*AQOAP,
792:      +                                                   A( 1, q ), 1,
793:      +                                                   A( 1, p ), 1 )
794:                                              D( p ) = D( p ) / CS
795:                                              D( q ) = D( q )*CS
796:                                              IF( RSVEC ) THEN
797:                                                 CALL DAXPY( MVL,
798:      +                                               T*APOAQ, V( 1, p ),
799:      +                                               1, V( 1, q ), 1 )
800:                                                 CALL DAXPY( MVL,
801:      +                                               -CS*SN*AQOAP,
802:      +                                               V( 1, q ), 1,
803:      +                                               V( 1, p ), 1 )
804:                                              END IF
805:                                           END IF
806:                                        END IF
807:                                     END IF
808:                                  END IF
809: *
810:                               ELSE
811:                                  IF( AAPP.GT.AAQQ ) THEN
812:                                     CALL DCOPY( M, A( 1, p ), 1, WORK,
813:      +                                          1 )
814:                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
815:      +                                           M, 1, WORK, LDA, IERR )
816:                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
817:      +                                           M, 1, A( 1, q ), LDA,
818:      +                                           IERR )
819:                                     TEMP1 = -AAPQ*D( p ) / D( q )
820:                                     CALL DAXPY( M, TEMP1, WORK, 1,
821:      +                                          A( 1, q ), 1 )
822:                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
823:      +                                           M, 1, A( 1, q ), LDA,
824:      +                                           IERR )
825:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
826:      +                                         ONE-AAPQ*AAPQ ) )
827:                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
828:                                  ELSE
829:                                     CALL DCOPY( M, A( 1, q ), 1, WORK,
830:      +                                          1 )
831:                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
832:      +                                           M, 1, WORK, LDA, IERR )
833:                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
834:      +                                           M, 1, A( 1, p ), LDA,
835:      +                                           IERR )
836:                                     TEMP1 = -AAPQ*D( q ) / D( p )
837:                                     CALL DAXPY( M, TEMP1, WORK, 1,
838:      +                                          A( 1, p ), 1 )
839:                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,
840:      +                                           M, 1, A( 1, p ), LDA,
841:      +                                           IERR )
842:                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
843:      +                                         ONE-AAPQ*AAPQ ) )
844:                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
845:                                  END IF
846:                               END IF
847: *           END IF ROTOK THEN ... ELSE
848: *
849: *           In the case of cancellation in updating SVA(q)
850: *           .. recompute SVA(q)
851:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
852:      +                            THEN
853:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
854:      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
855:                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
856:      +                                         D( q )
857:                                  ELSE
858:                                     T = ZERO
859:                                     AAQQ = ZERO
860:                                     CALL DLASSQ( M, A( 1, q ), 1, T,
861:      +                                           AAQQ )
862:                                     SVA( q ) = T*DSQRT( AAQQ )*D( q )
863:                                  END IF
864:                               END IF
865:                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
866:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
867:      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN
868:                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
869:      +                                     D( p )
870:                                  ELSE
871:                                     T = ZERO
872:                                     AAPP = ZERO
873:                                     CALL DLASSQ( M, A( 1, p ), 1, T,
874:      +                                           AAPP )
875:                                     AAPP = T*DSQRT( AAPP )*D( p )
876:                                  END IF
877:                                  SVA( p ) = AAPP
878:                               END IF
879: *              end of OK rotation
880:                            ELSE
881:                               NOTROT = NOTROT + 1
882:                               PSKIPPED = PSKIPPED + 1
883:                               IJBLSK = IJBLSK + 1
884:                            END IF
885:                         ELSE
886:                            NOTROT = NOTROT + 1
887:                            PSKIPPED = PSKIPPED + 1
888:                            IJBLSK = IJBLSK + 1
889:                         END IF
890: *
891:                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
892:      +                      THEN
893:                            SVA( p ) = AAPP
894:                            NOTROT = 0
895:                            GO TO 2011
896:                         END IF
897:                         IF( ( i.LE.SWBAND ) .AND.
898:      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
899:                            AAPP = -AAPP
900:                            NOTROT = 0
901:                            GO TO 2203
902:                         END IF
903: *
904:  2200                CONTINUE
905: *        end of the q-loop
906:  2203                CONTINUE
907: *
908:                      SVA( p ) = AAPP
909: *
910:                   ELSE
911:                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
912:      +                   MIN0( jgl+KBL-1, N ) - jgl + 1
913:                      IF( AAPP.LT.ZERO )NOTROT = 0
914:                   END IF
915: 
916:  2100          CONTINUE
917: *     end of the p-loop
918:  2010       CONTINUE
919: *     end of the jbc-loop
920:  2011       CONTINUE
921: *2011 bailed out of the jbc-loop
922:             DO 2012 p = igl, MIN0( igl+KBL-1, N )
923:                SVA( p ) = DABS( SVA( p ) )
924:  2012       CONTINUE
925: *
926:  2000    CONTINUE
927: *2000 :: end of the ibr-loop
928: *
929: *     .. update SVA(N)
930:          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
931:      +       THEN
932:             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N )
933:          ELSE
934:             T = ZERO
935:             AAPP = ZERO
936:             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
937:             SVA( N ) = T*DSQRT( AAPP )*D( N )
938:          END IF
939: *
940: *     Additional steering devices
941: *
942:          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
943:      +       ( ISWROT.LE.N ) ) )SWBAND = i
944: *
945:          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND.
946:      +       ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
947:             GO TO 1994
948:          END IF
949: *
950:          IF( NOTROT.GE.EMPTSW )GO TO 1994
951: 
952:  1993 CONTINUE
953: *     end i=1:NSWEEP loop
954: * #:) Reaching this point means that the procedure has comleted the given
955: *     number of iterations.
956:       INFO = NSWEEP - 1
957:       GO TO 1995
958:  1994 CONTINUE
959: * #:) Reaching this point means that during the i-th sweep all pivots were
960: *     below the given tolerance, causing early exit.
961: *
962:       INFO = 0
963: * #:) INFO = 0 confirms successful iterations.
964:  1995 CONTINUE
965: *
966: *     Sort the vector D.
967:       DO 5991 p = 1, N - 1
968:          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
969:          IF( p.NE.q ) THEN
970:             TEMP1 = SVA( p )
971:             SVA( p ) = SVA( q )
972:             SVA( q ) = TEMP1
973:             TEMP1 = D( p )
974:             D( p ) = D( q )
975:             D( q ) = TEMP1
976:             CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
977:             IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
978:          END IF
979:  5991 CONTINUE
980: *
981:       RETURN
982: *     ..
983: *     .. END OF DGSVJ0
984: *     ..
985:       END
986: