zhetrd3.h File Reference

#include <complex>

Go to the source code of this file.

Functions

void zhetrd3 (std::complex< double > A[3][3], std::complex< double > Q[3][3], double d[3], double e[2])

Function Documentation

void zhetrd3 ( std::complex< double >  A[3][3],
std::complex< double >  Q[3][3],
double  d[3],
double  e[2] 
)

Definition at line 30 of file zhetrd3.cxx.

References Munits::g, n, and SQR_ABS.

Referenced by zheevq3().

00034                                         :
00035 //            [ d[0]  e[0]       ]
00036 //    A = Q . [ e[0]  d[1]  e[1] ] . Q^T
00037 //            [       e[1]  d[2] ]
00038 // The function accesses only the diagonal and upper triangular parts of
00039 // A. The access is read-only.
00040 // ---------------------------------------------------------------------------
00041 {
00042   const int n = 3;
00043   std::complex<double> u[n], q[n];
00044   std::complex<double> omega, f;
00045   double K, h, g;
00046   
00047   // Initialize Q to the identitity matrix
00048 #ifndef EVALS_ONLY
00049   for (int i=0; i < n; i++)
00050   {
00051     Q[i][i] = 1.0;
00052     for (int j=0; j < i; j++)
00053       Q[i][j] = Q[j][i] = 0.0;
00054   }
00055 #endif
00056 
00057   // Bring first row and column to the desired form 
00058   h = SQR_ABS(A[0][1]) + SQR_ABS(A[0][2]);
00059   if (real(A[0][1]) > 0)
00060     g = -sqrt(h);
00061   else
00062     g = sqrt(h);
00063   e[0] = g;
00064   f    = g * A[0][1];
00065   u[1] = conj(A[0][1]) - g;
00066   u[2] = conj(A[0][2]);
00067   
00068   omega = h - f;
00069   if (real(omega) > 0.0)
00070   {
00071     omega = 0.5 * (1.0 + conj(omega)/omega) / real(omega);
00072     K = 0.0;
00073     for (int i=1; i < n; i++)
00074     {
00075       f    = conj(A[1][i]) * u[1] + A[i][2] * u[2];
00076       q[i] = omega * f;                  // p
00077       K   += real(conj(u[i]) * f);      // u* A u
00078     }
00079     K *= 0.5 * SQR_ABS(omega);
00080 
00081     for (int i=1; i < n; i++)
00082       q[i] = q[i] - K * u[i];
00083     
00084     d[0] = real(A[0][0]);
00085     d[1] = real(A[1][1]) - 2.0*real(q[1]*conj(u[1]));
00086     d[2] = real(A[2][2]) - 2.0*real(q[2]*conj(u[2]));
00087     
00088     // Store inverse Householder transformation in Q
00089 #ifndef EVALS_ONLY
00090     for (int j=1; j < n; j++)
00091     {
00092       f = omega * conj(u[j]);
00093       for (int i=1; i < n; i++)
00094         Q[i][j] = Q[i][j] - f*u[i];
00095     }
00096 #endif
00097 
00098     // Calculate updated A[1][2] and store it in f
00099     f = A[1][2] - q[1]*conj(u[2]) - u[1]*conj(q[2]);
00100   }
00101   else
00102   {
00103     for (int i=0; i < n; i++)
00104       d[i] = real(A[i][i]);
00105     f = A[1][2];
00106   }
00107 
00108   // Make (23) element real
00109   e[1] = abs(f);
00110 #ifndef EVALS_ONLY
00111   if (e[1] != 0.0)
00112   {
00113     f = conj(f) / e[1];
00114     for (int i=1; i < n; i++)
00115       Q[i][n-1] = Q[i][n-1] * f;
00116   }
00117 #endif
00118 }


Generated on 2 Nov 2017 for loon by  doxygen 1.6.1