Cantera  2.4.0
utilities.h
Go to the documentation of this file.
1 /**
2  * @file utilities.h
3  * Various templated functions that carry out common vector
4  * operations (see \ref utils).
5  */
6 
7 // This file is part of Cantera. See License.txt in the top-level directory or
8 // at http://www.cantera.org/license.txt for license and copyright information.
9 
10 /**
11  * @defgroup utils Templated Utility Functions
12  *
13  * These are templates to perform various simple operations on arrays. Note that
14  * the compiler will inline these, so using them carries no performance penalty.
15  */
16 
17 #ifndef CT_UTILITIES_H
18 #define CT_UTILITIES_H
19 
20 #include "ct_defs.h"
21 #include "global.h"
22 #include <stdexcept>
23 
24 #include <numeric>
25 
26 namespace Cantera
27 {
28 //! Unary operator to multiply the argument by a constant.
29 /*!
30  * The form of this operator is designed for use by std::transform.
31  * @see @ref scale().
32  */
33 template<class T> struct timesConstant : public std::unary_function<T, double> {
34  //! Constructor
35  /*!
36  * @param c Constant of templated type T that will be stored internally
37  * within the object and used in the multiplication operation
38  */
39  timesConstant(T c) : m_c(c) {}
40 
41  //! Parenthesis operator returning a double
42  /*!
43  * @param x Variable of templated type T that will be used in the
44  * multiplication operator
45  * @returns a value of type double from the internal multiplication
46  */
47  double operator()(T x) {
48  return m_c * x;
49  }
50 
51  //! Stored constant value of time T
52  T m_c;
53 };
54 
55 //! Templated Inner product of two vectors of length 4.
56 /*!
57  * If either \a x or \a y has length greater than 4, only the first 4 elements
58  * will be used.
59  *
60  * @param x first reference to the templated class V
61  * @param y second reference to the templated class V
62  * @return This class returns a hard-coded type, doublereal.
63  */
64 template<class V>
65 inline doublereal dot4(const V& x, const V& y)
66 {
67  return x[0]*y[0] + x[1]*y[1] + x[2]*y[2] + x[3]*y[3];
68 }
69 
70 //! Templated Inner product of two vectors of length 5
71 /*!
72  * If either \a x or \a y has length greater than 4, only the first 4 elements
73  * will be used.
74  *
75  * @param x first reference to the templated class V
76  * @param y second reference to the templated class V
77  * @return This class returns a hard-coded type, doublereal.
78  */
79 template<class V>
80 inline doublereal dot5(const V& x, const V& y)
81 {
82  return x[0]*y[0] + x[1]*y[1] + x[2]*y[2] + x[3]*y[3] +
83  x[4]*y[4];
84 }
85 
86 //! Function that calculates a templated inner product.
87 /*!
88  * This inner product is templated twice. The output variable is hard coded
89  * to return a doublereal.
90  *
91  * template<class InputIter, class InputIter2>
92  *
93  * @code
94  * double x[8], y[8];
95  * doublereal dsum = dot<double *,double *>(x, &x+7, y);
96  * @endcode
97  *
98  * @param x_begin Iterator pointing to the beginning, belonging to the
99  * iterator class InputIter.
100  * @param x_end Iterator pointing to the end, belonging to the
101  * iterator class InputIter.
102  * @param y_begin Iterator pointing to the beginning of y, belonging to the
103  * iterator class InputIter2.
104  * @return The return is hard-coded to return a double.
105  */
106 template<class InputIter, class InputIter2>
107 inline doublereal dot(InputIter x_begin, InputIter x_end,
108  InputIter2 y_begin)
109 {
110  return inner_product(x_begin, x_end, y_begin, 0.0);
111 }
112 
113 //! Multiply elements of an array by a scale factor.
114 /*!
115  * \code
116  * vector_fp in(8, 1.0), out(8);
117  * scale(in.begin(), in.end(), out.begin(), factor);
118  * \endcode
119  *
120  * @param begin Iterator pointing to the beginning, belonging to the
121  * iterator class InputIter.
122  * @param end Iterator pointing to the end, belonging to the
123  * iterator class InputIter.
124  * @param out Iterator pointing to the beginning of out, belonging to the
125  * iterator class OutputIter. This is the output variable
126  * for this routine.
127  * @param scale_factor input scale factor belonging to the class S.
128  */
129 template<class InputIter, class OutputIter, class S>
130 inline void scale(InputIter begin, InputIter end,
131  OutputIter out, S scale_factor)
132 {
133  std::transform(begin, end, out, timesConstant<S>(scale_factor));
134 }
135 
136 //! Multiply each entry in x by the corresponding entry in y.
137 /*!
138  * The template arguments are: template<class InputIter, class OutputIter>
139  *
140  * Simple code Equivalent:
141  * \code
142  * double x[10], y[10]
143  * for (n = 0; n < 10; n++) {
144  * x[n] *= y[n];
145  * }
146  * \endcode
147  * Example of function call usage to implement the simple code example:
148  * \code
149  * double x[10], y[10]
150  * multiply_each(x, x+10, y);
151  * \endcode
152  *
153  * @param x_begin Iterator pointing to the beginning of the vector x,
154  * belonging to the iterator class InputIter.
155  * @param x_end Iterator pointing to the end of the vector x, belonging to
156  * the iterator class InputIter. The difference between end and
157  * begin determines the loop length
158  * @param y_begin Iterator pointing to the beginning of the vector y,
159  * belonging to the iterator class outputIter.
160  */
161 template<class InputIter, class OutputIter>
162 inline void multiply_each(OutputIter x_begin, OutputIter x_end,
163  InputIter y_begin)
164 {
165  for (; x_begin != x_end; ++x_begin, ++y_begin) {
166  *x_begin *= *y_begin;
167  }
168 }
169 
170 //! The maximum absolute value (templated version)
171 /*!
172  * The template arguments are: template<class InputIter>
173  *
174  * Simple code Equivalent:
175  * \code
176  * double x[10] amax = 0.0;
177  * for (int n = 0; n < 10; n++) {
178  * if (fabs(x[n]) > amax) amax = fabs(x[10]);
179  * }
180  * return amax;
181  * \endcode
182  * Example of function call usage to implement the simple code example:
183  * \code
184  * double x[10]
185  * double amax = absmax(x, x+10);
186  * \endcode
187  *
188  * @param begin Iterator pointing to the beginning of the x vector,
189  * belonging to the iterator class InputIter.
190  * @param end Iterator pointing to the end of the x vector, belonging to
191  * the iterator class InputIter. The difference between end and
192  * begin determines the loop length
193  */
194 template<class InputIter>
195 inline doublereal absmax(InputIter begin, InputIter end)
196 {
197  doublereal amax = 0.0;
198  for (; begin != end; ++begin) {
199  amax = std::max(fabs(*begin), amax);
200  }
201  return amax;
202 }
203 
204 //! Normalize the values in a sequence, such that they sum to 1.0 (templated
205 //! version)
206 /*!
207  * The template arguments are: template<class InputIter, class OutputIter>
208  *
209  * Simple Equivalent:
210  * \code
211  * double x[10], y[10], sum = 0.0;
212  * for (int n = 0; n < 10; n++) {
213  * sum += x[10];
214  * }
215  * for (int n = 0; n < 10; n++) {
216  * y[n] = x[n]/sum;
217  * }
218  * \endcode
219  * Example of function call usage:
220  * \code
221  * double x[10], y[10];
222  * normalize(x, x+10, y);
223  * \endcode
224  *
225  * @param begin Iterator pointing to the beginning of the x vector,
226  * belonging to the iterator class InputIter.
227  * @param end Iterator pointing to the end of the x vector, belonging to
228  * the iterator class InputIter. The difference between end and
229  * begin determines the loop length
230  * @param out Iterator pointing to the beginning of the output vector,
231  * belonging to the iterator class OutputIter.
232  */
233 template<class InputIter, class OutputIter>
234 inline void normalize(InputIter begin, InputIter end,
235  OutputIter out)
236 {
237  doublereal sum = accumulate(begin, end, 0.0);
238  for (; begin != end; ++begin, ++out) {
239  *out = *begin/sum;
240  }
241 }
242 
243 //! Templated divide of each element of \a x by the corresponding element of \a y.
244 /*!
245  * The template arguments are: template<class InputIter, class OutputIter>
246  *
247  * Simple Equivalent:
248  * \code
249  * double x[10], y[10];
250  * for (n = 0; n < 10; n++) {
251  * x[n] /= y[n];
252  * }
253  * \endcode
254  * Example of code usage:
255  * \code
256  * double x[10], y[10];
257  * divide_each(x, x+10, y);
258  * \endcode
259  *
260  * @param x_begin Iterator pointing to the beginning of the x vector,
261  * belonging to the iterator class OutputIter.
262  * @param x_end Iterator pointing to the end of the x vector, belonging to
263  * the iterator class OutputIter. The difference between end
264  * and begin determines the number of inner iterations.
265  * @param y_begin Iterator pointing to the beginning of the yvector, belonging
266  * to the iterator class InputIter.
267  */
268 template<class InputIter, class OutputIter>
269 inline void divide_each(OutputIter x_begin, OutputIter x_end,
270  InputIter y_begin)
271 {
272  for (; x_begin != x_end; ++x_begin, ++y_begin) {
273  *x_begin /= *y_begin;
274  }
275 }
276 
277 //! Increment each entry in \a x by the corresponding entry in \a y.
278 /*!
279  * The template arguments are: template<class InputIter, class OutputIter>
280  *
281  * @param x_begin Iterator pointing to the beginning of the x vector,
282  * belonging to the iterator class OutputIter.
283  * @param x_end Iterator pointing to the end of the x vector, belonging to
284  * the iterator class OutputIter. The difference between end
285  * and begin determines the number of inner iterations.
286  * @param y_begin Iterator pointing to the beginning of the yvector, belonging
287  * to the iterator class InputIter.
288  */
289 template<class InputIter, class OutputIter>
290 inline void sum_each(OutputIter x_begin, OutputIter x_end,
291  InputIter y_begin)
292 {
293  for (; x_begin != x_end; ++x_begin, ++y_begin) {
294  *x_begin += *y_begin;
295  }
296 }
297 
298 //! Copies a contiguous range in a sequence to indexed
299 //! positions in another sequence.
300 /*!
301  * The template arguments are: template<class InputIter, class OutputIter, class IndexIter>
302  *
303  * Example:
304  *
305  * \code
306  * vector_fp x(3), y(20);
307  * vector_int index(3);
308  * index[0] = 9;
309  * index[1] = 2;
310  * index[3] = 16;
311  * scatter_copy(x.begin(), x.end(), y.begin(), index.begin());
312  * \endcode
313  *
314  * This routine is templated 3 times.
315  * InputIter is an iterator for the source vector
316  * OutputIter is an iterator for the destination vector
317  * IndexIter is an iterator for the index into the destination vector.
318  *
319  * @param begin Iterator pointing to the beginning of the source vector,
320  * belonging to the iterator class InputIter.
321  * @param end Iterator pointing to the end of the source vector, belonging
322  * to the iterator class InputIter. The difference between end
323  * and begin determines the number of inner iterations.
324  * @param result Iterator pointing to the beginning of the output vector,
325  * belonging to the iterator class outputIter.
326  * @param index Iterator pointing to the beginning of the index vector, belonging to the
327  * iterator class IndexIter.
328  */
329 template<class InputIter, class OutputIter, class IndexIter>
330 inline void scatter_copy(InputIter begin, InputIter end,
331  OutputIter result, IndexIter index)
332 {
333  for (; begin != end; ++begin, ++index) {
334  *(result + *index) = *begin;
335  }
336 }
337 
338 //! Multiply selected elements in an array by a contiguous sequence of
339 //! multipliers.
340 /*!
341  * The template arguments are: template<class InputIter, class RandAccessIter, class IndexIter>
342  *
343  * Example:
344  * \code
345  * double multipliers[] = {8.9, -2.0, 5.6};
346  * int index[] = {7, 4, 13};
347  * vector_fp data(20);
348  * ...
349  * // Multiply elements 7, 4, and 13 in data by multipliers[0], multipliers[1],and multipliers[2],
350  * // respectively
351  * scatter_mult(multipliers, multipliers + 3, data.begin(), index);
352  * \endcode
353  *
354  * @param mult_begin Iterator pointing to the beginning of the multiplier
355  * vector, belonging to the iterator class InputIter.
356  * @param mult_end Iterator pointing to the end of the multiplier vector,
357  * belonging to the iterator class InputIter. The difference
358  * between end and begin determines the number of inner
359  * iterations.
360  * @param data Iterator pointing to the beginning of the output vector,
361  * belonging to the iterator class RandAccessIter, that will
362  * be selectively multiplied.
363  * @param index Iterator pointing to the beginning of the index vector,
364  * belonging to the iterator class IndexIter.
365  */
366 template<class InputIter, class RandAccessIter, class IndexIter>
367 inline void scatter_mult(InputIter mult_begin, InputIter mult_end,
368  RandAccessIter data, IndexIter index)
369 {
370  for (; mult_begin != mult_end; ++mult_begin, ++index) {
371  *(data + *index) *= *mult_begin;
372  }
373 }
374 
375 //! Compute \f[ \sum_k x_k \log x_k. \f].
376 /*!
377  * The template arguments are: template<class InputIter>
378  *
379  * A small number (1.0E-20) is added before taking the log. This templated
380  * class does the indicated sun. The template must be an iterator.
381  *
382  * @param begin Iterator pointing to the beginning, belonging to the
383  * iterator class InputIter.
384  * @param end Iterator pointing to the end, belonging to the
385  * iterator class InputIter.
386  * @return The return from this class is a double.
387  */
388 template<class InputIter>
389 inline doublereal sum_xlogx(InputIter begin, InputIter end)
390 {
391  doublereal sum = 0.0;
392  for (; begin != end; ++begin) {
393  sum += (*begin) * std::log(*begin + Tiny);
394  }
395  return sum;
396 }
397 
398 //! Compute \f[ \sum_k x_k \log Q_k. \f].
399 /*!
400  * The template arguments are: template<class InputIter1, class InputIter2>
401  *
402  * This class is templated twice. The first template, InputIter1 is the iterator
403  * that points to $x_k$. The second iterator InputIter2, point to $Q_k$. A small
404  * number (1.0E-20) is added before taking the log.
405  *
406  * @param begin Iterator pointing to the beginning, belonging to the
407  * iterator class InputIter1.
408  * @param end Iterator pointing to the end, belonging to the
409  * iterator class InputIter1.
410  * @param Q_begin Iterator pointing to the beginning of Q_k, belonging to the
411  * iterator class InputIter2.
412  * @return The return from this class is hard coded to a doublereal.
413  */
414 template<class InputIter1, class InputIter2>
415 inline doublereal sum_xlogQ(InputIter1 begin, InputIter1 end,
416  InputIter2 Q_begin)
417 {
418  doublereal sum = 0.0;
419  for (; begin != end; ++begin, ++Q_begin) {
420  sum += (*begin) * std::log(*Q_begin + Tiny);
421  }
422  return sum;
423 }
424 
425 //! Templated evaluation of a polynomial of order 6
426 /*!
427  * @param x Value of the independent variable - First template parameter
428  * @param c Pointer to the polynomial - Second template parameter
429  */
430 template<class D, class R>
431 R poly6(D x, R* c)
432 {
433  return ((((((c[6]*x + c[5])*x + c[4])*x + c[3])*x +
434  c[2])*x + c[1])*x + c[0]);
435 }
436 
437 //! Templated evaluation of a polynomial of order 8
438 /*!
439  * @param x Value of the independent variable - First template parameter
440  * @param c Pointer to the polynomial - Second template parameter
441  */
442 template<class D, class R>
443 R poly8(D x, R* c)
444 {
445  return ((((((((c[8]*x + c[7])*x + c[6])*x + c[5])*x + c[4])*x + c[3])*x +
446  c[2])*x + c[1])*x + c[0]);
447 }
448 
449 //! Templated evaluation of a polynomial of order 5
450 /*!
451  * @param x Value of the independent variable - First template parameter
452  * @param c Pointer to the polynomial - Second template parameter
453  */
454 template<class D, class R>
455 R poly5(D x, R* c)
456 {
457  return (((((c[5]*x + c[4])*x + c[3])*x +
458  c[2])*x + c[1])*x + c[0]);
459 }
460 
461 //! Evaluates a polynomial of order 4.
462 /*!
463  * @param x Value of the independent variable.
464  * @param c Pointer to the polynomial coefficient array.
465  */
466 template<class D, class R>
467 R poly4(D x, R* c)
468 {
469  return ((((c[4]*x + c[3])*x +
470  c[2])*x + c[1])*x + c[0]);
471 }
472 
473 //! Templated evaluation of a polynomial of order 3
474 /*!
475  * @param x Value of the independent variable - First template parameter
476  * @param c Pointer to the polynomial - Second template parameter
477  */
478 template<class D, class R>
479 R poly3(D x, R* c)
480 {
481  return (((c[3]*x + c[2])*x + c[1])*x + c[0]);
482 }
483 
484 //@}
485 
486 //! Check to see that a number is finite (not NaN, +Inf or -Inf)
487 void checkFinite(const double tmp);
488 
489 //! Check to see that all elements in an array are finite
490 /*!
491  * Throws an exception if any element is NaN, +Inf, or -Inf
492  * @param name Name to be used in the exception message if the check fails
493  * @param values Array of *N* values to be checked
494  * @param N Number of elements in *values*
495  */
496 void checkFinite(const std::string& name, double* values, size_t N);
497 
498 //! Const accessor for a value in a std::map.
499 /*
500  * Similar to std::map.at(key), but returns *default_val* if the key is not
501  * found instead of throwing an exception.
502  */
503 template <class T, class U>
504 const U& getValue(const std::map<T, U>& m, const T& key, const U& default_val) {
505  typename std::map<T,U>::const_iterator iter = m.find(key);
506  return (iter == m.end()) ? default_val : iter->second;
507 }
508 
509 }
510 
511 #endif
R poly3(D x, R *c)
Templated evaluation of a polynomial of order 3.
Definition: utilities.h:479
void divide_each(OutputIter x_begin, OutputIter x_end, InputIter y_begin)
Templated divide of each element of x by the corresponding element of y.
Definition: utilities.h:269
timesConstant(T c)
Constructor.
Definition: utilities.h:39
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
Definition: checkFinite.cpp:15
doublereal dot4(const V &x, const V &y)
Templated Inner product of two vectors of length 4.
Definition: utilities.h:65
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
Definition: utilities.h:455
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
Unary operator to multiply the argument by a constant.
Definition: utilities.h:33
double operator()(T x)
Parenthesis operator returning a double.
Definition: utilities.h:47
void multiply_each(OutputIter x_begin, OutputIter x_end, InputIter y_begin)
Multiply each entry in x by the corresponding entry in y.
Definition: utilities.h:162
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition: utilities.h:467
R poly8(D x, R *c)
Templated evaluation of a polynomial of order 8.
Definition: utilities.h:443
void normalize(InputIter begin, InputIter end, OutputIter out)
Normalize the values in a sequence, such that they sum to 1.0 (templated version) ...
Definition: utilities.h:234
R poly6(D x, R *c)
Templated evaluation of a polynomial of order 6.
Definition: utilities.h:431
void scatter_mult(InputIter mult_begin, InputIter mult_end, RandAccessIter data, IndexIter index)
Multiply selected elements in an array by a contiguous sequence of multipliers.
Definition: utilities.h:367
doublereal absmax(InputIter begin, InputIter end)
The maximum absolute value (templated version)
Definition: utilities.h:195
const U & getValue(const std::map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a std::map.
Definition: utilities.h:504
doublereal sum_xlogQ(InputIter1 begin, InputIter1 end, InputIter2 Q_begin)
Compute .
Definition: utilities.h:415
doublereal sum_xlogx(InputIter begin, InputIter end)
Compute .
Definition: utilities.h:389
doublereal dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition: utilities.h:107
void sum_each(OutputIter x_begin, OutputIter x_end, InputIter y_begin)
Increment each entry in x by the corresponding entry in y.
Definition: utilities.h:290
doublereal dot5(const V &x, const V &y)
Templated Inner product of two vectors of length 5.
Definition: utilities.h:80
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:130
const doublereal Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:143
void scatter_copy(InputIter begin, InputIter end, OutputIter result, IndexIter index)
Copies a contiguous range in a sequence to indexed positions in another sequence. ...
Definition: utilities.h:330
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
T m_c
Stored constant value of time T.
Definition: utilities.h:52