zheevh3.h File Reference

#include <complex>

Go to the source code of this file.

Functions

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

Function Documentation

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.

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

Referenced by OscProb::PMNS_NSI::SolveHam(), and OscProb::PMNS_Fast::SolveHam().

00041              :
00042 //   A: The hermitian input matrix
00043 //   Q: Storage buffer for eigenvectors
00044 //   w: Storage buffer for eigenvalues
00045 // ----------------------------------------------------------------------------
00046 // Return value:
00047 //   0: Success
00048 //  -1: Error
00049 // ----------------------------------------------------------------------------
00050 // Dependencies:
00051 //   zheevc3(), zhetrd3(), zheevq3()
00052 // ----------------------------------------------------------------------------
00053 // Version history:
00054 //   v1.1: Simplified fallback condition --> speed-up
00055 //   v1.0: First released version
00056 // ----------------------------------------------------------------------------
00057 {
00058 #ifndef EVALS_ONLY
00059   double norm;          // Squared norm or inverse norm of current eigenvector
00060 //  double n0, n1;        // Norm of first and second columns of A
00061   double error;         // Estimated maximum roundoff error
00062   double t, u;          // Intermediate storage
00063   int j;                // Loop counter
00064 #endif
00065 
00066   // Calculate eigenvalues
00067   zheevc3(A, w);
00068 
00069 #ifndef EVALS_ONLY
00070 //  n0 = SQR(real(A[0][0])) + SQR_ABS(A[0][1]) + SQR_ABS(A[0][2]);
00071 //  n1 = SQR_ABS(A[0][1]) + SQR(real(A[1][1])) + SQR_ABS(A[1][2]);
00072   
00073   t = fabs(w[0]);
00074   if ((u=fabs(w[1])) > t)
00075     t = u;
00076   if ((u=fabs(w[2])) > t)
00077     t = u;
00078   if (t < 1.0)
00079     u = t;
00080   else
00081     u = SQR(t);
00082   error = 256.0 * DBL_EPSILON * SQR(u);
00083 //  error = 256.0 * DBL_EPSILON * (n0 + u) * (n1 + u);
00084 
00085   Q[0][1] = A[0][1]*A[1][2] - A[0][2]*real(A[1][1]);
00086   Q[1][1] = A[0][2]*conj(A[0][1]) - A[1][2]*real(A[0][0]);
00087   Q[2][1] = SQR_ABS(A[0][1]);
00088 
00089   // Calculate first eigenvector by the formula
00090   //   v[0] = conj( (A - w[0]).e1 x (A - w[0]).e2 )
00091   Q[0][0] = Q[0][1] + A[0][2]*w[0];
00092   Q[1][0] = Q[1][1] + A[1][2]*w[0];
00093   Q[2][0] = (real(A[0][0]) - w[0]) * (real(A[1][1]) - w[0]) - Q[2][1];
00094   norm    = SQR_ABS(Q[0][0]) + SQR_ABS(Q[1][0]) + SQR(real(Q[2][0]));
00095 
00096   // If vectors are nearly linearly dependent, or if there might have
00097   // been large cancellations in the calculation of A(I,I) - W(1), fall
00098   // back to QL algorithm
00099   // Note that this simultaneously ensures that multiple eigenvalues do
00100   // not cause problems: If W(1) = W(2), then A - W(1) * I has rank 1,
00101   // i.e. all columns of A - W(1) * I are linearly dependent.
00102   if (norm <= error)
00103     return zheevq3(A, Q, w);
00104   else                      // This is the standard branch
00105   {
00106     norm = sqrt(1.0 / norm);
00107     for (j=0; j < 3; j++)
00108       Q[j][0] = Q[j][0] * norm;
00109   }
00110   
00111   // Calculate second eigenvector by the formula
00112   //   v[1] = conj( (A - w[1]).e1 x (A - w[1]).e2 )
00113   Q[0][1]  = Q[0][1] + A[0][2]*w[1];
00114   Q[1][1]  = Q[1][1] + A[1][2]*w[1];
00115   Q[2][1]  = (real(A[0][0]) - w[1]) * (real(A[1][1]) - w[1]) - real(Q[2][1]);
00116   norm     = SQR_ABS(Q[0][1]) + SQR_ABS(Q[1][1]) + SQR(real(Q[2][1]));
00117   if (norm <= error)
00118     return zheevq3(A, Q, w);
00119   else
00120   {
00121     norm = sqrt(1.0 / norm);
00122     for (j=0; j < 3; j++)
00123       Q[j][1] = Q[j][1] * norm;
00124   }
00125   
00126   // Calculate third eigenvector according to
00127   //   v[2] = conj(v[0] x v[1])
00128   Q[0][2] = conj(Q[1][0]*Q[2][1] - Q[2][0]*Q[1][1]);
00129   Q[1][2] = conj(Q[2][0]*Q[0][1] - Q[0][0]*Q[2][1]);
00130   Q[2][2] = conj(Q[0][0]*Q[1][1] - Q[1][0]*Q[0][1]);
00131 #endif
00132 
00133   return 0;
00134 }


Generated on 22 Nov 2017 for loon by  doxygen 1.6.1