zheevq3.h
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // Numerical diagonalization of 3x3 matrcies
3 // Copyright (C) 2006 Joachim Kopp
4 // ----------------------------------------------------------------------------
5 // This library is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU Lesser General Public
7 // License as published by the Free Software Foundation; either
8 // version 2.1 of the License, or (at your option) any later version.
9 //
10 // This library is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with this library; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 // ----------------------------------------------------------------------------
19 #ifndef __ZHEEVQ3_H
20 #define __ZHEEVQ3_H
21 
22 #include <complex>
23 
24 #include <stdio.h>
25 #include <math.h>
26 #include <complex>
27 #include "zhetrd3.h"
28 
29 // ----------------------------------------------------------------------------
30 template <typename T>
31 int zheevq3(std::complex<T> A[3][3], std::complex<T> Q[3][3], T w[3])
32 // ----------------------------------------------------------------------------
33 // Calculates the eigenvalues and normalized eigenvectors of a hermitian 3x3
34 // matrix A using the QL algorithm with implicit shifts, preceded by a
35 // Householder reduction to real tridiagonal form.
36 // The function accesses only the diagonal and upper triangular parts of A.
37 // The access is read-only.
38 // ----------------------------------------------------------------------------
39 // Parameters:
40 // A: The hermitian input matrix
41 // Q: Storage buffer for eigenvectors
42 // w: Storage buffer for eigenvalues
43 // ----------------------------------------------------------------------------
44 // Return value:
45 // 0: Success
46 // -1: Error (no convergence)
47 // ----------------------------------------------------------------------------
48 // Dependencies:
49 // zhetrd3()
50 // ----------------------------------------------------------------------------
51 {
52  const int n = 3;
53  T e[3]; // The third element is used only as temporary workspace
54  T g, r, p, f, b, s, c; // Intermediate storage
55  std::complex<T> t(0,0);
56  int nIter;
57  int m;
58 
59  // Transform A to real tridiagonal form by the Householder method
60  zhetrd3(A, Q, w, e);
61 
62  // Calculate eigensystem of the remaining real symmetric tridiagonal matrix
63  // with the QL method
64  //
65  // Loop over all off-diagonal elements
66  for (int l=0; l < n-1; l++)
67  {
68  nIter = 0;
69  while (1)
70  {
71  // Check for convergence and exit iteration loop if off-diagonal
72  // element e(l) is zero
73  for (m=l; m <= n-2; m++)
74  {
75  g = fabs(w[m])+fabs(w[m+1]);
76  if (fabs(e[m]) + g == g)
77  break;
78  }
79  if (m == l)
80  break;
81 
82  if (nIter++ >= 30)
83  return -1;
84 
85  // Calculate g = d_m - k
86  g = (w[l+1] - w[l]) / (e[l] + e[l]);
87  r = sqrt(SQR(g) + 1.0);
88  if (g > 0)
89  g = w[m] - w[l] + e[l]/(g + r);
90  else
91  g = w[m] - w[l] + e[l]/(g - r);
92 
93  s = c = 1.0;
94  p = 0.0;
95  for (int i=m-1; i >= l; i--)
96  {
97  f = s * e[i];
98  b = c * e[i];
99  if (fabs(f) > fabs(g))
100  {
101  c = g / f;
102  r = sqrt(SQR(c) + 1.0);
103  e[i+1] = f * r;
104  c *= (s = 1.0/r);
105  }
106  else
107  {
108  s = f / g;
109  r = sqrt(SQR(s) + 1.0);
110  e[i+1] = g * r;
111  s *= (c = 1.0/r);
112  }
113 
114  g = w[i+1] - p;
115  r = (w[i] - g)*s + 2.0*c*b;
116  p = s * r;
117  w[i+1] = g + p;
118  g = c*r - b;
119 
120  // Form eigenvectors
121 #ifndef EVALS_ONLY
122  for (int k=0; k < n; k++)
123  {
124  t = Q[k][i+1];
125  Q[k][i+1] = s*Q[k][i] + c*t;
126  Q[k][i] = c*Q[k][i] - s*t;
127  }
128 #endif
129  }
130  w[l] -= p;
131  e[l] = g;
132  e[m] = 0.0;
133  }
134  }
135 
136  return 0;
137 }
138 
139 #endif
int zheevq3(std::complex< T > A[3][3], std::complex< T > Q[3][3], T w[3])
Definition: zheevq3.h:31
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
const XML_Char * s
Definition: expat.h:262
std::void_t< T > n
void zhetrd3(std::complex< T > A[3][3], std::complex< T > Q[3][3], T d[3], T e[2])
Definition: zhetrd3.h:40
static const double A
Definition: Units.h:82
#define SQR(x)
Definition: TrackBasis.cxx:13
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
double T
Definition: Xdiff_gwt.C:5
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20