herwig is hosted by Hepforge, IPPP Durham
Herwig 7.3.0
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
45namespace boost { namespace numeric { namespace ublas {
46
47template<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