Chemical Data Processing Library C++ API - Version 1.1.1
JacobiDiagonalization.hpp
Go to the documentation of this file.
1 /*
2  * JacobiDiagonalization.hpp
3  *
4  * This file is part of the Chemical Data Processing Toolkit
5  *
6  * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
7  *
8  * This library is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2 of the License, or (at your option) any later version.
12  *
13  * This library is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  * Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public License
19  * along with this library; see the file COPYING. If not, write to
20  * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
21  * Boston, MA 02111-1307, USA.
22  */
23 
30 #ifndef CDPL_MATH_JACOBIDIAGONALIZATION_HPP
31 #define CDPL_MATH_JACOBIDIAGONALIZATION_HPP
32 
33 #include <cstddef>
34 
35 #include "CDPL/Math/Vector.hpp"
36 #include "CDPL/Math/Matrix.hpp"
37 #include "CDPL/Math/CommonType.hpp"
38 #include "CDPL/Math/TypeTraits.hpp"
39 #include "CDPL/Base/Exceptions.hpp"
40 
41 
42 namespace CDPL
43 {
44 
45  namespace Math
46  {
47 
66  template <typename M1, typename V, typename M2>
67  bool jacobiDiagonalize(MatrixExpression<M1>& a, VectorExpression<V>& d, MatrixExpression<M2>& v, std::size_t max_iter = 50);
68  } // namespace Math
69 } // namespace CDPL
70 
71 
72 // Implementation
73 
74 namespace
75 {
76 
77  template <typename M>
78  void jacobiRotate(M& a, typename M::SizeType i, typename M::SizeType j,
79  typename M::SizeType k, typename M::SizeType l,
80  typename M::ValueType tau, typename M::ValueType s)
81  {
82  typedef typename M::ValueType ValueType;
83 
84  ValueType g = a(i, j);
85  ValueType h = a(k, l);
86 
87  a(i, j) = g - s * (h + g * tau);
88  a(k, l) = h + s * (g - h * tau);
89  }
90 } // namespace
91 
92 
93 template <typename M1, typename V, typename M2>
95 {
97  typename M2::ValueType>::Type ValueType;
98 
100  typename M2::SizeType>::Type SizeType;
101 
102  SizeType n = a().getSize1();
103 
104  CDPL_MATH_CHECK(n == a().getSize2() && n > 0 && SizeType(d().getSize()) >= n, "Preconditions violated", Base::SizeError);
105 
106  typename VectorTemporaryTraits<V>::Type b(n);
107  typename VectorTemporaryTraits<V>::Type z(n);
108 
109  v().assign(IdentityMatrix<ValueType>(n, n));
110 
111  for (SizeType ip = 0; ip < n; ip++) {
112  b(ip) = d()(ip) = a()(ip, ip); // Initialize b and d to the diagonal of a
113  z(ip) = ValueType(); // This vector will accumulate terms of the form tapq as in equation (11.1.14).
114  }
115 
116  for (std::size_t i = 0; i < max_iter; i++) {
117  ValueType sm = ValueType();
118 
119  for (SizeType ip = 0; ip < n - 1; ip++) // Sum off-diagonal elements.
120  for (SizeType iq = ip + 1; iq < n; iq++)
121  sm += TypeTraits<ValueType>::abs(a()(ip, iq));
122 
123  if (sm == ValueType()) // The normal return, which relies on quadratic convergence to machine underflow.
124  return true;
125 
126  ValueType tresh;
127 
128  if (i < 3)
129  tresh = ValueType(0.2) * sm / (n * n); // ...on the first three sweeps.
130  else
131  tresh = ValueType(); // ...thereafter.
132 
133  for (SizeType ip = 0; ip < n - 1; ip++) {
134  for (SizeType iq = ip + 1; iq < n; iq++) {
135  ValueType g = 100 * TypeTraits<ValueType>::abs(a()(ip, iq));
136 
137  if (i > 3 && (TypeTraits<ValueType>::abs(d()(ip)) + g) == TypeTraits<ValueType>::abs(d()(ip)) // After four sweeps, skip the rotation if the off-diagonal element is small.
138  && (TypeTraits<ValueType>::abs(d()(iq)) + g) == TypeTraits<ValueType>::abs(d()(iq)))
139  a()(ip, iq) = ValueType();
140 
141  else if (TypeTraits<ValueType>::abs(a()(ip, iq)) > tresh) {
142  ValueType h = d()(iq) - d()(ip);
143  ValueType t;
144 
146  t = (a()(ip, iq)) / h;
147 
148  else {
149  ValueType theta = 0.5 * h / a()(ip, iq);
150  t = 1 / (TypeTraits<ValueType>::abs(theta) + TypeTraits<ValueType>::sqrt(1 + theta * theta));
151 
152  if (theta < ValueType())
153  t = -t;
154  }
155 
156  ValueType c = 1 / TypeTraits<ValueType>::sqrt(1 + t * t);
157  ValueType s = t * c;
158  ValueType tau = s / (1 + c);
159  h = t * a()(ip, iq);
160  z(ip) -= h;
161  z(iq) += h;
162  d()(ip) -= h;
163  d()(iq) += h;
164  a()(ip, iq) = ValueType();
165 
166  for (SizeType j = 0; j < ip; j++) // Case of rotations 1 < j < p.
167  jacobiRotate(a(), j, ip, j, iq, tau, s);
168 
169  for (SizeType j = ip + 1; j < iq; j++) // Case of rotations p < j < q.
170  jacobiRotate(a(), ip, j, j, iq, tau, s);
171 
172  for (SizeType j = iq + 1; j < n; j++) // Case of rotations q < j < n.
173  jacobiRotate(a(), ip, j, iq, j, tau, s);
174 
175  for (SizeType j = 0; j < n; j++)
176  jacobiRotate(v(), j, ip, j, iq, tau, s);
177  }
178  }
179  }
180 
181  for (SizeType ip = 0; ip < n; ip++) {
182  b(ip) += z(ip);
183  d()(ip) = b(ip); // Update d with the sum of tapq,
184  z(ip) = ValueType(); // and reinitialize z.
185  }
186  }
187 
188  return false;
189 }
190 
191 #endif // CDPL_MATH_JACOBIDIAGONALIZATION_HPP
CDPL::Math::TypeTraits
Definition: TypeTraits.hpp:171
CommonType.hpp
Common type deduction.
CDPL::Math::VectorExpression
Definition: Expression.hpp:54
CDPL::Math::MatrixExpression
Definition: Expression.hpp:76
CDPL_MATH_CHECK
#define CDPL_MATH_CHECK(expr, msg, e)
Definition: Check.hpp:36
CDPL::Math::CommonType
Definition: CommonType.hpp:41
CDPL::Base::SizeError
Thrown to indicate that the size of a (multidimensional) array is not correct.
Definition: Base/Exceptions.hpp:133
CDPL::Chem::AtomType::M
const unsigned int M
A generic type that covers any element that is a metal.
Definition: AtomType.hpp:637
TypeTraits.hpp
Definition of type traits.
Exceptions.hpp
Definition of exception classes.
CDPL::Chem::CIPDescriptor::s
const unsigned int s
Specifies that the stereocenter has s configuration.
Definition: CIPDescriptor.hpp:81
CDPL
The namespace of the Chemical Data Processing Library.
CDPL::Math::IdentityMatrix
Definition: Matrix.hpp:1607
CDPL::Math::jacobiDiagonalize
bool jacobiDiagonalize(MatrixExpression< M1 > &a, VectorExpression< V > &d, MatrixExpression< M2 > &v, std::size_t max_iter=50)
Computes all eigenvalues and eigenvectors of a real symmetric matrix an using Jacobi's algorithm [WJA...
Definition: JacobiDiagonalization.hpp:94
Matrix.hpp
Definition of matrix data types.
CDPL::Math::VectorTemporaryTraits::Type
V::VectorTemporaryType Type
Definition: TypeTraits.hpp:181
Vector.hpp
Definition of vector data types.