zheevq3.h File Reference

#include <complex>

Go to the source code of this file.

Functions

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

Function Documentation

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

Definition at line 31 of file zheevq3.cxx.

References MuELoss::e, Munits::g, Munits::m, n, SQR, and zhetrd3().

Referenced by zheevh3().

00039              :
00040 //   A: The hermitian input matrix
00041 //   Q: Storage buffer for eigenvectors
00042 //   w: Storage buffer for eigenvalues
00043 // ----------------------------------------------------------------------------
00044 // Return value:
00045 //   0: Success
00046 //  -1: Error (no convergence)
00047 // ----------------------------------------------------------------------------
00048 // Dependencies:
00049 //   zhetrd3()
00050 // ----------------------------------------------------------------------------
00051 {
00052   const int n = 3;
00053   double e[3];                 // The third element is used only as temporary workspace
00054   double g, r, p, f, b, s, c;  // Intermediate storage
00055   std::complex<double> t;
00056   int nIter;
00057   int m;
00058 
00059   // Transform A to real tridiagonal form by the Householder method
00060   zhetrd3(A, Q, w, e);
00061   
00062   // Calculate eigensystem of the remaining real symmetric tridiagonal matrix
00063   // with the QL method
00064   //
00065   // Loop over all off-diagonal elements
00066   for (int l=0; l < n-1; l++)
00067   {
00068     nIter = 0;
00069     while (1)
00070     {
00071       // Check for convergence and exit iteration loop if off-diagonal
00072       // element e(l) is zero
00073       for (m=l; m <= n-2; m++)
00074       {
00075         g = fabs(w[m])+fabs(w[m+1]);
00076         if (fabs(e[m]) + g == g)
00077           break;
00078       }
00079       if (m == l)
00080         break;
00081       
00082       if (nIter++ >= 30)
00083         return -1;
00084 
00085       // Calculate g = d_m - k
00086       g = (w[l+1] - w[l]) / (e[l] + e[l]);
00087       r = sqrt(SQR(g) + 1.0);
00088       if (g > 0)
00089         g = w[m] - w[l] + e[l]/(g + r);
00090       else
00091         g = w[m] - w[l] + e[l]/(g - r);
00092 
00093       s = c = 1.0;
00094       p = 0.0;
00095       for (int i=m-1; i >= l; i--)
00096       {
00097         f = s * e[i];
00098         b = c * e[i];
00099         if (fabs(f) > fabs(g))
00100         {
00101           c      = g / f;
00102           r      = sqrt(SQR(c) + 1.0);
00103           e[i+1] = f * r;
00104           c     *= (s = 1.0/r);
00105         }
00106         else
00107         {
00108           s      = f / g;
00109           r      = sqrt(SQR(s) + 1.0);
00110           e[i+1] = g * r;
00111           s     *= (c = 1.0/r);
00112         }
00113         
00114         g = w[i+1] - p;
00115         r = (w[i] - g)*s + 2.0*c*b;
00116         p = s * r;
00117         w[i+1] = g + p;
00118         g = c*r - b;
00119 
00120         // Form eigenvectors
00121 #ifndef EVALS_ONLY
00122         for (int k=0; k < n; k++)
00123         {
00124           t = Q[k][i+1];
00125           Q[k][i+1] = s*Q[k][i] + c*t;
00126           Q[k][i]   = c*Q[k][i] - s*t;
00127         }
00128 #endif 
00129       }
00130       w[l] -= p;
00131       e[l]  = g;
00132       e[m]  = 0.0;
00133     }
00134   }
00135 
00136   return 0;
00137 }


Generated on 22 Nov 2017 for loon by  doxygen 1.6.1