zheevc3.cxx File Reference

#include <stdio.h>
#include <cmath>
#include <complex>
#include "zheevc3.h"

Go to the source code of this file.

Macros
#define	M_SQRT3   1.73205080756887729352744634151
#define	SQR(x)   ((x)*(x))
#define	SQR_ABS(x)   (SQR(real(x)) + SQR(imag(x)))

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

Macro Definition Documentation

M_SQRT3

#define M_SQRT3   1.73205080756887729352744634151

Definition at line 25 of file zheevc3.cxx.

SQR

#define SQR

(

x

)

((x)*(x))

Definition at line 28 of file zheevc3.cxx.

SQR_ABS

#define SQR_ABS

(

x

)

(SQR(real(x)) + SQR(imag(x)))

Definition at line 29 of file zheevc3.cxx.

Function Documentation

zheevc3()

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

Definition at line 33 of file zheevc3.cxx.

{
  double m, c1, c0;
 
  // Determine coefficients of characteristic poynomial. We write
  //       | a   d   f  |
  //  A =  | d*  b   e  |
  //       | f*  e*  c  |
  std::complex<double> de = A[0][1] * A[1][2];                            // d * e
  double dd = SQR_ABS(A[0][1]);                                  // d * conj(d)
  double ee = SQR_ABS(A[1][2]);                                  // e * conj(e)
  double ff = SQR_ABS(A[0][2]);                                  // f * conj(f)
  m  = real(A[0][0]) + real(A[1][1]) + real(A[2][2]);
  c1 = (real(A[0][0])*real(A[1][1])  // a*b + a*c + b*c - d*conj(d) - e*conj(e) - f*conj(f)
          + real(A[0][0])*real(A[2][2])
          + real(A[1][1])*real(A[2][2]))
          - (dd + ee + ff);
  c0 = real(A[2][2])*dd + real(A[0][0])*ee + real(A[1][1])*ff
            - real(A[0][0])*real(A[1][1])*real(A[2][2])
            - 2.0 * (real(A[0][2])*real(de) + imag(A[0][2])*imag(de));
                             // c*d*conj(d) + a*e*conj(e) + b*f*conj(f) - a*b*c - 2*Re(conj(f)*d*e)
 
  double p, sqrt_p, q, c, s, phi;
  p = SQR(m) - 3.0*c1;
  q = m*(p - (3.0/2.0)*c1) - (27.0/2.0)*c0;
  sqrt_p = sqrt(fabs(p));
 
  phi = 27.0 * ( 0.25*SQR(c1)*(p - c1) + c0*(q + 27.0/4.0*c0));
  phi = (1.0/3.0) * atan2(sqrt(fabs(phi)), q);
 
  c = sqrt_p*cos(phi);
  s = (1.0/M_SQRT3)*sqrt_p*sin(phi);
 
  w[1]  = (1.0/3.0)*(m - c);
  w[2]  = w[1] + s;
  w[0]  = w[1] + c;
  w[1] -= s;
 
  return 0;
}

References M_SQRT3, SQR, and SQR_ABS.

Referenced by zheevh3().