Annotation of rpl/lapack/lapack/dlasd7.f, revision 1.14
1.11 bertrand 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.
1.8 bertrand 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: *
1.11 bertrand 267: *> \date September 2012
1.8 bertrand 268: *
269: *> \ingroup auxOTHERauxiliary
270: *
271: *> \par Contributors:
272: * ==================
273: *>
274: *> Ming Gu and Huan Ren, Computer Science Division, University of
275: *> California at Berkeley, USA
276: *>
277: * =====================================================================
1.1 bertrand 278: SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
279: $ VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
280: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
281: $ C, S, INFO )
282: *
1.11 bertrand 283: * -- LAPACK auxiliary routine (version 3.4.2) --
1.1 bertrand 284: * -- LAPACK is a software package provided by Univ. of Tennessee, --
285: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 286: * September 2012
1.1 bertrand 287: *
288: * .. Scalar Arguments ..
289: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
290: $ NR, SQRE
291: DOUBLE PRECISION ALPHA, BETA, C, S
292: * ..
293: * .. Array Arguments ..
294: INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
295: $ IDXQ( * ), PERM( * )
296: DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
297: $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
298: $ ZW( * )
299: * ..
300: *
301: * =====================================================================
302: *
303: * .. Parameters ..
304: DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
305: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
306: $ EIGHT = 8.0D+0 )
307: * ..
308: * .. Local Scalars ..
309: *
310: INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
311: $ NLP1, NLP2
312: DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
313: * ..
314: * .. External Subroutines ..
315: EXTERNAL DCOPY, DLAMRG, DROT, XERBLA
316: * ..
317: * .. External Functions ..
318: DOUBLE PRECISION DLAMCH, DLAPY2
319: EXTERNAL DLAMCH, DLAPY2
320: * ..
321: * .. Intrinsic Functions ..
322: INTRINSIC ABS, MAX
323: * ..
324: * .. Executable Statements ..
325: *
326: * Test the input parameters.
327: *
328: INFO = 0
329: N = NL + NR + 1
330: M = N + SQRE
331: *
332: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
333: INFO = -1
334: ELSE IF( NL.LT.1 ) THEN
335: INFO = -2
336: ELSE IF( NR.LT.1 ) THEN
337: INFO = -3
338: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
339: INFO = -4
340: ELSE IF( LDGCOL.LT.N ) THEN
341: INFO = -22
342: ELSE IF( LDGNUM.LT.N ) THEN
343: INFO = -24
344: END IF
345: IF( INFO.NE.0 ) THEN
346: CALL XERBLA( 'DLASD7', -INFO )
347: RETURN
348: END IF
349: *
350: NLP1 = NL + 1
351: NLP2 = NL + 2
352: IF( ICOMPQ.EQ.1 ) THEN
353: GIVPTR = 0
354: END IF
355: *
356: * Generate the first part of the vector Z and move the singular
357: * values in the first part of D one position backward.
358: *
359: Z1 = ALPHA*VL( NLP1 )
360: VL( NLP1 ) = ZERO
361: TAU = VF( NLP1 )
362: DO 10 I = NL, 1, -1
363: Z( I+1 ) = ALPHA*VL( I )
364: VL( I ) = ZERO
365: VF( I+1 ) = VF( I )
366: D( I+1 ) = D( I )
367: IDXQ( I+1 ) = IDXQ( I ) + 1
368: 10 CONTINUE
369: VF( 1 ) = TAU
370: *
371: * Generate the second part of the vector Z.
372: *
373: DO 20 I = NLP2, M
374: Z( I ) = BETA*VF( I )
375: VF( I ) = ZERO
376: 20 CONTINUE
377: *
378: * Sort the singular values into increasing order
379: *
380: DO 30 I = NLP2, N
381: IDXQ( I ) = IDXQ( I ) + NLP1
382: 30 CONTINUE
383: *
384: * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
385: *
386: DO 40 I = 2, N
387: DSIGMA( I ) = D( IDXQ( I ) )
388: ZW( I ) = Z( IDXQ( I ) )
389: VFW( I ) = VF( IDXQ( I ) )
390: VLW( I ) = VL( IDXQ( I ) )
391: 40 CONTINUE
392: *
393: CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
394: *
395: DO 50 I = 2, N
396: IDXI = 1 + IDX( I )
397: D( I ) = DSIGMA( IDXI )
398: Z( I ) = ZW( IDXI )
399: VF( I ) = VFW( IDXI )
400: VL( I ) = VLW( IDXI )
401: 50 CONTINUE
402: *
403: * Calculate the allowable deflation tolerence
404: *
405: EPS = DLAMCH( 'Epsilon' )
406: TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
407: TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
408: *
409: * There are 2 kinds of deflation -- first a value in the z-vector
410: * is small, second two (or more) singular values are very close
411: * together (their difference is small).
412: *
413: * If the value in the z-vector is small, we simply permute the
414: * array so that the corresponding singular value is moved to the
415: * end.
416: *
417: * If two values in the D-vector are close, we perform a two-sided
418: * rotation designed to make one of the corresponding z-vector
419: * entries zero, and then permute the array so that the deflated
420: * singular value is moved to the end.
421: *
422: * If there are multiple singular values then the problem deflates.
423: * Here the number of equal singular values are found. As each equal
424: * singular value is found, an elementary reflector is computed to
425: * rotate the corresponding singular subspace so that the
426: * corresponding components of Z are zero in this new basis.
427: *
428: K = 1
429: K2 = N + 1
430: DO 60 J = 2, N
431: IF( ABS( Z( J ) ).LE.TOL ) THEN
432: *
433: * Deflate due to small z component.
434: *
435: K2 = K2 - 1
436: IDXP( K2 ) = J
437: IF( J.EQ.N )
438: $ GO TO 100
439: ELSE
440: JPREV = J
441: GO TO 70
442: END IF
443: 60 CONTINUE
444: 70 CONTINUE
445: J = JPREV
446: 80 CONTINUE
447: J = J + 1
448: IF( J.GT.N )
449: $ GO TO 90
450: IF( ABS( Z( J ) ).LE.TOL ) THEN
451: *
452: * Deflate due to small z component.
453: *
454: K2 = K2 - 1
455: IDXP( K2 ) = J
456: ELSE
457: *
458: * Check if singular values are close enough to allow deflation.
459: *
460: IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
461: *
462: * Deflation is possible.
463: *
464: S = Z( JPREV )
465: C = Z( J )
466: *
467: * Find sqrt(a**2+b**2) without overflow or
468: * destructive underflow.
469: *
470: TAU = DLAPY2( C, S )
471: Z( J ) = TAU
472: Z( JPREV ) = ZERO
473: C = C / TAU
474: S = -S / TAU
475: *
476: * Record the appropriate Givens rotation
477: *
478: IF( ICOMPQ.EQ.1 ) THEN
479: GIVPTR = GIVPTR + 1
480: IDXJP = IDXQ( IDX( JPREV )+1 )
481: IDXJ = IDXQ( IDX( J )+1 )
482: IF( IDXJP.LE.NLP1 ) THEN
483: IDXJP = IDXJP - 1
484: END IF
485: IF( IDXJ.LE.NLP1 ) THEN
486: IDXJ = IDXJ - 1
487: END IF
488: GIVCOL( GIVPTR, 2 ) = IDXJP
489: GIVCOL( GIVPTR, 1 ) = IDXJ
490: GIVNUM( GIVPTR, 2 ) = C
491: GIVNUM( GIVPTR, 1 ) = S
492: END IF
493: CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
494: CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
495: K2 = K2 - 1
496: IDXP( K2 ) = JPREV
497: JPREV = J
498: ELSE
499: K = K + 1
500: ZW( K ) = Z( JPREV )
501: DSIGMA( K ) = D( JPREV )
502: IDXP( K ) = JPREV
503: JPREV = J
504: END IF
505: END IF
506: GO TO 80
507: 90 CONTINUE
508: *
509: * Record the last singular value.
510: *
511: K = K + 1
512: ZW( K ) = Z( JPREV )
513: DSIGMA( K ) = D( JPREV )
514: IDXP( K ) = JPREV
515: *
516: 100 CONTINUE
517: *
518: * Sort the singular values into DSIGMA. The singular values which
519: * were not deflated go into the first K slots of DSIGMA, except
520: * that DSIGMA(1) is treated separately.
521: *
522: DO 110 J = 2, N
523: JP = IDXP( J )
524: DSIGMA( J ) = D( JP )
525: VFW( J ) = VF( JP )
526: VLW( J ) = VL( JP )
527: 110 CONTINUE
528: IF( ICOMPQ.EQ.1 ) THEN
529: DO 120 J = 2, N
530: JP = IDXP( J )
531: PERM( J ) = IDXQ( IDX( JP )+1 )
532: IF( PERM( J ).LE.NLP1 ) THEN
533: PERM( J ) = PERM( J ) - 1
534: END IF
535: 120 CONTINUE
536: END IF
537: *
538: * The deflated singular values go back into the last N - K slots of
539: * D.
540: *
541: CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
542: *
543: * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
544: * VL(M).
545: *
546: DSIGMA( 1 ) = ZERO
547: HLFTOL = TOL / TWO
548: IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
549: $ DSIGMA( 2 ) = HLFTOL
550: IF( M.GT.N ) THEN
551: Z( 1 ) = DLAPY2( Z1, Z( M ) )
552: IF( Z( 1 ).LE.TOL ) THEN
553: C = ONE
554: S = ZERO
555: Z( 1 ) = TOL
556: ELSE
557: C = Z1 / Z( 1 )
558: S = -Z( M ) / Z( 1 )
559: END IF
560: CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
561: CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
562: ELSE
563: IF( ABS( Z1 ).LE.TOL ) THEN
564: Z( 1 ) = TOL
565: ELSE
566: Z( 1 ) = Z1
567: END IF
568: END IF
569: *
570: * Restore Z, VF, and VL.
571: *
572: CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
573: CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
574: CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
575: *
576: RETURN
577: *
578: * End of DLASD7
579: *
580: 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