zheevc3.h File Reference

#include <complex>

Go to the source code of this file.

Functions

int zheevc3 (std::complex< double > A[3][3], double w[3])

Function Documentation

int zheevc3 ( std::complex< double >  A[3][3],
double  w[3] 
)

Definition at line 33 of file zheevc3.cxx.

References Munits::m, M_SQRT3, SQR, and SQR_ABS.

Referenced by zheevh3().

00040              :
00041 //   A: The hermitian input matrix
00042 //   w: Storage buffer for eigenvalues
00043 // ----------------------------------------------------------------------------
00044 // Return value:
00045 //   0: Success
00046 //  -1: Error
00047 // ----------------------------------------------------------------------------
00048 {
00049   double m, c1, c0;
00050   
00051   // Determine coefficients of characteristic poynomial. We write
00052   //       | a   d   f  |
00053   //  A =  | d*  b   e  |
00054   //       | f*  e*  c  |
00055   std::complex<double> de = A[0][1] * A[1][2];                            // d * e
00056   double dd = SQR_ABS(A[0][1]);                                  // d * conj(d)
00057   double ee = SQR_ABS(A[1][2]);                                  // e * conj(e)
00058   double ff = SQR_ABS(A[0][2]);                                  // f * conj(f)
00059   m  = real(A[0][0]) + real(A[1][1]) + real(A[2][2]);
00060   c1 = (real(A[0][0])*real(A[1][1])  // a*b + a*c + b*c - d*conj(d) - e*conj(e) - f*conj(f)
00061           + real(A[0][0])*real(A[2][2])
00062           + real(A[1][1])*real(A[2][2]))
00063           - (dd + ee + ff);
00064   c0 = real(A[2][2])*dd + real(A[0][0])*ee + real(A[1][1])*ff
00065             - real(A[0][0])*real(A[1][1])*real(A[2][2])
00066             - 2.0 * (real(A[0][2])*real(de) + imag(A[0][2])*imag(de));
00067                              // c*d*conj(d) + a*e*conj(e) + b*f*conj(f) - a*b*c - 2*Re(conj(f)*d*e)
00068 
00069   double p, sqrt_p, q, c, s, phi;
00070   p = SQR(m) - 3.0*c1;
00071   q = m*(p - (3.0/2.0)*c1) - (27.0/2.0)*c0;
00072   sqrt_p = sqrt(fabs(p));
00073 
00074   phi = 27.0 * ( 0.25*SQR(c1)*(p - c1) + c0*(q + 27.0/4.0*c0));
00075   phi = (1.0/3.0) * atan2(sqrt(fabs(phi)), q);
00076   
00077   c = sqrt_p*cos(phi);
00078   s = (1.0/M_SQRT3)*sqrt_p*sin(phi);
00079 
00080   w[1]  = (1.0/3.0)*(m - c);
00081   w[2]  = w[1] + s;
00082   w[0]  = w[1] + c;
00083   w[1] -= s;
00084 
00085   return 0;
00086 }


Generated on 22 Nov 2017 for loon by  doxygen 1.6.1