1*> \brief \b SLASD5 computes the square root of the i-th eigenvalue of a positive symmetric rank-one modification of a 2-by-2 diagonal matrix. 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 SLASD5 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd5.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd5.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd5.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
22*
23*       .. Scalar Arguments ..
24*       INTEGER            I
25*       REAL               DSIGMA, RHO
26*       ..
27*       .. Array Arguments ..
28*       REAL               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 REAL 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 REAL array, dimension (2)
69*>         The components of the updating vector.
70*> \endverbatim
71*>
72*> \param[out] DELTA
73*> \verbatim
74*>          DELTA is REAL 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 REAL
83*>         The scalar in the symmetric updating formula.
84*> \endverbatim
85*>
86*> \param[out] DSIGMA
87*> \verbatim
88*>          DSIGMA is REAL
89*>         The computed sigma_I, the I-th updated eigenvalue.
90*> \endverbatim
91*>
92*> \param[out] WORK
93*> \verbatim
94*>          WORK is REAL 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*> \ingroup OTHERauxiliary
107*
108*> \par Contributors:
109*  ==================
110*>
111*>     Ren-Cang Li, Computer Science Division, University of California
112*>     at Berkeley, USA
113*>
114*  =====================================================================
115      SUBROUTINE SLASD5( I, D, Z, DELTA, RHO, DSIGMA, WORK )
116*
117*  -- LAPACK auxiliary routine --
118*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
119*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120*
121*     .. Scalar Arguments ..
122      INTEGER            I
123      REAL               DSIGMA, RHO
124*     ..
125*     .. Array Arguments ..
126      REAL               D( 2 ), DELTA( 2 ), WORK( 2 ), Z( 2 )
127*     ..
128*
129*  =====================================================================
130*
131*     .. Parameters ..
132      REAL               ZERO, ONE, TWO, THREE, FOUR
133      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0,
134     $                   THREE = 3.0E+0, FOUR = 4.0E+0 )
135*     ..
136*     .. Local Scalars ..
137      REAL               B, C, DEL, DELSQ, TAU, W
138*     ..
139*     .. Intrinsic Functions ..
140      INTRINSIC          ABS, SQRT
141*     ..
142*     .. Executable Statements ..
143*
144      DEL = D( 2 ) - D( 1 )
145      DELSQ = DEL*( D( 2 )+D( 1 ) )
146      IF( I.EQ.1 ) THEN
147         W = ONE + FOUR*RHO*( Z( 2 )*Z( 2 ) / ( D( 1 )+THREE*D( 2 ) )-
148     $       Z( 1 )*Z( 1 ) / ( THREE*D( 1 )+D( 2 ) ) ) / DEL
149         IF( W.GT.ZERO ) THEN
150            B = DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
151            C = RHO*Z( 1 )*Z( 1 )*DELSQ
152*
153*           B > ZERO, always
154*
155*           The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 )
156*
157            TAU = TWO*C / ( B+SQRT( ABS( B*B-FOUR*C ) ) )
158*
159*           The following TAU is DSIGMA - D( 1 )
160*
161            TAU = TAU / ( D( 1 )+SQRT( D( 1 )*D( 1 )+TAU ) )
162            DSIGMA = D( 1 ) + TAU
163            DELTA( 1 ) = -TAU
164            DELTA( 2 ) = DEL - TAU
165            WORK( 1 ) = TWO*D( 1 ) + TAU
166            WORK( 2 ) = ( D( 1 )+TAU ) + D( 2 )
167*           DELTA( 1 ) = -Z( 1 ) / TAU
168*           DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )
169         ELSE
170            B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
171            C = RHO*Z( 2 )*Z( 2 )*DELSQ
172*
173*           The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
174*
175            IF( B.GT.ZERO ) THEN
176               TAU = -TWO*C / ( B+SQRT( B*B+FOUR*C ) )
177            ELSE
178               TAU = ( B-SQRT( B*B+FOUR*C ) ) / TWO
179            END IF
180*
181*           The following TAU is DSIGMA - D( 2 )
182*
183            TAU = TAU / ( D( 2 )+SQRT( ABS( D( 2 )*D( 2 )+TAU ) ) )
184            DSIGMA = D( 2 ) + TAU
185            DELTA( 1 ) = -( DEL+TAU )
186            DELTA( 2 ) = -TAU
187            WORK( 1 ) = D( 1 ) + TAU + D( 2 )
188            WORK( 2 ) = TWO*D( 2 ) + TAU
189*           DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
190*           DELTA( 2 ) = -Z( 2 ) / TAU
191         END IF
192*        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
193*        DELTA( 1 ) = DELTA( 1 ) / TEMP
194*        DELTA( 2 ) = DELTA( 2 ) / TEMP
195      ELSE
196*
197*        Now I=2
198*
199         B = -DELSQ + RHO*( Z( 1 )*Z( 1 )+Z( 2 )*Z( 2 ) )
200         C = RHO*Z( 2 )*Z( 2 )*DELSQ
201*
202*        The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )
203*
204         IF( B.GT.ZERO ) THEN
205            TAU = ( B+SQRT( B*B+FOUR*C ) ) / TWO
206         ELSE
207            TAU = TWO*C / ( -B+SQRT( B*B+FOUR*C ) )
208         END IF
209*
210*        The following TAU is DSIGMA - D( 2 )
211*
212         TAU = TAU / ( D( 2 )+SQRT( D( 2 )*D( 2 )+TAU ) )
213         DSIGMA = D( 2 ) + TAU
214         DELTA( 1 ) = -( DEL+TAU )
215         DELTA( 2 ) = -TAU
216         WORK( 1 ) = D( 1 ) + TAU + D( 2 )
217         WORK( 2 ) = TWO*D( 2 ) + TAU
218*        DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )
219*        DELTA( 2 ) = -Z( 2 ) / TAU
220*        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )
221*        DELTA( 1 ) = DELTA( 1 ) / TEMP
222*        DELTA( 2 ) = DELTA( 2 ) / TEMP
223      END IF
224      RETURN
225*
226*     End of SLASD5
227*
228      END
229