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