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)                                    --
005: *
006: *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --
007: *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --
008: *  -- November 2008                                                   --
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: *     Scalar Arguments
019: *
020:       IMPLICIT    NONE
021:       INTEGER     INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
022:       DOUBLE PRECISION  EPS, SFMIN, TOL
023:       CHARACTER*1 JOBV
024: *
025: *     Array Arguments
026: *
027:       DOUBLE PRECISION A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
028:      &                 WORK( LWORK )
029: *     ..
030: *
031: *  Purpose
032: *  ~~~~~~~
033: *  DGSVJ0 is called from DGESVJ as a pre-processor and that is its main
034: *  purpose. It applies Jacobi rotations in the same way as DGESVJ 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: *  DGSVJ0 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 DABS(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: *     Local Parameters
144:       DOUBLE PRECISION   ZERO,  HALF,         ONE,         TWO
145:       PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0 )
146: 
147: *     Local Scalars
148:       DOUBLE PRECISION AAPP,  AAPP0, AAPQ,   AAQQ,   APOAQ,   AQOAP,
149:      &          BIG,   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
150:      &          ROOTSFMIN, ROOTTOL,  SMALL,  SN, T,  TEMP1,   THETA,
151:      &          THSIGN
152:       INTEGER   BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, ISWROT,
153:      &          jbc, jgl, KBL, LKAHEAD, MVL, NBL, NOTROT, p, PSKIPPED,
154:      &          q, ROWSKIP, SWBAND
155:       LOGICAL   APPLV, ROTOK, RSVEC
156: 
157: *     Local Arrays
158: *
159:       DOUBLE PRECISION  FASTR(5)
160: *
161: *     Intrinsic Functions
162: *
163:       INTRINSIC DABS, DMAX1, DBLE, MIN0, DSIGN, DSQRT
164: *
165: *     External Functions
166: *
167:       DOUBLE PRECISION  DDOT, DNRM2
168:       INTEGER   IDAMAX
169:       LOGICAL   LSAME
170:       EXTERNAL  IDAMAX, LSAME, DDOT, DNRM2
171: *
172: *     External Subroutines
173: *
174:       EXTERNAL  DAXPY, DCOPY, DLASCL, DLASSQ, DROTM, DSWAP
175: *
176: *     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~|
177: *
178:       APPLV = LSAME(JOBV,'A')
179:       RSVEC = LSAME(JOBV,'V')
180:       IF ( .NOT.( RSVEC .OR. APPLV .OR. LSAME(JOBV,'N'))) THEN
181:          INFO = -1
182:       ELSE IF ( M .LT. 0 ) THEN
183:          INFO = -2
184:       ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M )) THEN
185:          INFO = -3
186:       ELSE IF ( LDA .LT. M ) THEN
187:          INFO = -5
188:       ELSE IF ( MV .LT. 0 ) THEN
189:          INFO = -8
190:       ELSE IF ( LDV .LT. M ) THEN
191:          INFO = -10
192:       ELSE IF ( TOL .LE. EPS ) THEN
193:          INFO = -13
194:       ELSE IF ( NSWEEP .LT. 0 ) THEN
195:          INFO = -14
196:       ELSE IF ( LWORK .LT. M ) THEN
197:          INFO = -16
198:       ELSE
199:          INFO = 0
200:       END IF
201: *
202: *     #:(
203:       IF ( INFO .NE. 0 ) THEN
204:          CALL XERBLA( 'DGSVJ0', -INFO )
205:          RETURN
206:       END IF
207: *
208:       IF ( RSVEC ) THEN
209:          MVL = N
210:       ELSE IF ( APPLV ) THEN
211:          MVL = MV
212:       END IF
213:       RSVEC = RSVEC .OR. APPLV
214: 
215:       ROOTEPS     = DSQRT(EPS)
216:       ROOTSFMIN   = DSQRT(SFMIN)
217:       SMALL       = SFMIN  / EPS
218:       BIG         = ONE   / SFMIN
219:       ROOTBIG     = ONE  / ROOTSFMIN
220:       BIGTHETA    = ONE  / ROOTEPS
221:       ROOTTOL     = DSQRT(TOL)
222: *
223: *
224: *     -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#-
225: *
226:       EMPTSW   = ( N * ( N - 1 ) ) / 2
227:       NOTROT   = 0
228:       FASTR(1) = ZERO
229: *
230: *     -#- Row-cyclic pivot strategy with de Rijk's pivoting -#-
231: *
232: 
233:       SWBAND = 0
234: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
235: *     if SGESVJ is used as a computational routine in the preconditioned
236: *     Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
237: *     ......
238: 
239:       KBL = MIN0( 8, N )
240: *[TP] KBL is a tuning parameter that defines the tile size in the
241: *     tiling of the p-q loops of pivot pairs. In general, an optimal
242: *     value of KBL depends on the matrix dimensions and on the
243: *     parameters of the computer's memory.
244: *
245:       NBL = N / KBL
246:       IF ( ( NBL * KBL ) .NE. N ) NBL = NBL + 1
247: 
248:       BLSKIP = ( KBL**2 ) + 1
249: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
250: 
251:       ROWSKIP = MIN0( 5, KBL )
252: *[TP] ROWSKIP is a tuning parameter.
253: 
254:       LKAHEAD = 1
255: *[TP] LKAHEAD is a tuning parameter.
256:       SWBAND = 0
257:       PSKIPPED = 0
258: *
259:       DO 1993 i = 1, NSWEEP
260: *     .. go go go ...
261: *
262:       MXAAPQ = ZERO
263:       MXSINJ = ZERO
264:       ISWROT = 0
265: *
266:       NOTROT = 0
267:       PSKIPPED = 0
268: *
269:       DO 2000 ibr = 1, NBL
270: 
271:       igl = ( ibr - 1 ) * KBL + 1
272: *
273:       DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL - ibr )
274: *
275:       igl = igl + ir1 * KBL
276: *
277:       DO 2001 p = igl, MIN0( igl + KBL - 1, N - 1)
278: 
279: *     .. de Rijk's pivoting
280:       q   = IDAMAX( N-p+1, SVA(p), 1 ) + p - 1
281:       IF ( p .NE. q ) THEN
282:          CALL DSWAP( M, A(1,p), 1, A(1,q), 1 )
283:          IF ( RSVEC ) CALL DSWAP( MVL, V(1,p), 1, V(1,q), 1 )
284:          TEMP1   = SVA(p)
285:          SVA(p)  = SVA(q)
286:          SVA(q)  = TEMP1
287:          TEMP1   = D(p)
288:          D(p) = D(q)
289:          D(q) = TEMP1
290:       END IF
291: *
292:       IF ( ir1 .EQ. 0 ) THEN
293: *
294: *        Column norms are periodically updated by explicit
295: *        norm computation.
296: *        Caveat:
297: *        Some BLAS implementations compute DNRM2(M,A(1,p),1)
298: *        as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may result in
299: *        overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and
300: *        undeflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
301: *        Hence, DNRM2 cannot be trusted, not even in the case when
302: *        the true norm is far from the under(over)flow boundaries.
303: *        If properly implemented DNRM2 is available, the IF-THEN-ELSE
304: *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * D(p)".
305: *
306:          IF ((SVA(p) .LT. ROOTBIG) .AND. (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,q), 1 ) *
335:      &                  D(p) * D(q) / AAQQ ) / AAPP
336:             ELSE
337:                CALL DCOPY( M, A(1,p), 1, WORK, 1 )
338:                CALL DLASCL( 'G', 0, 0, AAPP, D(p), M,
339:      &              1, WORK, LDA, IERR )
340:                AAPQ = DDOT( M, WORK,1, A(1,q),1 )*D(q) / AAQQ
341:             END IF
342:          ELSE
343:             ROTOK  = AAPP .LE. ( AAQQ / SMALL )
344:             IF ( AAPP .GT.  ( SMALL / AAQQ ) ) THEN
345:                AAPQ = ( DDOT( M, A(1,p), 1, A(1,q), 1 ) *
346:      &               D(p) * D(q) / AAQQ ) / AAPP
347:             ELSE
348:                CALL DCOPY( M, A(1,q), 1, WORK, 1 )
349:                CALL DLASCL( 'G', 0, 0, AAQQ, D(q), M,
350:      &              1, WORK, LDA, IERR )
351:                AAPQ = DDOT( M, WORK,1, A(1,p),1 )*D(p) / AAPP
352:             END IF
353:          END IF
354: *
355:          MXAAPQ = DMAX1( MXAAPQ, DABS(AAPQ) )
356: *
357: *        TO rotate or NOT to rotate, THAT is the question ...
358: *
359:          IF ( DABS( AAPQ ) .GT. TOL ) THEN
360: *
361: *           .. rotate
362: *           ROTATED = ROTATED + ONE
363: *
364:             IF ( ir1 .EQ. 0 ) THEN
365:                NOTROT   = 0
366:                PSKIPPED = 0
367:                ISWROT   = ISWROT  + 1
368:             END IF
369: *
370:             IF ( ROTOK ) THEN
371: *
372:                AQOAP = AAQQ / AAPP
373:                APOAQ = AAPP / AAQQ
374:                THETA = - HALF * DABS( AQOAP - APOAQ ) / AAPQ
375: *
376:                IF ( DABS( THETA ) .GT. BIGTHETA ) THEN
377: *
378:                   T        = HALF / THETA
379:                   FASTR(3) =   T * D(p) / D(q)
380:                   FASTR(4) = - T * D(q) / D(p)
381:                   CALL DROTM( M,   A(1,p), 1, A(1,q), 1, FASTR )
382:                   IF ( RSVEC )
383:      &            CALL DROTM( MVL, V(1,p), 1, V(1,q), 1, FASTR )
384:                   SVA(q) = AAQQ*DSQRT( DMAX1(ZERO,ONE + T*APOAQ*AAPQ) )
385:                   AAPP   = AAPP*DSQRT( ONE - T*AQOAP*AAPQ )
386:                   MXSINJ = DMAX1( MXSINJ, DABS(T) )
387: *
388:                ELSE
389: *
390: *                 .. choose correct signum for THETA and rotate
391: *
392:                   THSIGN =  - DSIGN(ONE,AAPQ)
393:                   T  = ONE / ( THETA + THSIGN*DSQRT(ONE+THETA*THETA) )
394:                   CS = DSQRT( ONE / ( ONE + T*T ) )
395:                   SN = T * CS
396: *
397:                   MXSINJ = DMAX1( MXSINJ, DABS(SN) )
398:                   SVA(q) = AAQQ*DSQRT( DMAX1(ZERO, ONE+T*APOAQ*AAPQ) )
399:                   AAPP   = AAPP*DSQRT( DMAX1(ZERO, ONE-T*AQOAP*AAPQ) )
400: *
401:                   APOAQ = D(p) / D(q)
402:                   AQOAP = D(q) / D(p)
403:                   IF ( D(p) .GE. ONE ) THEN
404:                      IF ( D(q) .GE.  ONE ) THEN
405:                         FASTR(3) =   T * APOAQ
406:                         FASTR(4) = - T * AQOAP
407:                         D(p)  = D(p) * CS
408:                         D(q)  = D(q) * CS
409:                         CALL DROTM( M,   A(1,p),1, A(1,q),1, FASTR )
410:                         IF ( RSVEC )
411:      &                  CALL DROTM( MVL, V(1,p),1, V(1,q),1, FASTR )
412:                      ELSE
413:                         CALL DAXPY( M,    -T*AQOAP, A(1,q),1, A(1,p),1 )
414:                         CALL DAXPY( M, CS*SN*APOAQ, A(1,p),1, A(1,q),1 )
415:                         D(p) = D(p) * CS
416:                         D(q) = D(q) / CS
417:                         IF ( RSVEC ) THEN
418:                         CALL DAXPY(MVL,   -T*AQOAP, V(1,q),1,V(1,p),1)
419:                         CALL DAXPY(MVL,CS*SN*APOAQ, V(1,p),1,V(1,q),1)
420:                         END IF
421:                      END IF
422:                   ELSE
423:                      IF ( D(q) .GE. ONE ) THEN
424:                         CALL DAXPY( M,     T*APOAQ, A(1,p),1, A(1,q),1 )
425:                         CALL DAXPY( M,-CS*SN*AQOAP, A(1,q),1, A(1,p),1 )
426:                         D(p) = D(p) / CS
427:                         D(q) = D(q) * CS
428:                         IF ( RSVEC ) THEN
429:                         CALL DAXPY(MVL, T*APOAQ,    V(1,p),1,V(1,q),1)
430:                         CALL DAXPY(MVL,-CS*SN*AQOAP,V(1,q),1,V(1,p),1)
431:                         END IF
432:                      ELSE
433:                         IF ( D(p) .GE. D(q) ) THEN
434:                            CALL DAXPY( M,-T*AQOAP,   A(1,q),1,A(1,p),1 )
435:                            CALL DAXPY( M,CS*SN*APOAQ,A(1,p),1,A(1,q),1 )
436:                            D(p) = D(p) * CS
437:                            D(q) = D(q) / CS
438:                            IF ( RSVEC ) THEN
439:                            CALL DAXPY(MVL, -T*AQOAP,  V(1,q),1,V(1,p),1)
440:                            CALL DAXPY(MVL,CS*SN*APOAQ,V(1,p),1,V(1,q),1)
441:                            END IF
442:                         ELSE
443:                            CALL DAXPY( M, T*APOAQ,    A(1,p),1,A(1,q),1)
444:                            CALL DAXPY( M,-CS*SN*AQOAP,A(1,q),1,A(1,p),1)
445:                            D(p) = D(p) / CS
446:                            D(q) = D(q) * CS
447:                           IF ( RSVEC ) THEN
448:                           CALL DAXPY(MVL, T*APOAQ,    V(1,p),1,V(1,q),1)
449:                           CALL DAXPY(MVL,-CS*SN*AQOAP,V(1,q),1,V(1,p),1)
450:                           END IF
451:                         END IF
452:                      END IF
453:                   ENDIF
454:                END IF
455: *
456:             ELSE
457: *              .. have to use modified Gram-Schmidt like transformation
458:                CALL DCOPY( M, A(1,p), 1, WORK, 1 )
459:                CALL DLASCL( 'G',0,0,AAPP,ONE,M,1,WORK,LDA,IERR )
460:                CALL DLASCL( 'G',0,0,AAQQ,ONE,M,1,   A(1,q),LDA,IERR )
461:                TEMP1 = -AAPQ * D(p) / D(q)
462:                CALL DAXPY ( M, TEMP1, WORK, 1, A(1,q), 1 )
463:                CALL DLASCL( 'G',0,0,ONE,AAQQ,M,1,   A(1,q),LDA,IERR )
464:                SVA(q) = AAQQ*DSQRT( DMAX1( ZERO, ONE - AAPQ*AAPQ ) )
465:                MXSINJ = DMAX1( MXSINJ, SFMIN )
466:             END IF
467: *           END IF ROTOK THEN ... ELSE
468: *
469: *           In the case of cancellation in updating SVA(q), SVA(p)
470: *           recompute SVA(q), SVA(p).
471:             IF ( (SVA(q) / AAQQ )**2 .LE. ROOTEPS  ) THEN
472:                IF ((AAQQ .LT. ROOTBIG).AND.(AAQQ .GT. ROOTSFMIN)) THEN
473:                   SVA(q) = DNRM2( M, A(1,q), 1 ) * D(q)
474:                ELSE
475:                   T    = ZERO
476:                   AAQQ = ZERO
477:                   CALL DLASSQ( M, A(1,q), 1, T, AAQQ )
478:                   SVA(q) = T * DSQRT(AAQQ) * D(q)
479:                END IF
480:             END IF
481:             IF ( ( AAPP / AAPP0) .LE. ROOTEPS  ) THEN
482:                IF ((AAPP .LT. ROOTBIG).AND.(AAPP .GT. ROOTSFMIN)) THEN
483:                   AAPP = DNRM2( M, A(1,p), 1 ) * D(p)
484:                ELSE
485:                   T    = ZERO
486:                   AAPP = ZERO
487:                   CALL DLASSQ( M, A(1,p), 1, T, AAPP )
488:                   AAPP = T * DSQRT(AAPP) * D(p)
489:                END IF
490:                SVA(p) = AAPP
491:             END IF
492: *
493:          ELSE
494: *        A(:,p) and A(:,q) already numerically orthogonal
495:             IF ( ir1 .EQ. 0 ) NOTROT   = NOTROT + 1
496:             PSKIPPED = PSKIPPED + 1
497:          END IF
498:       ELSE
499: *        A(:,q) is zero column
500:          IF ( ir1. EQ. 0 ) NOTROT = NOTROT + 1
501:          PSKIPPED = PSKIPPED + 1
502:       END IF
503: *
504:       IF ( ( i .LE. SWBAND ) .AND. ( PSKIPPED .GT. ROWSKIP ) ) THEN
505:          IF ( ir1 .EQ. 0 ) AAPP = - AAPP
506:          NOTROT = 0
507:          GO TO 2103
508:       END IF
509: *
510:  2002 CONTINUE
511: *     END q-LOOP
512: *
513:  2103 CONTINUE
514: *     bailed out of q-loop
515: 
516:       SVA(p) = AAPP
517: 
518:       ELSE
519:          SVA(p) = AAPP
520:          IF ( ( ir1 .EQ. 0 ) .AND. (AAPP .EQ. ZERO) )
521:      &        NOTROT=NOTROT+MIN0(igl+KBL-1,N)-p
522:       END IF
523: *
524:  2001 CONTINUE
525: *     end of the p-loop
526: *     end of doing the block ( ibr, ibr )
527:  1002 CONTINUE
528: *     end of ir1-loop
529: *
530: *........................................................
531: * ... go to the off diagonal blocks
532: *
533:       igl = ( ibr - 1 ) * KBL + 1
534: *
535:       DO 2010 jbc = ibr + 1, NBL
536: *
537:          jgl = ( jbc - 1 ) * KBL + 1
538: *
539: *        doing the block at ( ibr, jbc )
540: *
541:          IJBLSK = 0
542:          DO 2100 p = igl, MIN0( igl + KBL - 1, N )
543: *
544:          AAPP = SVA(p)
545: *
546:          IF ( AAPP .GT. ZERO ) THEN
547: *
548:          PSKIPPED = 0
549: *
550:          DO 2200 q = jgl, MIN0( jgl + KBL - 1, N )
551: *
552:          AAQQ = SVA(q)
553: *
554:          IF ( AAQQ .GT. ZERO ) THEN
555:             AAPP0 = AAPP
556: *
557: *     -#- M x 2 Jacobi SVD -#-
558: *
559: *        -#- Safe Gram matrix computation -#-
560: *
561:          IF ( AAQQ .GE. ONE ) THEN
562:             IF ( AAPP .GE. AAQQ ) THEN
563:                ROTOK = ( SMALL*AAPP ) .LE. AAQQ
564:             ELSE
565:                ROTOK = ( SMALL*AAQQ ) .LE. AAPP
566:             END IF
567:             IF ( AAPP .LT. ( BIG / AAQQ ) ) THEN
568:                AAPQ = ( DDOT(M, A(1,p), 1, A(1,q), 1 ) *
569:      &                  D(p) * D(q) / AAQQ ) / AAPP
570:             ELSE
571:                CALL DCOPY( M, A(1,p), 1, WORK, 1 )
572:                CALL DLASCL( 'G', 0, 0, AAPP, D(p), M,
573:      &              1, WORK, LDA, IERR )
574:                AAPQ = DDOT( M, WORK, 1, A(1,q), 1 ) *
575:      &                D(q) / AAQQ
576:             END IF
577:          ELSE
578:             IF ( AAPP .GE. AAQQ ) THEN
579:                ROTOK = AAPP .LE. ( AAQQ / SMALL )
580:             ELSE
581:                ROTOK = AAQQ .LE. ( AAPP / SMALL )
582:             END IF
583:             IF ( AAPP .GT.  ( SMALL / AAQQ ) ) THEN
584:                AAPQ = ( DDOT( M, A(1,p), 1, A(1,q), 1 ) *
585:      &               D(p) * D(q) / AAQQ ) / AAPP
586:             ELSE
587:                CALL DCOPY( M, A(1,q), 1, WORK, 1 )
588:                CALL DLASCL( 'G', 0, 0, AAQQ, D(q), M, 1,
589:      &              WORK, LDA, IERR )
590:                AAPQ = DDOT(M,WORK,1,A(1,p),1) * D(p) / AAPP
591:             END IF
592:          END IF
593: *
594:          MXAAPQ = DMAX1( MXAAPQ, DABS(AAPQ) )
595: *
596: *        TO rotate or NOT to rotate, THAT is the question ...
597: *
598:          IF ( DABS( AAPQ ) .GT. TOL ) THEN
599:             NOTROT   = 0
600: *           ROTATED  = ROTATED + 1
601:             PSKIPPED = 0
602:             ISWROT   = ISWROT  + 1
603: *
604:             IF ( ROTOK ) THEN
605: *
606:                AQOAP = AAQQ / AAPP
607:                APOAQ = AAPP / AAQQ
608:                THETA = - HALF * DABS( AQOAP - APOAQ ) / AAPQ
609:                IF ( AAQQ .GT. AAPP0 ) THETA = - THETA
610: *
611:                IF ( DABS( THETA ) .GT. BIGTHETA ) THEN
612:                   T = HALF / THETA
613:                   FASTR(3) =  T * D(p) / D(q)
614:                   FASTR(4) = -T * D(q) / D(p)
615:                   CALL DROTM( M,   A(1,p), 1, A(1,q), 1, FASTR )
616:                   IF ( RSVEC )
617:      &            CALL DROTM( MVL, V(1,p), 1, V(1,q), 1, FASTR )
618:                   SVA(q) = AAQQ*DSQRT( DMAX1(ZERO,ONE + T*APOAQ*AAPQ) )
619:                   AAPP   = AAPP*DSQRT( DMAX1(ZERO,ONE - T*AQOAP*AAPQ) )
620:                   MXSINJ = DMAX1( MXSINJ, DABS(T) )
621:                ELSE
622: *
623: *                 .. choose correct signum for THETA and rotate
624: *
625:                   THSIGN = - DSIGN(ONE,AAPQ)
626:                   IF ( AAQQ .GT. AAPP0 ) THSIGN = - THSIGN
627:                   T  = ONE / ( THETA + THSIGN*DSQRT(ONE+THETA*THETA) )
628:                   CS = DSQRT( ONE / ( ONE + T*T ) )
629:                   SN = T * CS
630:                   MXSINJ = DMAX1( MXSINJ, DABS(SN) )
631:                   SVA(q) = AAQQ*DSQRT( DMAX1(ZERO, ONE+T*APOAQ*AAPQ) )
632:                   AAPP   = AAPP*DSQRT( ONE - T*AQOAP*AAPQ)
633: *
634:                   APOAQ = D(p) / D(q)
635:                   AQOAP = D(q) / D(p)
636:                   IF ( D(p) .GE. ONE ) THEN
637: *
638:                      IF ( D(q) .GE.  ONE ) THEN
639:                         FASTR(3) =   T * APOAQ
640:                         FASTR(4) = - T * AQOAP
641:                         D(p)  = D(p) * CS
642:                         D(q)  = D(q) * CS
643:                         CALL DROTM( M,   A(1,p),1, A(1,q),1, FASTR )
644:                         IF ( RSVEC )
645:      &                  CALL DROTM( MVL, V(1,p),1, V(1,q),1, FASTR )
646:                      ELSE
647:                         CALL DAXPY( M,    -T*AQOAP, A(1,q),1, A(1,p),1 )
648:                         CALL DAXPY( M, CS*SN*APOAQ, A(1,p),1, A(1,q),1 )
649:                         IF ( RSVEC ) THEN
650:                         CALL DAXPY( MVL, -T*AQOAP,  V(1,q),1, V(1,p),1 )
651:                         CALL DAXPY( MVL,CS*SN*APOAQ,V(1,p),1, V(1,q),1 )
652:                         END IF
653:                         D(p) = D(p) * CS
654:                         D(q) = D(q) / CS
655:                      END IF
656:                   ELSE
657:                      IF ( D(q) .GE. ONE ) THEN
658:                         CALL DAXPY( M,     T*APOAQ, A(1,p),1, A(1,q),1 )
659:                         CALL DAXPY( M,-CS*SN*AQOAP, A(1,q),1, A(1,p),1 )
660:                         IF ( RSVEC ) THEN
661:                         CALL DAXPY(MVL,T*APOAQ,     V(1,p),1, V(1,q),1 )
662:                         CALL DAXPY(MVL,-CS*SN*AQOAP,V(1,q),1, V(1,p),1 )
663:                         END IF
664:                         D(p) = D(p) / CS
665:                         D(q) = D(q) * CS
666:                      ELSE
667:                         IF ( D(p) .GE. D(q) ) THEN
668:                            CALL DAXPY( M,-T*AQOAP,   A(1,q),1,A(1,p),1 )
669:                            CALL DAXPY( M,CS*SN*APOAQ,A(1,p),1,A(1,q),1 )
670:                            D(p) = D(p) * CS
671:                            D(q) = D(q) / CS
672:                            IF ( RSVEC ) THEN
673:                            CALL DAXPY( MVL, -T*AQOAP, V(1,q),1,V(1,p),1)
674:                            CALL DAXPY(MVL,CS*SN*APOAQ,V(1,p),1,V(1,q),1)
675:                            END IF
676:                         ELSE
677:                            CALL DAXPY(M, T*APOAQ,    A(1,p),1,A(1,q),1)
678:                            CALL DAXPY(M,-CS*SN*AQOAP,A(1,q),1,A(1,p),1)
679:                            D(p) = D(p) / CS
680:                            D(q) = D(q) * CS
681:                           IF ( RSVEC ) THEN
682:                           CALL DAXPY(MVL, T*APOAQ,    V(1,p),1,V(1,q),1)
683:                           CALL DAXPY(MVL,-CS*SN*AQOAP,V(1,q),1,V(1,p),1)
684:                           END IF
685:                         END IF
686:                      END IF
687:                   ENDIF
688:                END IF
689: *
690:             ELSE
691:                IF ( AAPP .GT. AAQQ ) THEN
692:                   CALL DCOPY( M, A(1,p), 1, WORK, 1 )
693:                   CALL DLASCL('G',0,0,AAPP,ONE,M,1,WORK,LDA,IERR)
694:                   CALL DLASCL('G',0,0,AAQQ,ONE,M,1,   A(1,q),LDA,IERR)
695:                   TEMP1 = -AAPQ * D(p) / D(q)
696:                   CALL DAXPY(M,TEMP1,WORK,1,A(1,q),1)
697:                   CALL DLASCL('G',0,0,ONE,AAQQ,M,1,A(1,q),LDA,IERR)
698:                   SVA(q) = AAQQ*DSQRT(DMAX1(ZERO, ONE - AAPQ*AAPQ))
699:                   MXSINJ = DMAX1( MXSINJ, SFMIN )
700:                ELSE
701:                   CALL DCOPY( M, A(1,q), 1, WORK, 1 )
702:                   CALL DLASCL('G',0,0,AAQQ,ONE,M,1,WORK,LDA,IERR)
703:                   CALL DLASCL('G',0,0,AAPP,ONE,M,1,   A(1,p),LDA,IERR)
704:                   TEMP1 = -AAPQ * D(q) / D(p)
705:                   CALL DAXPY(M,TEMP1,WORK,1,A(1,p),1)
706:                   CALL DLASCL('G',0,0,ONE,AAPP,M,1,A(1,p),LDA,IERR)
707:                   SVA(p) = AAPP*DSQRT(DMAX1(ZERO, ONE - AAPQ*AAPQ))
708:                   MXSINJ = DMAX1( MXSINJ, SFMIN )
709:                END IF
710:             END IF
711: *           END IF ROTOK THEN ... ELSE
712: *
713: *           In the case of cancellation in updating SVA(q)
714: *           .. recompute SVA(q)
715:             IF ( (SVA(q) / AAQQ )**2 .LE. ROOTEPS  ) THEN
716:                IF ((AAQQ .LT. ROOTBIG).AND.(AAQQ .GT. ROOTSFMIN)) THEN
717:                   SVA(q) = DNRM2( M, A(1,q), 1 ) * D(q)
718:                ELSE
719:                   T    = ZERO
720:                   AAQQ = ZERO
721:                   CALL DLASSQ( M, A(1,q), 1, T, AAQQ )
722:                   SVA(q) = T * DSQRT(AAQQ) * D(q)
723:                END IF
724:             END IF
725:             IF ( (AAPP / AAPP0 )**2 .LE. ROOTEPS  ) THEN
726:                IF ((AAPP .LT. ROOTBIG).AND.(AAPP .GT. ROOTSFMIN)) THEN
727:                   AAPP = DNRM2( M, A(1,p), 1 ) * D(p)
728:                ELSE
729:                   T    = ZERO
730:                   AAPP = ZERO
731:                   CALL DLASSQ( M, A(1,p), 1, T, AAPP )
732:                   AAPP = T * DSQRT(AAPP) * D(p)
733:                END IF
734:                SVA(p) = AAPP
735:             END IF
736: *              end of OK rotation
737:          ELSE
738:             NOTROT   = NOTROT   + 1
739:             PSKIPPED = PSKIPPED + 1
740:             IJBLSK   = IJBLSK   + 1
741:          END IF
742:       ELSE
743:          NOTROT   = NOTROT   + 1
744:          PSKIPPED = PSKIPPED + 1
745:          IJBLSK   = IJBLSK   + 1
746:       END IF
747: *
748:       IF ( ( i .LE. SWBAND ) .AND. ( IJBLSK .GE. BLSKIP ) ) THEN
749:          SVA(p) = AAPP
750:          NOTROT = 0
751:          GO TO 2011
752:       END IF
753:       IF ( ( i .LE. SWBAND ) .AND. ( PSKIPPED .GT. ROWSKIP ) ) THEN
754:          AAPP = -AAPP
755:          NOTROT = 0
756:          GO TO 2203
757:       END IF
758: *
759:  2200    CONTINUE
760: *        end of the q-loop
761:  2203    CONTINUE
762: *
763:          SVA(p) = AAPP
764: *
765:       ELSE
766:          IF ( AAPP .EQ. ZERO ) NOTROT=NOTROT+MIN0(jgl+KBL-1,N)-jgl+1
767:          IF ( AAPP .LT. ZERO ) NOTROT = 0
768:       END IF
769: 
770:  2100 CONTINUE
771: *     end of the p-loop
772:  2010 CONTINUE
773: *     end of the jbc-loop
774:  2011 CONTINUE
775: *2011 bailed out of the jbc-loop
776:       DO 2012 p = igl, MIN0( igl + KBL - 1, N )
777:          SVA(p) = DABS(SVA(p))
778:  2012 CONTINUE
779: *
780:  2000 CONTINUE
781: *2000 :: end of the ibr-loop
782: *
783: *     .. update SVA(N)
784:       IF ((SVA(N) .LT. ROOTBIG).AND.(SVA(N) .GT. ROOTSFMIN)) THEN
785:          SVA(N) = DNRM2( M, A(1,N), 1 ) * D(N)
786:       ELSE
787:          T    = ZERO
788:          AAPP = ZERO
789:          CALL DLASSQ( M, A(1,N), 1, T, AAPP )
790:          SVA(N) = T * DSQRT(AAPP) * D(N)
791:       END IF
792: *
793: *     Additional steering devices
794: *
795:       IF ( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
796:      &     ( ISWROT .LE. N ) ) )
797:      &   SWBAND = i
798: *
799:       IF ((i.GT.SWBAND+1).AND. (MXAAPQ.LT.DBLE(N)*TOL).AND.
800:      &   (DBLE(N)*MXAAPQ*MXSINJ.LT.TOL))THEN
801:         GO TO 1994
802:       END IF
803: *
804:       IF ( NOTROT .GE. EMPTSW ) GO TO 1994
805: 
806:  1993 CONTINUE
807: *     end i=1:NSWEEP loop
808: * #:) Reaching this point means that the procedure has comleted the given
809: *     number of iterations.
810:       INFO = NSWEEP - 1
811:       GO TO 1995
812:  1994 CONTINUE
813: * #:) Reaching this point means that during the i-th sweep all pivots were
814: *     below the given tolerance, causing early exit.
815: *
816:       INFO = 0
817: * #:) INFO = 0 confirms successful iterations.
818:  1995 CONTINUE
819: *
820: *     Sort the vector D.
821:       DO 5991 p = 1, N - 1
822:          q = IDAMAX( N-p+1, SVA(p), 1 ) + p - 1
823:          IF ( p .NE. q ) THEN
824:             TEMP1  = SVA(p)
825:             SVA(p) = SVA(q)
826:             SVA(q) = TEMP1
827:             TEMP1   = D(p)
828:             D(p) = D(q)
829:             D(q) = TEMP1
830:             CALL DSWAP( M, A(1,p), 1, A(1,q), 1 )
831:             IF ( RSVEC ) CALL DSWAP( MVL, V(1,p), 1, V(1,q), 1 )
832:          END IF
833:  5991 CONTINUE
834: *
835:       RETURN
836: *     ..
837: *     .. END OF DGSVJ0
838: *     ..
839:       END
840: *
841: