1: *> \brief \b DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASD7 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
22: * VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
23: * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
24: * C, S, INFO )
25: *
26: * .. Scalar Arguments ..
27: * INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
28: * $ NR, SQRE
29: * DOUBLE PRECISION ALPHA, BETA, C, S
30: * ..
31: * .. Array Arguments ..
32: * INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
33: * $ IDXQ( * ), PERM( * )
34: * DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
35: * $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
36: * $ ZW( * )
37: * ..
38: *
39: *
40: *> \par Purpose:
41: * =============
42: *>
43: *> \verbatim
44: *>
45: *> DLASD7 merges the two sets of singular values together into a single
46: *> sorted set. Then it tries to deflate the size of the problem. There
47: *> are two ways in which deflation can occur: when two or more singular
48: *> values are close together or if there is a tiny entry in the Z
49: *> vector. For each such occurrence the order of the related
50: *> secular equation problem is reduced by one.
51: *>
52: *> DLASD7 is called from DLASD6.
53: *> \endverbatim
54: *
55: * Arguments:
56: * ==========
57: *
58: *> \param[in] ICOMPQ
59: *> \verbatim
60: *> ICOMPQ is INTEGER
61: *> Specifies whether singular vectors are to be computed
62: *> in compact form, as follows:
63: *> = 0: Compute singular values only.
64: *> = 1: Compute singular vectors of upper
65: *> bidiagonal matrix in compact form.
66: *> \endverbatim
67: *>
68: *> \param[in] NL
69: *> \verbatim
70: *> NL is INTEGER
71: *> The row dimension of the upper block. NL >= 1.
72: *> \endverbatim
73: *>
74: *> \param[in] NR
75: *> \verbatim
76: *> NR is INTEGER
77: *> The row dimension of the lower block. NR >= 1.
78: *> \endverbatim
79: *>
80: *> \param[in] SQRE
81: *> \verbatim
82: *> SQRE is INTEGER
83: *> = 0: the lower block is an NR-by-NR square matrix.
84: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
85: *>
86: *> The bidiagonal matrix has
87: *> N = NL + NR + 1 rows and
88: *> M = N + SQRE >= N columns.
89: *> \endverbatim
90: *>
91: *> \param[out] K
92: *> \verbatim
93: *> K is INTEGER
94: *> Contains the dimension of the non-deflated matrix, this is
95: *> the order of the related secular equation. 1 <= K <=N.
96: *> \endverbatim
97: *>
98: *> \param[in,out] D
99: *> \verbatim
100: *> D is DOUBLE PRECISION array, dimension ( N )
101: *> On entry D contains the singular values of the two submatrices
102: *> to be combined. On exit D contains the trailing (N-K) updated
103: *> singular values (those which were deflated) sorted into
104: *> increasing order.
105: *> \endverbatim
106: *>
107: *> \param[out] Z
108: *> \verbatim
109: *> Z is DOUBLE PRECISION array, dimension ( M )
110: *> On exit Z contains the updating row vector in the secular
111: *> equation.
112: *> \endverbatim
113: *>
114: *> \param[out] ZW
115: *> \verbatim
116: *> ZW is DOUBLE PRECISION array, dimension ( M )
117: *> Workspace for Z.
118: *> \endverbatim
119: *>
120: *> \param[in,out] VF
121: *> \verbatim
122: *> VF is DOUBLE PRECISION array, dimension ( M )
123: *> On entry, VF(1:NL+1) contains the first components of all
124: *> right singular vectors of the upper block; and VF(NL+2:M)
125: *> contains the first components of all right singular vectors
126: *> of the lower block. On exit, VF contains the first components
127: *> of all right singular vectors of the bidiagonal matrix.
128: *> \endverbatim
129: *>
130: *> \param[out] VFW
131: *> \verbatim
132: *> VFW is DOUBLE PRECISION array, dimension ( M )
133: *> Workspace for VF.
134: *> \endverbatim
135: *>
136: *> \param[in,out] VL
137: *> \verbatim
138: *> VL is DOUBLE PRECISION array, dimension ( M )
139: *> On entry, VL(1:NL+1) contains the last components of all
140: *> right singular vectors of the upper block; and VL(NL+2:M)
141: *> contains the last components of all right singular vectors
142: *> of the lower block. On exit, VL contains the last components
143: *> of all right singular vectors of the bidiagonal matrix.
144: *> \endverbatim
145: *>
146: *> \param[out] VLW
147: *> \verbatim
148: *> VLW is DOUBLE PRECISION array, dimension ( M )
149: *> Workspace for VL.
150: *> \endverbatim
151: *>
152: *> \param[in] ALPHA
153: *> \verbatim
154: *> ALPHA is DOUBLE PRECISION
155: *> Contains the diagonal element associated with the added row.
156: *> \endverbatim
157: *>
158: *> \param[in] BETA
159: *> \verbatim
160: *> BETA is DOUBLE PRECISION
161: *> Contains the off-diagonal element associated with the added
162: *> row.
163: *> \endverbatim
164: *>
165: *> \param[out] DSIGMA
166: *> \verbatim
167: *> DSIGMA is DOUBLE PRECISION array, dimension ( N )
168: *> Contains a copy of the diagonal elements (K-1 singular values
169: *> and one zero) in the secular equation.
170: *> \endverbatim
171: *>
172: *> \param[out] IDX
173: *> \verbatim
174: *> IDX is INTEGER array, dimension ( N )
175: *> This will contain the permutation used to sort the contents of
176: *> D into ascending order.
177: *> \endverbatim
178: *>
179: *> \param[out] IDXP
180: *> \verbatim
181: *> IDXP is INTEGER array, dimension ( N )
182: *> This will contain the permutation used to place deflated
183: *> values of D at the end of the array. On output IDXP(2:K)
184: *> points to the nondeflated D-values and IDXP(K+1:N)
185: *> points to the deflated singular values.
186: *> \endverbatim
187: *>
188: *> \param[in] IDXQ
189: *> \verbatim
190: *> IDXQ is INTEGER array, dimension ( N )
191: *> This contains the permutation which separately sorts the two
192: *> sub-problems in D into ascending order. Note that entries in
193: *> the first half of this permutation must first be moved one
194: *> position backward; and entries in the second half
195: *> must first have NL+1 added to their values.
196: *> \endverbatim
197: *>
198: *> \param[out] PERM
199: *> \verbatim
200: *> PERM is INTEGER array, dimension ( N )
201: *> The permutations (from deflation and sorting) to be applied
202: *> to each singular block. Not referenced if ICOMPQ = 0.
203: *> \endverbatim
204: *>
205: *> \param[out] GIVPTR
206: *> \verbatim
207: *> GIVPTR is INTEGER
208: *> The number of Givens rotations which took place in this
209: *> subproblem. Not referenced if ICOMPQ = 0.
210: *> \endverbatim
211: *>
212: *> \param[out] GIVCOL
213: *> \verbatim
214: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
215: *> Each pair of numbers indicates a pair of columns to take place
216: *> in a Givens rotation. Not referenced if ICOMPQ = 0.
217: *> \endverbatim
218: *>
219: *> \param[in] LDGCOL
220: *> \verbatim
221: *> LDGCOL is INTEGER
222: *> The leading dimension of GIVCOL, must be at least N.
223: *> \endverbatim
224: *>
225: *> \param[out] GIVNUM
226: *> \verbatim
227: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
228: *> Each number indicates the C or S value to be used in the
229: *> corresponding Givens rotation. Not referenced if ICOMPQ = 0.
230: *> \endverbatim
231: *>
232: *> \param[in] LDGNUM
233: *> \verbatim
234: *> LDGNUM is INTEGER
235: *> The leading dimension of GIVNUM, must be at least N.
236: *> \endverbatim
237: *>
238: *> \param[out] C
239: *> \verbatim
240: *> C is DOUBLE PRECISION
241: *> C contains garbage if SQRE =0 and the C-value of a Givens
242: *> rotation related to the right null space if SQRE = 1.
243: *> \endverbatim
244: *>
245: *> \param[out] S
246: *> \verbatim
247: *> S is DOUBLE PRECISION
248: *> S contains garbage if SQRE =0 and the S-value of a Givens
249: *> rotation related to the right null space if SQRE = 1.
250: *> \endverbatim
251: *>
252: *> \param[out] INFO
253: *> \verbatim
254: *> INFO is INTEGER
255: *> = 0: successful exit.
256: *> < 0: if INFO = -i, the i-th argument had an illegal value.
257: *> \endverbatim
258: *
259: * Authors:
260: * ========
261: *
262: *> \author Univ. of Tennessee
263: *> \author Univ. of California Berkeley
264: *> \author Univ. of Colorado Denver
265: *> \author NAG Ltd.
266: *
267: *> \ingroup OTHERauxiliary
268: *
269: *> \par Contributors:
270: * ==================
271: *>
272: *> Ming Gu and Huan Ren, Computer Science Division, University of
273: *> California at Berkeley, USA
274: *>
275: * =====================================================================
276: SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
277: $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
278: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
279: $ C, S, INFO )
280: *
281: * -- LAPACK auxiliary routine --
282: * -- LAPACK is a software package provided by Univ. of Tennessee, --
283: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
284: *
285: * .. Scalar Arguments ..
286: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
287: $ NR, SQRE
288: DOUBLE PRECISION ALPHA, BETA, C, S
289: * ..
290: * .. Array Arguments ..
291: INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
292: $ IDXQ( * ), PERM( * )
293: DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
294: $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
295: $ ZW( * )
296: * ..
297: *
298: * =====================================================================
299: *
300: * .. Parameters ..
301: DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
302: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
303: $ EIGHT = 8.0D+0 )
304: * ..
305: * .. Local Scalars ..
306: *
307: INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
308: $ NLP1, NLP2
309: DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
310: * ..
311: * .. External Subroutines ..
312: EXTERNAL DCOPY, DLAMRG, DROT, XERBLA
313: * ..
314: * .. External Functions ..
315: DOUBLE PRECISION DLAMCH, DLAPY2
316: EXTERNAL DLAMCH, DLAPY2
317: * ..
318: * .. Intrinsic Functions ..
319: INTRINSIC ABS, MAX
320: * ..
321: * .. Executable Statements ..
322: *
323: * Test the input parameters.
324: *
325: INFO = 0
326: N = NL + NR + 1
327: M = N + SQRE
328: *
329: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
330: INFO = -1
331: ELSE IF( NL.LT.1 ) THEN
332: INFO = -2
333: ELSE IF( NR.LT.1 ) THEN
334: INFO = -3
335: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
336: INFO = -4
337: ELSE IF( LDGCOL.LT.N ) THEN
338: INFO = -22
339: ELSE IF( LDGNUM.LT.N ) THEN
340: INFO = -24
341: END IF
342: IF( INFO.NE.0 ) THEN
343: CALL XERBLA( 'DLASD7', -INFO )
344: RETURN
345: END IF
346: *
347: NLP1 = NL + 1
348: NLP2 = NL + 2
349: IF( ICOMPQ.EQ.1 ) THEN
350: GIVPTR = 0
351: END IF
352: *
353: * Generate the first part of the vector Z and move the singular
354: * values in the first part of D one position backward.
355: *
356: Z1 = ALPHA*VL( NLP1 )
357: VL( NLP1 ) = ZERO
358: TAU = VF( NLP1 )
359: DO 10 I = NL, 1, -1
360: Z( I+1 ) = ALPHA*VL( I )
361: VL( I ) = ZERO
362: VF( I+1 ) = VF( I )
363: D( I+1 ) = D( I )
364: IDXQ( I+1 ) = IDXQ( I ) + 1
365: 10 CONTINUE
366: VF( 1 ) = TAU
367: *
368: * Generate the second part of the vector Z.
369: *
370: DO 20 I = NLP2, M
371: Z( I ) = BETA*VF( I )
372: VF( I ) = ZERO
373: 20 CONTINUE
374: *
375: * Sort the singular values into increasing order
376: *
377: DO 30 I = NLP2, N
378: IDXQ( I ) = IDXQ( I ) + NLP1
379: 30 CONTINUE
380: *
381: * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
382: *
383: DO 40 I = 2, N
384: DSIGMA( I ) = D( IDXQ( I ) )
385: ZW( I ) = Z( IDXQ( I ) )
386: VFW( I ) = VF( IDXQ( I ) )
387: VLW( I ) = VL( IDXQ( I ) )
388: 40 CONTINUE
389: *
390: CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
391: *
392: DO 50 I = 2, N
393: IDXI = 1 + IDX( I )
394: D( I ) = DSIGMA( IDXI )
395: Z( I ) = ZW( IDXI )
396: VF( I ) = VFW( IDXI )
397: VL( I ) = VLW( IDXI )
398: 50 CONTINUE
399: *
400: * Calculate the allowable deflation tolerance
401: *
402: EPS = DLAMCH( 'Epsilon' )
403: TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
404: TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
405: *
406: * There are 2 kinds of deflation -- first a value in the z-vector
407: * is small, second two (or more) singular values are very close
408: * together (their difference is small).
409: *
410: * If the value in the z-vector is small, we simply permute the
411: * array so that the corresponding singular value is moved to the
412: * end.
413: *
414: * If two values in the D-vector are close, we perform a two-sided
415: * rotation designed to make one of the corresponding z-vector
416: * entries zero, and then permute the array so that the deflated
417: * singular value is moved to the end.
418: *
419: * If there are multiple singular values then the problem deflates.
420: * Here the number of equal singular values are found. As each equal
421: * singular value is found, an elementary reflector is computed to
422: * rotate the corresponding singular subspace so that the
423: * corresponding components of Z are zero in this new basis.
424: *
425: K = 1
426: K2 = N + 1
427: DO 60 J = 2, N
428: IF( ABS( Z( J ) ).LE.TOL ) THEN
429: *
430: * Deflate due to small z component.
431: *
432: K2 = K2 - 1
433: IDXP( K2 ) = J
434: IF( J.EQ.N )
435: $ GO TO 100
436: ELSE
437: JPREV = J
438: GO TO 70
439: END IF
440: 60 CONTINUE
441: 70 CONTINUE
442: J = JPREV
443: 80 CONTINUE
444: J = J + 1
445: IF( J.GT.N )
446: $ GO TO 90
447: IF( ABS( Z( J ) ).LE.TOL ) THEN
448: *
449: * Deflate due to small z component.
450: *
451: K2 = K2 - 1
452: IDXP( K2 ) = J
453: ELSE
454: *
455: * Check if singular values are close enough to allow deflation.
456: *
457: IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
458: *
459: * Deflation is possible.
460: *
461: S = Z( JPREV )
462: C = Z( J )
463: *
464: * Find sqrt(a**2+b**2) without overflow or
465: * destructive underflow.
466: *
467: TAU = DLAPY2( C, S )
468: Z( J ) = TAU
469: Z( JPREV ) = ZERO
470: C = C / TAU
471: S = -S / TAU
472: *
473: * Record the appropriate Givens rotation
474: *
475: IF( ICOMPQ.EQ.1 ) THEN
476: GIVPTR = GIVPTR + 1
477: IDXJP = IDXQ( IDX( JPREV )+1 )
478: IDXJ = IDXQ( IDX( J )+1 )
479: IF( IDXJP.LE.NLP1 ) THEN
480: IDXJP = IDXJP - 1
481: END IF
482: IF( IDXJ.LE.NLP1 ) THEN
483: IDXJ = IDXJ - 1
484: END IF
485: GIVCOL( GIVPTR, 2 ) = IDXJP
486: GIVCOL( GIVPTR, 1 ) = IDXJ
487: GIVNUM( GIVPTR, 2 ) = C
488: GIVNUM( GIVPTR, 1 ) = S
489: END IF
490: CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
491: CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
492: K2 = K2 - 1
493: IDXP( K2 ) = JPREV
494: JPREV = J
495: ELSE
496: K = K + 1
497: ZW( K ) = Z( JPREV )
498: DSIGMA( K ) = D( JPREV )
499: IDXP( K ) = JPREV
500: JPREV = J
501: END IF
502: END IF
503: GO TO 80
504: 90 CONTINUE
505: *
506: * Record the last singular value.
507: *
508: K = K + 1
509: ZW( K ) = Z( JPREV )
510: DSIGMA( K ) = D( JPREV )
511: IDXP( K ) = JPREV
512: *
513: 100 CONTINUE
514: *
515: * Sort the singular values into DSIGMA. The singular values which
516: * were not deflated go into the first K slots of DSIGMA, except
517: * that DSIGMA(1) is treated separately.
518: *
519: DO 110 J = 2, N
520: JP = IDXP( J )
521: DSIGMA( J ) = D( JP )
522: VFW( J ) = VF( JP )
523: VLW( J ) = VL( JP )
524: 110 CONTINUE
525: IF( ICOMPQ.EQ.1 ) THEN
526: DO 120 J = 2, N
527: JP = IDXP( J )
528: PERM( J ) = IDXQ( IDX( JP )+1 )
529: IF( PERM( J ).LE.NLP1 ) THEN
530: PERM( J ) = PERM( J ) - 1
531: END IF
532: 120 CONTINUE
533: END IF
534: *
535: * The deflated singular values go back into the last N - K slots of
536: * D.
537: *
538: CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
539: *
540: * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
541: * VL(M).
542: *
543: DSIGMA( 1 ) = ZERO
544: HLFTOL = TOL / TWO
545: IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
546: $ DSIGMA( 2 ) = HLFTOL
547: IF( M.GT.N ) THEN
548: Z( 1 ) = DLAPY2( Z1, Z( M ) )
549: IF( Z( 1 ).LE.TOL ) THEN
550: C = ONE
551: S = ZERO
552: Z( 1 ) = TOL
553: ELSE
554: C = Z1 / Z( 1 )
555: S = -Z( M ) / Z( 1 )
556: END IF
557: CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
558: CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
559: ELSE
560: IF( ABS( Z1 ).LE.TOL ) THEN
561: Z( 1 ) = TOL
562: ELSE
563: Z( 1 ) = Z1
564: END IF
565: END IF
566: *
567: * Restore Z, VF, and VL.
568: *
569: CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
570: CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
571: CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
572: *
573: RETURN
574: *
575: * End of DLASD7
576: *
577: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd7.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack