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

zhetrd3()

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.

{
  const int n = 3;
  std::complex<double> u[n], q[n];
  std::complex<double> omega, f;
  double K, h, g;
 
  // Initialize Q to the identitity matrix
#ifndef EVALS_ONLY
  for (int i=0; i < n; i++)
  {
    Q[i][i] = 1.0;
    for (int j=0; j < i; j++)
      Q[i][j] = Q[j][i] = 0.0;
  }
#endif
 
  // Bring first row and column to the desired form
  h = SQR_ABS(A[0][1]) + SQR_ABS(A[0][2]);
  if (real(A[0][1]) > 0)
    g = -sqrt(h);
  else
    g = sqrt(h);
  e[0] = g;
  f    = g * A[0][1];
  u[1] = conj(A[0][1]) - g;
  u[2] = conj(A[0][2]);
 
  omega = h - f;
  if (real(omega) > 0.0)
  {
    omega = 0.5 * (1.0 + conj(omega)/omega) / real(omega);
    K = 0.0;
    for (int i=1; i < n; i++)
    {
      f    = conj(A[1][i]) * u[1] + A[i][2] * u[2];
      q[i] = omega * f;                  // p
      K   += real(conj(u[i]) * f);      // u* A u
    }
    K *= 0.5 * SQR_ABS(omega);
 
    for (int i=1; i < n; i++)
      q[i] = q[i] - K * u[i];
 
    d[0] = real(A[0][0]);
    d[1] = real(A[1][1]) - 2.0*real(q[1]*conj(u[1]));
    d[2] = real(A[2][2]) - 2.0*real(q[2]*conj(u[2]));
 
    // Store inverse Householder transformation in Q
#ifndef EVALS_ONLY
    for (int j=1; j < n; j++)
    {
      f = omega * conj(u[j]);
      for (int i=1; i < n; i++)
        Q[i][j] = Q[i][j] - f*u[i];
    }
#endif
 
    // Calculate updated A[1][2] and store it in f
    f = A[1][2] - q[1]*conj(u[2]) - u[1]*conj(q[2]);
  }
  else
  {
    for (int i=0; i < n; i++)
      d[i] = real(A[i][i]);
    f = A[1][2];
  }
 
  // Make (23) element real
  e[1] = abs(f);
#ifndef EVALS_ONLY
  if (e[1] != 0.0)
  {
    f = conj(f) / e[1];
    for (int i=1; i < n; i++)
      Q[i][n-1] = Q[i][n-1] * f;
  }
#endif
}

References SQR_ABS.

Referenced by zheevq3().