← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dlaev2.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linalg/dlaev2.f
 
1
*> \brief \b DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*> \htmlonly
 
9
*> Download DLAEV2 + dependencies 
 
10
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaev2.f"> 
 
11
*> [TGZ]</a> 
 
12
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaev2.f"> 
 
13
*> [ZIP]</a> 
 
14
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaev2.f"> 
 
15
*> [TXT]</a>
 
16
*> \endhtmlonly 
 
17
*
 
18
*  Definition:
 
19
*  ===========
 
20
*
 
21
*       SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
 
22
* 
 
23
*       .. Scalar Arguments ..
 
24
*       DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
 
25
*       ..
 
26
*  
 
27
*
 
28
*> \par Purpose:
 
29
*  =============
 
30
*>
 
31
*> \verbatim
 
32
*>
 
33
*> DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
 
34
*>    [  A   B  ]
 
35
*>    [  B   C  ].
 
36
*> On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
 
37
*> eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
 
38
*> eigenvector for RT1, giving the decomposition
 
39
*>
 
40
*>    [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
 
41
*>    [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
 
42
*> \endverbatim
 
43
*
 
44
*  Arguments:
 
45
*  ==========
 
46
*
 
47
*> \param[in] A
 
48
*> \verbatim
 
49
*>          A is DOUBLE PRECISION
 
50
*>          The (1,1) element of the 2-by-2 matrix.
 
51
*> \endverbatim
 
52
*>
 
53
*> \param[in] B
 
54
*> \verbatim
 
55
*>          B is DOUBLE PRECISION
 
56
*>          The (1,2) element and the conjugate of the (2,1) element of
 
57
*>          the 2-by-2 matrix.
 
58
*> \endverbatim
 
59
*>
 
60
*> \param[in] C
 
61
*> \verbatim
 
62
*>          C is DOUBLE PRECISION
 
63
*>          The (2,2) element of the 2-by-2 matrix.
 
64
*> \endverbatim
 
65
*>
 
66
*> \param[out] RT1
 
67
*> \verbatim
 
68
*>          RT1 is DOUBLE PRECISION
 
69
*>          The eigenvalue of larger absolute value.
 
70
*> \endverbatim
 
71
*>
 
72
*> \param[out] RT2
 
73
*> \verbatim
 
74
*>          RT2 is DOUBLE PRECISION
 
75
*>          The eigenvalue of smaller absolute value.
 
76
*> \endverbatim
 
77
*>
 
78
*> \param[out] CS1
 
79
*> \verbatim
 
80
*>          CS1 is DOUBLE PRECISION
 
81
*> \endverbatim
 
82
*>
 
83
*> \param[out] SN1
 
84
*> \verbatim
 
85
*>          SN1 is DOUBLE PRECISION
 
86
*>          The vector (CS1, SN1) is a unit right eigenvector for RT1.
 
87
*> \endverbatim
 
88
*
 
89
*  Authors:
 
90
*  ========
 
91
*
 
92
*> \author Univ. of Tennessee 
 
93
*> \author Univ. of California Berkeley 
 
94
*> \author Univ. of Colorado Denver 
 
95
*> \author NAG Ltd. 
 
96
*
 
97
*> \date September 2012
 
98
*
 
99
*> \ingroup auxOTHERauxiliary
 
100
*
 
101
*> \par Further Details:
 
102
*  =====================
 
103
*>
 
104
*> \verbatim
 
105
*>
 
106
*>  RT1 is accurate to a few ulps barring over/underflow.
 
107
*>
 
108
*>  RT2 may be inaccurate if there is massive cancellation in the
 
109
*>  determinant A*C-B*B; higher precision or correctly rounded or
 
110
*>  correctly truncated arithmetic would be needed to compute RT2
 
111
*>  accurately in all cases.
 
112
*>
 
113
*>  CS1 and SN1 are accurate to a few ulps barring over/underflow.
 
114
*>
 
115
*>  Overflow is possible only if RT1 is within a factor of 5 of overflow.
 
116
*>  Underflow is harmless if the input data is 0 or exceeds
 
117
*>     underflow_threshold / macheps.
 
118
*> \endverbatim
 
119
*>
 
120
*  =====================================================================
 
121
      SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
 
122
*
 
123
*  -- LAPACK auxiliary routine (version 3.4.2) --
 
124
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
125
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
126
*     September 2012
 
127
*
 
128
*     .. Scalar Arguments ..
 
129
      DOUBLE PRECISION   A, B, C, CS1, RT1, RT2, SN1
 
130
*     ..
 
131
*
 
132
* =====================================================================
 
133
*
 
134
*     .. Parameters ..
 
135
      DOUBLE PRECISION   ONE
 
136
      PARAMETER          ( ONE = 1.0D0 )
 
137
      DOUBLE PRECISION   TWO
 
138
      PARAMETER          ( TWO = 2.0D0 )
 
139
      DOUBLE PRECISION   ZERO
 
140
      PARAMETER          ( ZERO = 0.0D0 )
 
141
      DOUBLE PRECISION   HALF
 
142
      PARAMETER          ( HALF = 0.5D0 )
 
143
*     ..
 
144
*     .. Local Scalars ..
 
145
      INTEGER            SGN1, SGN2
 
146
      DOUBLE PRECISION   AB, ACMN, ACMX, ACS, ADF, CS, CT, DF, RT, SM,
 
147
     $                   TB, TN
 
148
*     ..
 
149
*     .. Intrinsic Functions ..
 
150
      INTRINSIC          ABS, SQRT
 
151
*     ..
 
152
*     .. Executable Statements ..
 
153
*
 
154
*     Compute the eigenvalues
 
155
*
 
156
      SM = A + C
 
157
      DF = A - C
 
158
      ADF = ABS( DF )
 
159
      TB = B + B
 
160
      AB = ABS( TB )
 
161
      IF( ABS( A ).GT.ABS( C ) ) THEN
 
162
         ACMX = A
 
163
         ACMN = C
 
164
      ELSE
 
165
         ACMX = C
 
166
         ACMN = A
 
167
      END IF
 
168
      IF( ADF.GT.AB ) THEN
 
169
         RT = ADF*SQRT( ONE+( AB / ADF )**2 )
 
170
      ELSE IF( ADF.LT.AB ) THEN
 
171
         RT = AB*SQRT( ONE+( ADF / AB )**2 )
 
172
      ELSE
 
173
*
 
174
*        Includes case AB=ADF=0
 
175
*
 
176
         RT = AB*SQRT( TWO )
 
177
      END IF
 
178
      IF( SM.LT.ZERO ) THEN
 
179
         RT1 = HALF*( SM-RT )
 
180
         SGN1 = -1
 
181
*
 
182
*        Order of execution important.
 
183
*        To get fully accurate smaller eigenvalue,
 
184
*        next line needs to be executed in higher precision.
 
185
*
 
186
         RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
 
187
      ELSE IF( SM.GT.ZERO ) THEN
 
188
         RT1 = HALF*( SM+RT )
 
189
         SGN1 = 1
 
190
*
 
191
*        Order of execution important.
 
192
*        To get fully accurate smaller eigenvalue,
 
193
*        next line needs to be executed in higher precision.
 
194
*
 
195
         RT2 = ( ACMX / RT1 )*ACMN - ( B / RT1 )*B
 
196
      ELSE
 
197
*
 
198
*        Includes case RT1 = RT2 = 0
 
199
*
 
200
         RT1 = HALF*RT
 
201
         RT2 = -HALF*RT
 
202
         SGN1 = 1
 
203
      END IF
 
204
*
 
205
*     Compute the eigenvector
 
206
*
 
207
      IF( DF.GE.ZERO ) THEN
 
208
         CS = DF + RT
 
209
         SGN2 = 1
 
210
      ELSE
 
211
         CS = DF - RT
 
212
         SGN2 = -1
 
213
      END IF
 
214
      ACS = ABS( CS )
 
215
      IF( ACS.GT.AB ) THEN
 
216
         CT = -TB / CS
 
217
         SN1 = ONE / SQRT( ONE+CT*CT )
 
218
         CS1 = CT*SN1
 
219
      ELSE
 
220
         IF( AB.EQ.ZERO ) THEN
 
221
            CS1 = ONE
 
222
            SN1 = ZERO
 
223
         ELSE
 
224
            TN = -CS / TB
 
225
            CS1 = ONE / SQRT( ONE+TN*TN )
 
226
            SN1 = TN*CS1
 
227
         END IF
 
228
      END IF
 
229
      IF( SGN1.EQ.SGN2 ) THEN
 
230
         TN = CS1
 
231
         CS1 = -SN1
 
232
         SN1 = TN
 
233
      END IF
 
234
      RETURN
 
235
*
 
236
*     End of DLAEV2
 
237
*
 
238
      END