zheevq3.cxx File Reference

#include <stdio.h>
#include <cmath>
#include <complex>
#include "zhetrd3.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	zheevq3 (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 26 of file zheevq3.cxx.

SQR_ABS

#define SQR_ABS

(

x

)

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

Definition at line 27 of file zheevq3.cxx.

Function Documentation

zheevq3()

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.

{
  const int n = 3;
  double e[3];                 // The third element is used only as temporary workspace
  double g, r, p, f, b, s, c;  // Intermediate storage
  std::complex<double> t;
  int nIter;
  int m;
 
  // Transform A to real tridiagonal form by the Householder method
  zhetrd3(A, Q, w, e);
 
  // Calculate eigensystem of the remaining real symmetric tridiagonal matrix
  // with the QL method
  //
  // Loop over all off-diagonal elements
  for (int l=0; l < n-1; l++)
  {
    nIter = 0;
    while (1)
    {
      // Check for convergence and exit iteration loop if off-diagonal
      // element e(l) is zero
      for (m=l; m <= n-2; m++)
      {
        g = fabs(w[m])+fabs(w[m+1]);
        if (fabs(e[m]) + g == g)
          break;
      }
      if (m == l)
        break;
 
      if (nIter++ >= 30)
        return -1;
 
      // Calculate g = d_m - k
      g = (w[l+1] - w[l]) / (e[l] + e[l]);
      r = sqrt(SQR(g) + 1.0);
      if (g > 0)
        g = w[m] - w[l] + e[l]/(g + r);
      else
        g = w[m] - w[l] + e[l]/(g - r);
 
      s = c = 1.0;
      p = 0.0;
      for (int i=m-1; i >= l; i--)
      {
        f = s * e[i];
        b = c * e[i];
        if (fabs(f) > fabs(g))
        {
          c      = g / f;
          r      = sqrt(SQR(c) + 1.0);
          e[i+1] = f * r;
          c     *= (s = 1.0/r);
        }
        else
        {
          s      = f / g;
          r      = sqrt(SQR(s) + 1.0);
          e[i+1] = g * r;
          s     *= (c = 1.0/r);
        }
 
        g = w[i+1] - p;
        r = (w[i] - g)*s + 2.0*c*b;
        p = s * r;
        w[i+1] = g + p;
        g = c*r - b;
 
        // Form eigenvectors
#ifndef EVALS_ONLY
        for (int k=0; k < n; k++)
        {
          t = Q[k][i+1];
          Q[k][i+1] = s*Q[k][i] + c*t;
          Q[k][i]   = c*Q[k][i] - s*t;
        }
#endif
      }
      w[l] -= p;
      e[l]  = g;
      e[m]  = 0.0;
    }
  }
 
  return 0;
}

References SQR, and zhetrd3().

Referenced by zheevh3().