Annotation of rpl/lapack/lapack/dlasd8.f, revision 1.3
1.1 bertrand 1: SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
2: $ DSIGMA, WORK, INFO )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * October 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER ICOMPQ, INFO, K, LDDIFR
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
14: $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
15: $ Z( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLASD8 finds the square roots of the roots of the secular equation,
22: * as defined by the values in DSIGMA and Z. It makes the appropriate
23: * calls to DLASD4, and stores, for each element in D, the distance
24: * to its two nearest poles (elements in DSIGMA). It also updates
25: * the arrays VF and VL, the first and last components of all the
26: * right singular vectors of the original bidiagonal matrix.
27: *
28: * DLASD8 is called from DLASD6.
29: *
30: * Arguments
31: * =========
32: *
33: * ICOMPQ (input) INTEGER
34: * Specifies whether singular vectors are to be computed in
35: * factored form in the calling routine:
36: * = 0: Compute singular values only.
37: * = 1: Compute singular vectors in factored form as well.
38: *
39: * K (input) INTEGER
40: * The number of terms in the rational function to be solved
41: * by DLASD4. K >= 1.
42: *
43: * D (output) DOUBLE PRECISION array, dimension ( K )
44: * On output, D contains the updated singular values.
45: *
46: * Z (input/output) DOUBLE PRECISION array, dimension ( K )
47: * On entry, the first K elements of this array contain the
48: * components of the deflation-adjusted updating row vector.
49: * On exit, Z is updated.
50: *
51: * VF (input/output) DOUBLE PRECISION array, dimension ( K )
52: * On entry, VF contains information passed through DBEDE8.
53: * On exit, VF contains the first K components of the first
54: * components of all right singular vectors of the bidiagonal
55: * matrix.
56: *
57: * VL (input/output) DOUBLE PRECISION array, dimension ( K )
58: * On entry, VL contains information passed through DBEDE8.
59: * On exit, VL contains the first K components of the last
60: * components of all right singular vectors of the bidiagonal
61: * matrix.
62: *
63: * DIFL (output) DOUBLE PRECISION array, dimension ( K )
64: * On exit, DIFL(I) = D(I) - DSIGMA(I).
65: *
66: * DIFR (output) DOUBLE PRECISION array,
67: * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
68: * dimension ( K ) if ICOMPQ = 0.
69: * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
70: * defined and will not be referenced.
71: *
72: * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
73: * normalizing factors for the right singular vector matrix.
74: *
75: * LDDIFR (input) INTEGER
76: * The leading dimension of DIFR, must be at least K.
77: *
78: * DSIGMA (input/output) DOUBLE PRECISION array, dimension ( K )
79: * On entry, the first K elements of this array contain the old
80: * roots of the deflated updating problem. These are the poles
81: * of the secular equation.
82: * On exit, the elements of DSIGMA may be very slightly altered
83: * in value.
84: *
85: * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K
86: *
87: * INFO (output) INTEGER
88: * = 0: successful exit.
89: * < 0: if INFO = -i, the i-th argument had an illegal value.
90: * > 0: if INFO = 1, an singular value did not converge
91: *
92: * Further Details
93: * ===============
94: *
95: * Based on contributions by
96: * Ming Gu and Huan Ren, Computer Science Division, University of
97: * California at Berkeley, USA
98: *
99: * =====================================================================
100: *
101: * .. Parameters ..
102: DOUBLE PRECISION ONE
103: PARAMETER ( ONE = 1.0D+0 )
104: * ..
105: * .. Local Scalars ..
106: INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
107: DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
108: * ..
109: * .. External Subroutines ..
110: EXTERNAL DCOPY, DLASCL, DLASD4, DLASET, XERBLA
111: * ..
112: * .. External Functions ..
113: DOUBLE PRECISION DDOT, DLAMC3, DNRM2
114: EXTERNAL DDOT, DLAMC3, DNRM2
115: * ..
116: * .. Intrinsic Functions ..
117: INTRINSIC ABS, SIGN, SQRT
118: * ..
119: * .. Executable Statements ..
120: *
121: * Test the input parameters.
122: *
123: INFO = 0
124: *
125: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
126: INFO = -1
127: ELSE IF( K.LT.1 ) THEN
128: INFO = -2
129: ELSE IF( LDDIFR.LT.K ) THEN
130: INFO = -9
131: END IF
132: IF( INFO.NE.0 ) THEN
133: CALL XERBLA( 'DLASD8', -INFO )
134: RETURN
135: END IF
136: *
137: * Quick return if possible
138: *
139: IF( K.EQ.1 ) THEN
140: D( 1 ) = ABS( Z( 1 ) )
141: DIFL( 1 ) = D( 1 )
142: IF( ICOMPQ.EQ.1 ) THEN
143: DIFL( 2 ) = ONE
144: DIFR( 1, 2 ) = ONE
145: END IF
146: RETURN
147: END IF
148: *
149: * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
150: * be computed with high relative accuracy (barring over/underflow).
151: * This is a problem on machines without a guard digit in
152: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
153: * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
154: * which on any of these machines zeros out the bottommost
155: * bit of DSIGMA(I) if it is 1; this makes the subsequent
156: * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
157: * occurs. On binary machines with a guard digit (almost all
158: * machines) it does not change DSIGMA(I) at all. On hexadecimal
159: * and decimal machines with a guard digit, it slightly
160: * changes the bottommost bits of DSIGMA(I). It does not account
161: * for hexadecimal or decimal machines without guard digits
162: * (we know of none). We use a subroutine call to compute
163: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
164: * this code.
165: *
166: DO 10 I = 1, K
167: DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I )
168: 10 CONTINUE
169: *
170: * Book keeping.
171: *
172: IWK1 = 1
173: IWK2 = IWK1 + K
174: IWK3 = IWK2 + K
175: IWK2I = IWK2 - 1
176: IWK3I = IWK3 - 1
177: *
178: * Normalize Z.
179: *
180: RHO = DNRM2( K, Z, 1 )
181: CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
182: RHO = RHO*RHO
183: *
184: * Initialize WORK(IWK3).
185: *
186: CALL DLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K )
187: *
188: * Compute the updated singular values, the arrays DIFL, DIFR,
189: * and the updated Z.
190: *
191: DO 40 J = 1, K
192: CALL DLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ),
193: $ WORK( IWK2 ), INFO )
194: *
195: * If the root finder fails, the computation is terminated.
196: *
197: IF( INFO.NE.0 ) THEN
198: RETURN
199: END IF
200: WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
201: DIFL( J ) = -WORK( J )
202: DIFR( J, 1 ) = -WORK( J+1 )
203: DO 20 I = 1, J - 1
204: WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
205: $ WORK( IWK2I+I ) / ( DSIGMA( I )-
206: $ DSIGMA( J ) ) / ( DSIGMA( I )+
207: $ DSIGMA( J ) )
208: 20 CONTINUE
209: DO 30 I = J + 1, K
210: WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
211: $ WORK( IWK2I+I ) / ( DSIGMA( I )-
212: $ DSIGMA( J ) ) / ( DSIGMA( I )+
213: $ DSIGMA( J ) )
214: 30 CONTINUE
215: 40 CONTINUE
216: *
217: * Compute updated Z.
218: *
219: DO 50 I = 1, K
220: Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
221: 50 CONTINUE
222: *
223: * Update VF and VL.
224: *
225: DO 80 J = 1, K
226: DIFLJ = DIFL( J )
227: DJ = D( J )
228: DSIGJ = -DSIGMA( J )
229: IF( J.LT.K ) THEN
230: DIFRJ = -DIFR( J, 1 )
231: DSIGJP = -DSIGMA( J+1 )
232: END IF
233: WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
234: DO 60 I = 1, J - 1
235: WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
236: $ / ( DSIGMA( I )+DJ )
237: 60 CONTINUE
238: DO 70 I = J + 1, K
239: WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
240: $ / ( DSIGMA( I )+DJ )
241: 70 CONTINUE
242: TEMP = DNRM2( K, WORK, 1 )
243: WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP
244: WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP
245: IF( ICOMPQ.EQ.1 ) THEN
246: DIFR( J, 2 ) = TEMP
247: END IF
248: 80 CONTINUE
249: *
250: CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 )
251: CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 )
252: *
253: RETURN
254: *
255: * End of DLASD8
256: *
257: END
258:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd8.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack