herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
expm-1.h
1 //
2 // Copyright (c) 2007
3 // Tsai, Dung-Bang
4 // National Taiwan University, Department of Physics
5 //
6 // E-Mail : dbtsai (at) gmail.com
7 // Begine : 2007/11/20
8 // Last modify : 2007/11/22
9 // Version : v0.1
10 //
11 // EXPGM_PAD computes the matrix exponential exp(H) for general matrixs,
12 // including complex and real matrixs using the irreducible (p,p) degree
13 // rational Pade approximation to the exponential
14 // exp(z) = r(z)=(+/-)( I+2*(Q(z)/P(z))).
15 //
16 // Usage :
17 //
18 // U = expm_pad(H)
19 // U = expm_pad(H, p)
20 //
21 // where p is internally set to 6 (recommended and gererally satisfactory).
22 //
23 // See also MATLAB supplied functions, EXPM and EXPM1.
24 //
25 // Reference :
26 // EXPOKIT, Software Package for Computing Matrix Exponentials.
27 // ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
28 //
29 // Permission to use, copy, modify, distribute and sell this software
30 // and its documentation for any purpose is hereby granted without fee,
31 // provided that the above copyright notice appear in all copies and
32 // that both that copyright notice and this permission notice appear
33 // in supporting documentation. The authors make no representations
34 // about the suitability of this software for any purpose.
35 // It is provided "as is" without express or implied warranty.
36 //
37 
38 #ifndef _BOOST_UBLAS_EXPM_
39 #define _BOOST_UBLAS_EXPM_
40 #include <complex>
41 #include <boost/numeric/ublas/vector.hpp>
42 #include <boost/numeric/ublas/matrix.hpp>
43 #include <boost/numeric/ublas/lu.hpp>
44 
45 namespace boost { namespace numeric { namespace ublas {
46 
47 template<typename MATRIX> MATRIX expm_pad(const MATRIX &H, const int p = 6) {
48  typedef typename MATRIX::value_type value_type;
49  typedef typename MATRIX::size_type size_type;
50  typedef double real_value_type; // Correct me. Need to modify.
51  assert(H.size1() == H.size2());
52  const size_type n = H.size1();
53  const identity_matrix<value_type> I(n);
54  matrix<value_type> U(n,n),H2(n,n),P(n,n),Q(n,n);
55  real_value_type norm = 0.0;
56 
57  // Calcuate Pade coefficients (1-based instead of 0-based as in the c vector)
58  vector<real_value_type> c(p+2);
59  c(1)=1;
60  for(size_type i = 1; i <= p; ++i)
61  c(i+1) = c(i) * ((p + 1.0 - i)/(i * (2.0 * p + 1 - i)));
62  // Calcuate the infinty norm of H, which is defined as the largest row sum of a matrix
63  for(size_type i=0; i<n; ++i)
64  {
65  real_value_type temp = 0.0;
66  for(size_type j=0;j<n;j++)
67  temp += std::abs<real_value_type>(H(j,i)); // Correct me, if H is complex, can I use that abs?
68  norm = std::max<real_value_type>(norm, temp);
69  }
70  if (norm == 0.0)
71  {
72  boost::throw_exception(boost::numeric::ublas::bad_argument());
73  std::cerr<<"Error! Null input in the routine EXPM_PAD.\n";
74  exit(0);
75  }
76  // Scaling, seek s such that || H*2^(-s) || < 1/2, and set scale = 2^(-s)
77  int s = 0;
78  real_value_type scale = 1.0;
79  if(norm > 0.5) {
80  s = std::max<int>(0, static_cast<int>((log(norm) / log(2.0) + 2.0)));
81  scale /= static_cast<real_value_type>(std::pow(2.0, s));
82  U.assign(scale * H); // Here U is used as temp value due to that H is const
83  }
84  else
85  U.assign(H);
86  // Horner evaluation of the irreducible fraction, see the following ref above.
87  // Initialise P (numerator) and Q (denominator)
88  H2.assign( prod(U, U) );
89  Q.assign( c(p+1)*I );
90  P.assign( c(p)*I );
91  size_type odd = 1;
92  for( size_type k = p - 1; k > 0; --k)
93  {
94  if( odd == 1)
95  {
96  Q = ( prod(Q, H2) + c(k) * I );
97  }
98  else
99  {
100  P = ( prod(P, H2) + c(k) * I );
101  }
102  odd = 1 - odd;
103  }
104  if( odd == 1)
105  {
106  Q = ( prod(Q, U) );
107  Q -= P ;
108  //U.assign( -(I + 2*(Q\P)));
109  }
110  else
111  {
112  P = (prod(P, U));
113  Q -= P;
114  //U.assign( I + 2*(Q\P));
115  }
116  // In origine expokit package, they use lapack ZGESV to obtain inverse matrix,
117  // and in that ZGESV routine, it uses LU decomposition for obtaing inverse matrix.
118  // Since in ublas, there is no matrix inversion template, I simply use the build-in
119  // LU decompostion package in ublas, and back substitute by myself.
120  //
122  permutation_matrix<size_type> pm(n);
123  int res = lu_factorize(Q, pm);
124  if( res != 0)
125  {
126  std::cerr << "Error in the matrix inversion in template expm_pad.\n";
127  exit(0);
128  }
129  H2 = I; // H2 is not needed anymore, so it is temporary used as identity matrix for substituting.
130 
131  lu_substitute(Q, pm, H2);
132  if( odd == 1)
133  U.assign( -(I + 2.0 * prod(H2, P)));
134  else
135  U.assign( I + 2.0 * prod(H2, P));
136  // Squaring
137  for(size_t i = 0; i < s; ++i)
138  {
139  U = (prod(U,U));
140  }
141  return U;
142  }
143 
144 }}}
145 
146 
147 #endif
Definition: expm-1.h:45