1*> \brief \b DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLANV2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlanv2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlanv2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlanv2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
22*
23*       .. Scalar Arguments ..
24*       DOUBLE PRECISION   A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
25*       ..
26*
27*
28*> \par Purpose:
29*  =============
30*>
31*> \verbatim
32*>
33*> DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
34*> matrix in standard form:
35*>
36*>      [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
37*>      [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
38*>
39*> where either
40*> 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
41*> 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
42*> conjugate eigenvalues.
43*> \endverbatim
44*
45*  Arguments:
46*  ==========
47*
48*> \param[in,out] A
49*> \verbatim
50*>          A is DOUBLE PRECISION
51*> \endverbatim
52*>
53*> \param[in,out] B
54*> \verbatim
55*>          B is DOUBLE PRECISION
56*> \endverbatim
57*>
58*> \param[in,out] C
59*> \verbatim
60*>          C is DOUBLE PRECISION
61*> \endverbatim
62*>
63*> \param[in,out] D
64*> \verbatim
65*>          D is DOUBLE PRECISION
66*>          On entry, the elements of the input matrix.
67*>          On exit, they are overwritten by the elements of the
68*>          standardised Schur form.
69*> \endverbatim
70*>
71*> \param[out] RT1R
72*> \verbatim
73*>          RT1R is DOUBLE PRECISION
74*> \endverbatim
75*>
76*> \param[out] RT1I
77*> \verbatim
78*>          RT1I is DOUBLE PRECISION
79*> \endverbatim
80*>
81*> \param[out] RT2R
82*> \verbatim
83*>          RT2R is DOUBLE PRECISION
84*> \endverbatim
85*>
86*> \param[out] RT2I
87*> \verbatim
88*>          RT2I is DOUBLE PRECISION
89*>          The real and imaginary parts of the eigenvalues. If the
90*>          eigenvalues are a complex conjugate pair, RT1I > 0.
91*> \endverbatim
92*>
93*> \param[out] CS
94*> \verbatim
95*>          CS is DOUBLE PRECISION
96*> \endverbatim
97*>
98*> \param[out] SN
99*> \verbatim
100*>          SN is DOUBLE PRECISION
101*>          Parameters of the rotation matrix.
102*> \endverbatim
103*
104*  Authors:
105*  ========
106*
107*> \author Univ. of Tennessee
108*> \author Univ. of California Berkeley
109*> \author Univ. of Colorado Denver
110*> \author NAG Ltd.
111*
112*> \date December 2016
113*
114*> \ingroup doubleOTHERauxiliary
115*
116*> \par Further Details:
117*  =====================
118*>
119*> \verbatim
120*>
121*>  Modified by V. Sima, Research Institute for Informatics, Bucharest,
122*>  Romania, to reduce the risk of cancellation errors,
123*>  when computing real eigenvalues, and to ensure, if possible, that
124*>  abs(RT1R) >= abs(RT2R).
125*> \endverbatim
126*>
127*  =====================================================================
128      SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
129*
130*  -- LAPACK auxiliary routine (version 3.7.0) --
131*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
132*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
133*     December 2016
134*
135*     .. Scalar Arguments ..
136      DOUBLE PRECISION   A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
137*     ..
138*
139*  =====================================================================
140*
141*     .. Parameters ..
142      DOUBLE PRECISION   ZERO, HALF, ONE, TWO
143      PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0,
144     $                     TWO = 2.0D0 )
145      DOUBLE PRECISION   MULTPL
146      PARAMETER          ( MULTPL = 4.0D+0 )
147*     ..
148*     .. Local Scalars ..
149      DOUBLE PRECISION   AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB,
150     $                   SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z, SAFMIN,
151     $                   SAFMN2, SAFMX2
152      INTEGER            COUNT
153*     ..
154*     .. External Functions ..
155      DOUBLE PRECISION   DLAMCH, DLAPY2
156      EXTERNAL           DLAMCH, DLAPY2
157*     ..
158*     .. Intrinsic Functions ..
159      INTRINSIC          ABS, MAX, MIN, SIGN, SQRT
160*     ..
161*     .. Executable Statements ..
162*
163      SAFMIN = DLAMCH( 'S' )
164      EPS = DLAMCH( 'P' )
165      SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) /
166     $            LOG( DLAMCH( 'B' ) ) / TWO )
167      SAFMX2 = ONE / SAFMN2
168      IF( C.EQ.ZERO ) THEN
169         CS = ONE
170         SN = ZERO
171*
172      ELSE IF( B.EQ.ZERO ) THEN
173*
174*        Swap rows and columns
175*
176         CS = ZERO
177         SN = ONE
178         TEMP = D
179         D = A
180         A = TEMP
181         B = -C
182         C = ZERO
183*
184      ELSE IF( ( A-D ).EQ.ZERO .AND. SIGN( ONE, B ).NE.SIGN( ONE, C ) )
185     $          THEN
186         CS = ONE
187         SN = ZERO
188*
189      ELSE
190*
191         TEMP = A - D
192         P = HALF*TEMP
193         BCMAX = MAX( ABS( B ), ABS( C ) )
194         BCMIS = MIN( ABS( B ), ABS( C ) )*SIGN( ONE, B )*SIGN( ONE, C )
195         SCALE = MAX( ABS( P ), BCMAX )
196         Z = ( P / SCALE )*P + ( BCMAX / SCALE )*BCMIS
197*
198*        If Z is of the order of the machine accuracy, postpone the
199*        decision on the nature of eigenvalues
200*
201         IF( Z.GE.MULTPL*EPS ) THEN
202*
203*           Real eigenvalues. Compute A and D.
204*
205            Z = P + SIGN( SQRT( SCALE )*SQRT( Z ), P )
206            A = D + Z
207            D = D - ( BCMAX / Z )*BCMIS
208*
209*           Compute B and the rotation matrix
210*
211            TAU = DLAPY2( C, Z )
212            CS = Z / TAU
213            SN = C / TAU
214            B = B - C
215            C = ZERO
216*
217         ELSE
218*
219*           Complex eigenvalues, or real (almost) equal eigenvalues.
220*           Make diagonal elements equal.
221*
222            COUNT = 0
223            SIGMA = B + C
224   10       CONTINUE
225            COUNT = COUNT + 1
226            SCALE = MAX( ABS(TEMP), ABS(SIGMA) )
227            IF( SCALE.GE.SAFMX2 ) THEN
228               SIGMA = SIGMA * SAFMN2
229               TEMP = TEMP * SAFMN2
230               IF (COUNT .LE. 20)
231     $            GOTO 10
232            END IF
233            IF( SCALE.LE.SAFMN2 ) THEN
234               SIGMA = SIGMA * SAFMX2
235               TEMP = TEMP * SAFMX2
236               IF (COUNT .LE. 20)
237     $            GOTO 10
238            END IF
239            P = HALF*TEMP
240            TAU = DLAPY2( SIGMA, TEMP )
241            CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
242            SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )
243*
244*           Compute [ AA  BB ] = [ A  B ] [ CS -SN ]
245*                   [ CC  DD ]   [ C  D ] [ SN  CS ]
246*
247            AA = A*CS + B*SN
248            BB = -A*SN + B*CS
249            CC = C*CS + D*SN
250            DD = -C*SN + D*CS
251*
252*           Compute [ A  B ] = [ CS  SN ] [ AA  BB ]
253*                   [ C  D ]   [-SN  CS ] [ CC  DD ]
254*
255            A = AA*CS + CC*SN
256            B = BB*CS + DD*SN
257            C = -AA*SN + CC*CS
258            D = -BB*SN + DD*CS
259*
260            TEMP = HALF*( A+D )
261            A = TEMP
262            D = TEMP
263*
264            IF( C.NE.ZERO ) THEN
265               IF( B.NE.ZERO ) THEN
266                  IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
267*
268*                    Real eigenvalues: reduce to upper triangular form
269*
270                     SAB = SQRT( ABS( B ) )
271                     SAC = SQRT( ABS( C ) )
272                     P = SIGN( SAB*SAC, C )
273                     TAU = ONE / SQRT( ABS( B+C ) )
274                     A = TEMP + P
275                     D = TEMP - P
276                     B = B - C
277                     C = ZERO
278                     CS1 = SAB*TAU
279                     SN1 = SAC*TAU
280                     TEMP = CS*CS1 - SN*SN1
281                     SN = CS*SN1 + SN*CS1
282                     CS = TEMP
283                  END IF
284               ELSE
285                  B = -C
286                  C = ZERO
287                  TEMP = CS
288                  CS = -SN
289                  SN = TEMP
290               END IF
291            END IF
292         END IF
293*
294      END IF
295*
296*     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
297*
298      RT1R = A
299      RT2R = D
300      IF( C.EQ.ZERO ) THEN
301         RT1I = ZERO
302         RT2I = ZERO
303      ELSE
304         RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
305         RT2I = -RT1I
306      END IF
307      RETURN
308*
309*     End of DLANV2
310*
311      END
312