1: *> \brief \b DLASD5
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASD5 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd5.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd5.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd5.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER I
25: * DOUBLE PRECISION DSIGMA, RHO
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> This subroutine computes the square root of the I-th eigenvalue
38: *> of a positive symmetric rank-one modification of a 2-by-2 diagonal
39: *> matrix
40: *>
41: *> diag( D ) * diag( D ) + RHO * Z * transpose(Z) .
42: *>
43: *> The diagonal entries in the array D are assumed to satisfy
44: *>
45: *> 0 <= D(i) < D(j) for i < j .
46: *>
47: *> We also assume RHO > 0 and that the Euclidean norm of the vector
48: *> Z is one.
49: *> \endverbatim
50: *
51: * Arguments:
52: * ==========
53: *
54: *> \param[in] I
55: *> \verbatim
56: *> I is INTEGER
57: *> The index of the eigenvalue to be computed. I = 1 or I = 2.
58: *> \endverbatim
59: *>
60: *> \param[in] D
61: *> \verbatim
62: *> D is DOUBLE PRECISION array, dimension ( 2 )
63: *> The original eigenvalues. We assume 0 <= D(1) < D(2).
64: *> \endverbatim
65: *>
66: *> \param[in] Z
67: *> \verbatim
68: *> Z is DOUBLE PRECISION array, dimension ( 2 )
69: *> The components of the updating vector.
70: *> \endverbatim
71: *>
72: *> \param[out] DELTA
73: *> \verbatim
74: *> DELTA is DOUBLE PRECISION array, dimension ( 2 )
75: *> Contains (D(j) - sigma_I) in its j-th component.
76: *> The vector DELTA contains the information necessary
77: *> to construct the eigenvectors.
78: *> \endverbatim
79: *>
80: *> \param[in] RHO
81: *> \verbatim
82: *> RHO is DOUBLE PRECISION
83: *> The scalar in the symmetric updating formula.
84: *> \endverbatim
85: *>
86: *> \param[out] DSIGMA
87: *> \verbatim
88: *> DSIGMA is DOUBLE PRECISION
89: *> The computed sigma_I, the I-th updated eigenvalue.
90: *> \endverbatim
91: *>
92: *> \param[out] WORK
93: *> \verbatim
94: *> WORK is DOUBLE PRECISION array, dimension ( 2 )
95: *> WORK contains (D(j) + sigma_I) in its j-th component.
96: *> \endverbatim
97: *
98: * Authors:
99: * ========
100: *
101: *> \author Univ. of Tennessee
102: *> \author Univ. of California Berkeley
103: *> \author Univ. of Colorado Denver
104: *> \author NAG Ltd.
105: *
106: *> \date November 2011
107: *
108: *> \ingroup auxOTHERauxiliary
109: *
110: *> \par Contributors:
111: * ==================
112: *>
113: *> Ren-Cang Li, Computer Science Division, University of California
114: *> at Berkeley, USA
115: *>
116: * =====================================================================
117: SUBROUTINE DLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
118: *
119: * -- LAPACK auxiliary routine (version 3.4.0) --
120: * -- LAPACK is a software package provided by Univ. of Tennessee, --
121: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
122: * November 2011
123: *
124: * .. Scalar Arguments ..
125: INTEGER I
126: DOUBLE PRECISION DSIGMA, RHO
127: * ..
128: * .. Array Arguments ..
129: DOUBLE PRECISION D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
130: * ..
131: *
132: * =====================================================================
133: *
134: * .. Parameters ..
135: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR
136: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
137: $ THREE = 3.0D+0, FOUR = 4.0D+0 )
138: * ..
139: * .. Local Scalars ..
140: DOUBLE PRECISION B, C, DEL, DELSQ, TAU, W
141: * ..
142: * .. Intrinsic Functions ..
143: INTRINSIC ABS, SQRT
144: * ..
145: * .. Executable Statements ..
146: *
147: DEL = D( 2 ) - D( 1 )
148: DELSQ = DEL*( D( 2 )+D( 1 ) )
149: IF( I.EQ.1 ) THEN
150: W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )-
151: $ Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL
152: IF( W.GT.ZERO ) THEN
153: B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
154: C = RHO*Z( 1 )*Z( 1 )*DELSQ
155: *
156: * B > ZERO, always
157: *
158: * The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 )
159: *
160: TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) )
161: *
162: * The following TAU is DSIGMA - D( 1 )
163: *
164: TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) )
165: DSIGMA = D( 1 ) + TAU
166: DELTA( 1 ) = -TAU
167: DELTA( 2 ) = DEL - TAU
168: WORK( 1 ) = TWO*D( 1 ) + TAU
169: WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 )
170: * DELTA( 1 ) = -Z( 1 ) / TAU
171: * DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )
172: ELSE
173: B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
174: C = RHO*Z( 2 )*Z( 2 )*DELSQ
175: *
176: * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
177: *
178: IF( B.GT.ZERO ) THEN
179: TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) )
180: ELSE
181: TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO
182: END IF
183: *
184: * The following TAU is DSIGMA - D( 2 )
185: *
186: TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) )
187: DSIGMA = D( 2 ) + TAU
188: DELTA( 1 ) = -( DEL+TAU )
189: DELTA( 2 ) = -TAU
190: WORK( 1 ) = D( 1 ) + TAU + D( 2 )
191: WORK( 2 ) = TWO*D( 2 ) + TAU
192: * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
193: * DELTA( 2 ) = -Z( 2 ) / TAU
194: END IF
195: * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
196: * DELTA( 1 ) = DELTA( 1 ) / TEMP
197: * DELTA( 2 ) = DELTA( 2 ) / TEMP
198: ELSE
199: *
200: * Now I=2
201: *
202: B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
203: C = RHO*Z( 2 )*Z( 2 )*DELSQ
204: *
205: * The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
206: *
207: IF( B.GT.ZERO ) THEN
208: TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO
209: ELSE
210: TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) )
211: END IF
212: *
213: * The following TAU is DSIGMA - D( 2 )
214: *
215: TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) )
216: DSIGMA = D( 2 ) + TAU
217: DELTA( 1 ) = -( DEL+TAU )
218: DELTA( 2 ) = -TAU
219: WORK( 1 ) = D( 1 ) + TAU + D( 2 )
220: WORK( 2 ) = TWO*D( 2 ) + TAU
221: * DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
222: * DELTA( 2 ) = -Z( 2 ) / TAU
223: * TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
224: * DELTA( 1 ) = DELTA( 1 ) / TEMP
225: * DELTA( 2 ) = DELTA( 2 ) / TEMP
226: END IF
227: RETURN
228: *
229: * End of DLASD5
230: *
231: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd5.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack