001:       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
002:      $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
003:      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
004:      $                   C, S, INFO )
005: *
006: *  -- LAPACK auxiliary routine (version 3.2) --
007: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
012:      $                   NR, SQRE
013:       DOUBLE PRECISION   ALPHA, BETA, C, S
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
017:      $                   IDXQ( * ), PERM( * )
018:       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
019:      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
020:      $                   ZW( * )
021: *     ..
022: *
023: *  Purpose
024: *  =======
025: *
026: *  DLASD7 merges the two sets of singular values together into a single
027: *  sorted set. Then it tries to deflate the size of the problem. There
028: *  are two ways in which deflation can occur:  when two or more singular
029: *  values are close together or if there is a tiny entry in the Z
030: *  vector. For each such occurrence the order of the related
031: *  secular equation problem is reduced by one.
032: *
033: *  DLASD7 is called from DLASD6.
034: *
035: *  Arguments
036: *  =========
037: *
038: *  ICOMPQ  (input) INTEGER
039: *          Specifies whether singular vectors are to be computed
040: *          in compact form, as follows:
041: *          = 0: Compute singular values only.
042: *          = 1: Compute singular vectors of upper
043: *               bidiagonal matrix in compact form.
044: *
045: *  NL     (input) INTEGER
046: *         The row dimension of the upper block. NL >= 1.
047: *
048: *  NR     (input) INTEGER
049: *         The row dimension of the lower block. NR >= 1.
050: *
051: *  SQRE   (input) INTEGER
052: *         = 0: the lower block is an NR-by-NR square matrix.
053: *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
054: *
055: *         The bidiagonal matrix has
056: *         N = NL + NR + 1 rows and
057: *         M = N + SQRE >= N columns.
058: *
059: *  K      (output) INTEGER
060: *         Contains the dimension of the non-deflated matrix, this is
061: *         the order of the related secular equation. 1 <= K <=N.
062: *
063: *  D      (input/output) DOUBLE PRECISION array, dimension ( N )
064: *         On entry D contains the singular values of the two submatrices
065: *         to be combined. On exit D contains the trailing (N-K) updated
066: *         singular values (those which were deflated) sorted into
067: *         increasing order.
068: *
069: *  Z      (output) DOUBLE PRECISION array, dimension ( M )
070: *         On exit Z contains the updating row vector in the secular
071: *         equation.
072: *
073: *  ZW     (workspace) DOUBLE PRECISION array, dimension ( M )
074: *         Workspace for Z.
075: *
076: *  VF     (input/output) DOUBLE PRECISION array, dimension ( M )
077: *         On entry, VF(1:NL+1) contains the first components of all
078: *         right singular vectors of the upper block; and VF(NL+2:M)
079: *         contains the first components of all right singular vectors
080: *         of the lower block. On exit, VF contains the first components
081: *         of all right singular vectors of the bidiagonal matrix.
082: *
083: *  VFW    (workspace) DOUBLE PRECISION array, dimension ( M )
084: *         Workspace for VF.
085: *
086: *  VL     (input/output) DOUBLE PRECISION array, dimension ( M )
087: *         On entry, VL(1:NL+1) contains the  last components of all
088: *         right singular vectors of the upper block; and VL(NL+2:M)
089: *         contains the last components of all right singular vectors
090: *         of the lower block. On exit, VL contains the last components
091: *         of all right singular vectors of the bidiagonal matrix.
092: *
093: *  VLW    (workspace) DOUBLE PRECISION array, dimension ( M )
094: *         Workspace for VL.
095: *
096: *  ALPHA  (input) DOUBLE PRECISION
097: *         Contains the diagonal element associated with the added row.
098: *
099: *  BETA   (input) DOUBLE PRECISION
100: *         Contains the off-diagonal element associated with the added
101: *         row.
102: *
103: *  DSIGMA (output) DOUBLE PRECISION array, dimension ( N )
104: *         Contains a copy of the diagonal elements (K-1 singular values
105: *         and one zero) in the secular equation.
106: *
107: *  IDX    (workspace) INTEGER array, dimension ( N )
108: *         This will contain the permutation used to sort the contents of
109: *         D into ascending order.
110: *
111: *  IDXP   (workspace) INTEGER array, dimension ( N )
112: *         This will contain the permutation used to place deflated
113: *         values of D at the end of the array. On output IDXP(2:K)
114: *         points to the nondeflated D-values and IDXP(K+1:N)
115: *         points to the deflated singular values.
116: *
117: *  IDXQ   (input) INTEGER array, dimension ( N )
118: *         This contains the permutation which separately sorts the two
119: *         sub-problems in D into ascending order.  Note that entries in
120: *         the first half of this permutation must first be moved one
121: *         position backward; and entries in the second half
122: *         must first have NL+1 added to their values.
123: *
124: *  PERM   (output) INTEGER array, dimension ( N )
125: *         The permutations (from deflation and sorting) to be applied
126: *         to each singular block. Not referenced if ICOMPQ = 0.
127: *
128: *  GIVPTR (output) INTEGER
129: *         The number of Givens rotations which took place in this
130: *         subproblem. Not referenced if ICOMPQ = 0.
131: *
132: *  GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 )
133: *         Each pair of numbers indicates a pair of columns to take place
134: *         in a Givens rotation. Not referenced if ICOMPQ = 0.
135: *
136: *  LDGCOL (input) INTEGER
137: *         The leading dimension of GIVCOL, must be at least N.
138: *
139: *  GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
140: *         Each number indicates the C or S value to be used in the
141: *         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
142: *
143: *  LDGNUM (input) INTEGER
144: *         The leading dimension of GIVNUM, must be at least N.
145: *
146: *  C      (output) DOUBLE PRECISION
147: *         C contains garbage if SQRE =0 and the C-value of a Givens
148: *         rotation related to the right null space if SQRE = 1.
149: *
150: *  S      (output) DOUBLE PRECISION
151: *         S contains garbage if SQRE =0 and the S-value of a Givens
152: *         rotation related to the right null space if SQRE = 1.
153: *
154: *  INFO   (output) INTEGER
155: *         = 0:  successful exit.
156: *         < 0:  if INFO = -i, the i-th argument had an illegal value.
157: *
158: *  Further Details
159: *  ===============
160: *
161: *  Based on contributions by
162: *     Ming Gu and Huan Ren, Computer Science Division, University of
163: *     California at Berkeley, USA
164: *
165: *  =====================================================================
166: *
167: *     .. Parameters ..
168:       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
169:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
170:      $                   EIGHT = 8.0D+0 )
171: *     ..
172: *     .. Local Scalars ..
173: *
174:       INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
175:      $                   NLP1, NLP2
176:       DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
177: *     ..
178: *     .. External Subroutines ..
179:       EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA
180: *     ..
181: *     .. External Functions ..
182:       DOUBLE PRECISION   DLAMCH, DLAPY2
183:       EXTERNAL           DLAMCH, DLAPY2
184: *     ..
185: *     .. Intrinsic Functions ..
186:       INTRINSIC          ABS, MAX
187: *     ..
188: *     .. Executable Statements ..
189: *
190: *     Test the input parameters.
191: *
192:       INFO = 0
193:       N = NL + NR + 1
194:       M = N + SQRE
195: *
196:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
197:          INFO = -1
198:       ELSE IF( NL.LT.1 ) THEN
199:          INFO = -2
200:       ELSE IF( NR.LT.1 ) THEN
201:          INFO = -3
202:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
203:          INFO = -4
204:       ELSE IF( LDGCOL.LT.N ) THEN
205:          INFO = -22
206:       ELSE IF( LDGNUM.LT.N ) THEN
207:          INFO = -24
208:       END IF
209:       IF( INFO.NE.0 ) THEN
210:          CALL XERBLA( 'DLASD7', -INFO )
211:          RETURN
212:       END IF
213: *
214:       NLP1 = NL + 1
215:       NLP2 = NL + 2
216:       IF( ICOMPQ.EQ.1 ) THEN
217:          GIVPTR = 0
218:       END IF
219: *
220: *     Generate the first part of the vector Z and move the singular
221: *     values in the first part of D one position backward.
222: *
223:       Z1 = ALPHA*VL( NLP1 )
224:       VL( NLP1 ) = ZERO
225:       TAU = VF( NLP1 )
226:       DO 10 I = NL, 1, -1
227:          Z( I+1 ) = ALPHA*VL( I )
228:          VL( I ) = ZERO
229:          VF( I+1 ) = VF( I )
230:          D( I+1 ) = D( I )
231:          IDXQ( I+1 ) = IDXQ( I ) + 1
232:    10 CONTINUE
233:       VF( 1 ) = TAU
234: *
235: *     Generate the second part of the vector Z.
236: *
237:       DO 20 I = NLP2, M
238:          Z( I ) = BETA*VF( I )
239:          VF( I ) = ZERO
240:    20 CONTINUE
241: *
242: *     Sort the singular values into increasing order
243: *
244:       DO 30 I = NLP2, N
245:          IDXQ( I ) = IDXQ( I ) + NLP1
246:    30 CONTINUE
247: *
248: *     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
249: *
250:       DO 40 I = 2, N
251:          DSIGMA( I ) = D( IDXQ( I ) )
252:          ZW( I ) = Z( IDXQ( I ) )
253:          VFW( I ) = VF( IDXQ( I ) )
254:          VLW( I ) = VL( IDXQ( I ) )
255:    40 CONTINUE
256: *
257:       CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
258: *
259:       DO 50 I = 2, N
260:          IDXI = 1 + IDX( I )
261:          D( I ) = DSIGMA( IDXI )
262:          Z( I ) = ZW( IDXI )
263:          VF( I ) = VFW( IDXI )
264:          VL( I ) = VLW( IDXI )
265:    50 CONTINUE
266: *
267: *     Calculate the allowable deflation tolerence
268: *
269:       EPS = DLAMCH( 'Epsilon' )
270:       TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
271:       TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
272: *
273: *     There are 2 kinds of deflation -- first a value in the z-vector
274: *     is small, second two (or more) singular values are very close
275: *     together (their difference is small).
276: *
277: *     If the value in the z-vector is small, we simply permute the
278: *     array so that the corresponding singular value is moved to the
279: *     end.
280: *
281: *     If two values in the D-vector are close, we perform a two-sided
282: *     rotation designed to make one of the corresponding z-vector
283: *     entries zero, and then permute the array so that the deflated
284: *     singular value is moved to the end.
285: *
286: *     If there are multiple singular values then the problem deflates.
287: *     Here the number of equal singular values are found.  As each equal
288: *     singular value is found, an elementary reflector is computed to
289: *     rotate the corresponding singular subspace so that the
290: *     corresponding components of Z are zero in this new basis.
291: *
292:       K = 1
293:       K2 = N + 1
294:       DO 60 J = 2, N
295:          IF( ABS( Z( J ) ).LE.TOL ) THEN
296: *
297: *           Deflate due to small z component.
298: *
299:             K2 = K2 - 1
300:             IDXP( K2 ) = J
301:             IF( J.EQ.N )
302:      $         GO TO 100
303:          ELSE
304:             JPREV = J
305:             GO TO 70
306:          END IF
307:    60 CONTINUE
308:    70 CONTINUE
309:       J = JPREV
310:    80 CONTINUE
311:       J = J + 1
312:       IF( J.GT.N )
313:      $   GO TO 90
314:       IF( ABS( Z( J ) ).LE.TOL ) THEN
315: *
316: *        Deflate due to small z component.
317: *
318:          K2 = K2 - 1
319:          IDXP( K2 ) = J
320:       ELSE
321: *
322: *        Check if singular values are close enough to allow deflation.
323: *
324:          IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
325: *
326: *           Deflation is possible.
327: *
328:             S = Z( JPREV )
329:             C = Z( J )
330: *
331: *           Find sqrt(a**2+b**2) without overflow or
332: *           destructive underflow.
333: *
334:             TAU = DLAPY2( C, S )
335:             Z( J ) = TAU
336:             Z( JPREV ) = ZERO
337:             C = C / TAU
338:             S = -S / TAU
339: *
340: *           Record the appropriate Givens rotation
341: *
342:             IF( ICOMPQ.EQ.1 ) THEN
343:                GIVPTR = GIVPTR + 1
344:                IDXJP = IDXQ( IDX( JPREV )+1 )
345:                IDXJ = IDXQ( IDX( J )+1 )
346:                IF( IDXJP.LE.NLP1 ) THEN
347:                   IDXJP = IDXJP - 1
348:                END IF
349:                IF( IDXJ.LE.NLP1 ) THEN
350:                   IDXJ = IDXJ - 1
351:                END IF
352:                GIVCOL( GIVPTR, 2 ) = IDXJP
353:                GIVCOL( GIVPTR, 1 ) = IDXJ
354:                GIVNUM( GIVPTR, 2 ) = C
355:                GIVNUM( GIVPTR, 1 ) = S
356:             END IF
357:             CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
358:             CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
359:             K2 = K2 - 1
360:             IDXP( K2 ) = JPREV
361:             JPREV = J
362:          ELSE
363:             K = K + 1
364:             ZW( K ) = Z( JPREV )
365:             DSIGMA( K ) = D( JPREV )
366:             IDXP( K ) = JPREV
367:             JPREV = J
368:          END IF
369:       END IF
370:       GO TO 80
371:    90 CONTINUE
372: *
373: *     Record the last singular value.
374: *
375:       K = K + 1
376:       ZW( K ) = Z( JPREV )
377:       DSIGMA( K ) = D( JPREV )
378:       IDXP( K ) = JPREV
379: *
380:   100 CONTINUE
381: *
382: *     Sort the singular values into DSIGMA. The singular values which
383: *     were not deflated go into the first K slots of DSIGMA, except
384: *     that DSIGMA(1) is treated separately.
385: *
386:       DO 110 J = 2, N
387:          JP = IDXP( J )
388:          DSIGMA( J ) = D( JP )
389:          VFW( J ) = VF( JP )
390:          VLW( J ) = VL( JP )
391:   110 CONTINUE
392:       IF( ICOMPQ.EQ.1 ) THEN
393:          DO 120 J = 2, N
394:             JP = IDXP( J )
395:             PERM( J ) = IDXQ( IDX( JP )+1 )
396:             IF( PERM( J ).LE.NLP1 ) THEN
397:                PERM( J ) = PERM( J ) - 1
398:             END IF
399:   120    CONTINUE
400:       END IF
401: *
402: *     The deflated singular values go back into the last N - K slots of
403: *     D.
404: *
405:       CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
406: *
407: *     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
408: *     VL(M).
409: *
410:       DSIGMA( 1 ) = ZERO
411:       HLFTOL = TOL / TWO
412:       IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
413:      $   DSIGMA( 2 ) = HLFTOL
414:       IF( M.GT.N ) THEN
415:          Z( 1 ) = DLAPY2( Z1, Z( M ) )
416:          IF( Z( 1 ).LE.TOL ) THEN
417:             C = ONE
418:             S = ZERO
419:             Z( 1 ) = TOL
420:          ELSE
421:             C = Z1 / Z( 1 )
422:             S = -Z( M ) / Z( 1 )
423:          END IF
424:          CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
425:          CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
426:       ELSE
427:          IF( ABS( Z1 ).LE.TOL ) THEN
428:             Z( 1 ) = TOL
429:          ELSE
430:             Z( 1 ) = Z1
431:          END IF
432:       END IF
433: *
434: *     Restore Z, VF, and VL.
435: *
436:       CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
437:       CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
438:       CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
439: *
440:       RETURN
441: *
442: *     End of DLASD7
443: *
444:       END
445: