Annotation of rpl/lapack/lapack/dlaev2.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLAEV2
! 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 November 2011
! 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: * =====================================================================
1.1 bertrand 121: SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 )
122: *
1.8 ! bertrand 123: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 124: * -- LAPACK is a software package provided by Univ. of Tennessee, --
125: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 126: * November 2011
1.1 bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaev2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack