zheevh3.cxx File Reference

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

Go to the source code of this file.

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

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

Macro Definition Documentation

SQR

#define SQR

(

x

)

((x)*(x))

Definition at line 27 of file zheevh3.cxx.

SQR_ABS

#define SQR_ABS

(

x

)

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

Definition at line 28 of file zheevh3.cxx.

Function Documentation

zheevh3()

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

Definition at line 31 of file zheevh3.cxx.

{
#ifndef EVALS_ONLY
  double norm;          // Squared norm or inverse norm of current eigenvector
//  double n0, n1;        // Norm of first and second columns of A
  double error;         // Estimated maximum roundoff error
  double t, u;          // Intermediate storage
  int j;                // Loop counter
#endif
 
  // Calculate eigenvalues
  zheevc3(A, w);
 
#ifndef EVALS_ONLY
//  n0 = SQR(real(A[0][0])) + SQR_ABS(A[0][1]) + SQR_ABS(A[0][2]);
//  n1 = SQR_ABS(A[0][1]) + SQR(real(A[1][1])) + SQR_ABS(A[1][2]);
 
  t = fabs(w[0]);
  if ((u=fabs(w[1])) > t)
    t = u;
  if ((u=fabs(w[2])) > t)
    t = u;
  if (t < 1.0)
    u = t;
  else
    u = SQR(t);
  error = 256.0 * DBL_EPSILON * SQR(u);
//  error = 256.0 * DBL_EPSILON * (n0 + u) * (n1 + u);
 
  Q[0][1] = A[0][1]*A[1][2] - A[0][2]*real(A[1][1]);
  Q[1][1] = A[0][2]*conj(A[0][1]) - A[1][2]*real(A[0][0]);
  Q[2][1] = SQR_ABS(A[0][1]);
 
  // Calculate first eigenvector by the formula
  //   v[0] = conj( (A - w[0]).e1 x (A - w[0]).e2 )
  Q[0][0] = Q[0][1] + A[0][2]*w[0];
  Q[1][0] = Q[1][1] + A[1][2]*w[0];
  Q[2][0] = (real(A[0][0]) - w[0]) * (real(A[1][1]) - w[0]) - Q[2][1];
  norm    = SQR_ABS(Q[0][0]) + SQR_ABS(Q[1][0]) + SQR(real(Q[2][0]));
 
  // If vectors are nearly linearly dependent, or if there might have
  // been large cancellations in the calculation of A(I,I) - W(1), fall
  // back to QL algorithm
  // Note that this simultaneously ensures that multiple eigenvalues do
  // not cause problems: If W(1) = W(2), then A - W(1) * I has rank 1,
  // i.e. all columns of A - W(1) * I are linearly dependent.
  if (norm <= error)
    return zheevq3(A, Q, w);
  else                      // This is the standard branch
  {
    norm = sqrt(1.0 / norm);
    for (j=0; j < 3; j++)
      Q[j][0] = Q[j][0] * norm;
  }
 
  // Calculate second eigenvector by the formula
  //   v[1] = conj( (A - w[1]).e1 x (A - w[1]).e2 )
  Q[0][1]  = Q[0][1] + A[0][2]*w[1];
  Q[1][1]  = Q[1][1] + A[1][2]*w[1];
  Q[2][1]  = (real(A[0][0]) - w[1]) * (real(A[1][1]) - w[1]) - real(Q[2][1]);
  norm     = SQR_ABS(Q[0][1]) + SQR_ABS(Q[1][1]) + SQR(real(Q[2][1]));
  if (norm <= error)
    return zheevq3(A, Q, w);
  else
  {
    norm = sqrt(1.0 / norm);
    for (j=0; j < 3; j++)
      Q[j][1] = Q[j][1] * norm;
  }
 
  // Calculate third eigenvector according to
  //   v[2] = conj(v[0] x v[1])
  Q[0][2] = conj(Q[1][0]*Q[2][1] - Q[2][0]*Q[1][1]);
  Q[1][2] = conj(Q[2][0]*Q[0][1] - Q[0][0]*Q[2][1]);
  Q[2][2] = conj(Q[0][0]*Q[1][1] - Q[1][0]*Q[0][1]);
#endif
 
  return 0;
}

References SQR, SQR_ABS, zheevc3(), and zheevq3().

Referenced by OscProb::PMNS_DensityMatrix::SetInitialRho(), and OscProb::PMNS_Fast::SolveHamMatter().