19 std::string error =
"Error: Matrix-Vector dimensions are not compatible: " +
20 boost::lexical_cast<std::string>(x_o.
NumRows()) +
21 " rows vs " + boost::lexical_cast<std::string>(n) +
27 double alpha = MaxAbsDiag(A);
38 x_t = GetSolutionVector(B, x_o, alpha, t, tol);
44double UniformTaylor::MaxAbsDiag(
const Matrix& A) {
46 double a_ii = A(1, 1);
51 double max_a_ii = a_ii;
54 for (
int i = 2; i <= n; ++i) {
58 if (a_ii > max_a_ii) {
68 double alpha,
double t,
double tol) {
70 long double alpha_t =
alpha * t;
71 long double expat = exp(-alpha_t);
75 "Error: exp(-alpha * t) exceeds the range of a long double.";
76 error +=
"\nThe Uniform Taylor method cannot solve the matrix exponential.";
77 throw ValueError(error);
86 Vector C_next = expat * x_o;
91 int maxTerms = MaxNumTerms(alpha_t, tol);
95 for (
int k = 1; k < maxTerms; ++k) {
112int UniformTaylor::MaxNumTerms(
long double alpha_t,
double epsilon) {
113 long double nextTerm;
116 long double prevTerm = 1;
117 long double sumTerms = 1;
120 long double lowerBound = exp(alpha_t);
123 if (lowerBound == HUGE_VAL) {
125 "Error: exp(alpha * t) exceeds the range of a long double";
126 error +=
"\nThe Uniform Taylor method cannot solve the matrix exponential.";
127 throw ValueError(error);
130 lowerBound = lowerBound * (1 - epsilon);
134 bool stopSum =
false;
136 while (stopSum !=
true) {
138 if (sumTerms >= lowerBound) {
142 nextTerm = (alpha_t / p) * prevTerm;
146 sumTerms += nextTerm;
For values that are too big, too small, etc.
taken directly from OsiSolverInterface.cpp on 2/17/14 from https://projects.coin-or....
struct pyne::alpha alpha
a struct matching the '/decay/alphas' table in nuc_data.h5.