GetFEM  5.4.2
getfem_generic_assembly_tree.h
Go to the documentation of this file.
1 /*===========================================================================
2 
3  Copyright (C) 2013-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /** @file getfem_generic_assembly_tree.h
32  @author Yves Renard <Yves.Renard@insa-lyon.fr>
33  @date November 18, 2013.
34  @brief Definition of the syntax tree and basic operations on it.
35  Internal header for the generic assembly language part.
36  */
37 
38 
39 #ifndef GETFEM_GENERIC_ASSEMBLY_TREE_H__
40 #define GETFEM_GENERIC_ASSEMBLY_TREE_H__
41 
42 #include "getfem/getfem_generic_assembly.h"
43 #include "getfem/getfem_models.h"
44 #include "gmm/gmm_blas.h"
45 #include <iomanip>
46 #include "getfem/getfem_omp.h"
47 #include "getfem/dal_singleton.h"
48 #include "getfem/bgeot_rtree.h"
49 #include "getfem/bgeot_geotrans_inv.h"
50 #include "getfem/getfem_copyable_ptr.h"
51 #ifndef _WIN32
52 extern "C"{
53 #include <unistd.h>
54 }
55 #endif
56 
57 #define GA_DEBUG_ASSERT(a, b) GMM_ASSERT1(a, b)
58 // #define GA_DEBUG_ASSERT(a, b)
59 
60 #if 1
61 # define GA_TIC
62 # define GA_TOC(a)
63 # define GA_TOCTIC(a)
64 #else
65 # define GA_TIC scalar_type _ga_time_ = gmm::uclock_sec();
66 # define GA_TOC(a) { cout <<(a)<<" : "<<gmm::uclock_sec()-_ga_time_<< endl; }
67 # define GA_TOCTIC(a) { GA_TOC(a); _ga_time_ = gmm::uclock_sec(); }
68 #endif
69 
70 namespace getfem {
71 
72  // Basic token types (basic language components)
73  enum GA_TOKEN_TYPE {
74  GA_INVALID = 0, // invalid token
75  GA_END, // string end
76  GA_NAME, // A variable or user defined nonlinear function name
77  GA_SCALAR, // A real number
78  GA_PLUS, // '+'
79  GA_MINUS, // '-'
80  GA_UNARY_MINUS, // '-'
81  GA_MULT, // '*'
82  GA_DIV, // '/'
83  GA_COLON, // ':'
84  GA_QUOTE, // ''' transpose
85  GA_COLON_EQ, // ':=' macro def
86  GA_DEF, // 'Def' macro def
87  GA_SYM, // 'Sym(M)' operator
88  GA_SKEW, // 'Skew(M)' operator
89  GA_TRACE, // 'Trace(M)' operator
90  GA_DEVIATOR, // 'Deviator' operator
91  GA_INTERPOLATE, // 'Interpolate' operation
92  GA_INTERPOLATE_FILTER, // 'Interpolate_filter' operation
93  GA_INTERPOLATE_DERIVATIVE, // 'Interpolate_derivative' operation
94  GA_ELEMENTARY, // 'Elementary' operation (operation at the element level)
95  GA_SECONDARY_DOMAIN, // For the integration on a product of two domains
96  GA_XFEM_PLUS, // Evaluation on the + side of a level-set for fem_level_set
97  GA_XFEM_MINUS, // Evaluation on the - side of a level-set for fem_level_set
98  GA_PRINT, // 'Print' Print the tensor
99  GA_DOT, // '.'
100  GA_DOTMULT, // '.*' componentwise multiplication
101  GA_DOTDIV, // './' componentwise division
102  GA_TMULT, // '@' tensor product
103  GA_COMMA, // ','
104  GA_DCOMMA, // ',,'
105  GA_SEMICOLON, // ';'
106  GA_DSEMICOLON, // ';;'
107  GA_LPAR, // '('
108  GA_RPAR, // ')'
109  GA_LBRACKET, // '['
110  GA_RBRACKET, // ']'
111  GA_NB_TOKEN_TYPE
112  };
113 
114  // Detects Grad_, Hess_ or Div_
115  size_type ga_parse_prefix_operator(std::string &name);
116  // Detects Test_ and Test2_
117  size_type ga_parse_prefix_test(std::string &name);
118 
119  // Types of nodes for the syntax tree
120  enum GA_NODE_TYPE {
121  GA_NODE_VOID = 0,
122  GA_NODE_OP,
123  GA_NODE_PREDEF_FUNC,
124  GA_NODE_SPEC_FUNC,
125  GA_NODE_OPERATOR,
126  GA_NODE_CONSTANT,
127  GA_NODE_NAME,
128  GA_NODE_MACRO_PARAM,
129  GA_NODE_PARAMS,
130  GA_NODE_RESHAPE,
131  GA_NODE_CROSS_PRODUCT,
132  GA_NODE_SWAP_IND,
133  GA_NODE_IND_MOVE_LAST,
134  GA_NODE_CONTRACT,
135  GA_NODE_ALLINDICES,
136  GA_NODE_C_MATRIX,
137  GA_NODE_X,
138  GA_NODE_ELT_SIZE,
139  GA_NODE_ELT_K,
140  GA_NODE_ELT_B,
141  GA_NODE_NORMAL,
142  GA_NODE_VAL,
143  GA_NODE_GRAD,
144  GA_NODE_HESS,
145  GA_NODE_DIVERG,
146  GA_NODE_VAL_TEST,
147  GA_NODE_GRAD_TEST,
148  GA_NODE_HESS_TEST,
149  GA_NODE_DIVERG_TEST,
150  GA_NODE_INTERPOLATE,
151  GA_NODE_INTERPOLATE_FILTER,
152  GA_NODE_INTERPOLATE_VAL,
153  GA_NODE_INTERPOLATE_GRAD,
154  GA_NODE_INTERPOLATE_HESS,
155  GA_NODE_INTERPOLATE_DIVERG,
156  GA_NODE_INTERPOLATE_VAL_TEST,
157  GA_NODE_INTERPOLATE_GRAD_TEST,
158  GA_NODE_INTERPOLATE_HESS_TEST,
159  GA_NODE_INTERPOLATE_DIVERG_TEST,
160  GA_NODE_INTERPOLATE_X,
161  GA_NODE_INTERPOLATE_ELT_K,
162  GA_NODE_INTERPOLATE_ELT_B,
163  GA_NODE_INTERPOLATE_NORMAL,
164  GA_NODE_INTERPOLATE_DERIVATIVE,
165  GA_NODE_ELEMENTARY,
166  GA_NODE_ELEMENTARY_VAL,
167  GA_NODE_ELEMENTARY_GRAD,
168  GA_NODE_ELEMENTARY_HESS,
169  GA_NODE_ELEMENTARY_DIVERG,
170  GA_NODE_ELEMENTARY_VAL_TEST,
171  GA_NODE_ELEMENTARY_GRAD_TEST,
172  GA_NODE_ELEMENTARY_HESS_TEST,
173  GA_NODE_ELEMENTARY_DIVERG_TEST,
174  GA_NODE_SECONDARY_DOMAIN,
175  GA_NODE_SECONDARY_DOMAIN_VAL,
176  GA_NODE_SECONDARY_DOMAIN_GRAD,
177  GA_NODE_SECONDARY_DOMAIN_HESS,
178  GA_NODE_SECONDARY_DOMAIN_DIVERG,
179  GA_NODE_SECONDARY_DOMAIN_VAL_TEST,
180  GA_NODE_SECONDARY_DOMAIN_GRAD_TEST,
181  GA_NODE_SECONDARY_DOMAIN_HESS_TEST,
182  GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST,
183  GA_NODE_SECONDARY_DOMAIN_X,
184  GA_NODE_SECONDARY_DOMAIN_NORMAL,
185  GA_NODE_XFEM_PLUS,
186  GA_NODE_XFEM_PLUS_VAL,
187  GA_NODE_XFEM_PLUS_GRAD,
188  GA_NODE_XFEM_PLUS_HESS,
189  GA_NODE_XFEM_PLUS_DIVERG,
190  GA_NODE_XFEM_PLUS_VAL_TEST,
191  GA_NODE_XFEM_PLUS_GRAD_TEST,
192  GA_NODE_XFEM_PLUS_HESS_TEST,
193  GA_NODE_XFEM_PLUS_DIVERG_TEST,
194  GA_NODE_XFEM_MINUS,
195  GA_NODE_XFEM_MINUS_VAL,
196  GA_NODE_XFEM_MINUS_GRAD,
197  GA_NODE_XFEM_MINUS_HESS,
198  GA_NODE_XFEM_MINUS_DIVERG,
199  GA_NODE_XFEM_MINUS_VAL_TEST,
200  GA_NODE_XFEM_MINUS_GRAD_TEST,
201  GA_NODE_XFEM_MINUS_HESS_TEST,
202  GA_NODE_XFEM_MINUS_DIVERG_TEST,
203  GA_NODE_ZERO};
204 
205  typedef std::shared_ptr<std::string> pstring;
206  // Print error message indicating the position in the assembly string
207  void ga_throw_error_msg(pstring expr, size_type pos,
208  const std::string &msg);
209 
210 # define ga_throw_error(expr, pos, msg) \
211  { std::stringstream ss; ss << msg; \
212  ga_throw_error_msg(expr, pos, ss.str()); \
213  GMM_ASSERT1(false, "Error in assembly string" ); \
214  }
215 
216  // Structure for the tensor associated with a tree node
217  struct assembly_tensor {
218  bool is_copied;
219  int sparsity_; // 0: plain, 1: vectorized base, 2: vectorised grad, ...
220  size_type qdim_; // Dimension of the vectorization for sparsity tensors
221  base_tensor t;
222  assembly_tensor *tensor_copied;
223 
224  const base_tensor &org_tensor() const
225  { return is_copied ? tensor_copied->org_tensor() : t; }
226  base_tensor &org_tensor()
227  { return is_copied ? tensor_copied->org_tensor() : t; }
228 
229  const base_tensor &tensor() const
230  { return (is_copied ? tensor_copied->tensor() : t); }
231 
232  base_tensor &tensor()
233  { return (is_copied ? tensor_copied->tensor() : t); }
234 
235  void set_sparsity(int sp, size_type q)
236  { sparsity_ = sp; qdim_ = q; }
237 
238  size_type qdim() { return is_copied ? tensor_copied->qdim() : qdim_; }
239 
240  int sparsity() const
241  { return is_copied ? tensor_copied->sparsity() : sparsity_; }
242 
243  inline void set_to_original() { is_copied = false; }
244  inline void set_to_copy(assembly_tensor &t_) {
245  is_copied = true; sparsity_ = t_.sparsity_; qdim_ = t_.qdim_;
246  t = t_.org_tensor(); tensor_copied = &(t_);
247  }
248 
249  inline void adjust_sizes(const bgeot::multi_index &ssizes)
250  { t.adjust_sizes(ssizes); }
251 
252  inline void adjust_sizes()
253  { if (t.sizes().size() || t.size() != 1) t.init(); }
254 
255  inline void adjust_sizes(size_type i)
256  { if (t.sizes().size() != 1 || t.sizes()[0] != i) t.init(i); }
257 
258  inline void adjust_sizes(size_type i, size_type j) {
259  if (t.sizes().size() != 2 || t.sizes()[0] != i || t.sizes()[1] != j)
260  t.init(i, j);
261  }
262 
263  inline void adjust_sizes(size_type i, size_type j, size_type k) {
264  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
265  || t.sizes()[2] != k)
266  t.init(i, j, k);
267  }
268  inline void adjust_sizes(size_type i, size_type j,
269  size_type k, size_type l) {
270  if (t.sizes().size() != 3 || t.sizes()[0] != i || t.sizes()[1] != j
271  || t.sizes()[2] != k || t.sizes()[3] != l)
272  t.init(i, j, k, l);
273  }
274 
275  void init_scalar_tensor(scalar_type v)
276  { set_to_original(); t.adjust_sizes(); t[0] = v; }
277 
278  void init_vector_tensor(size_type d)
279  { set_to_original(); t.adjust_sizes(d); }
280 
281  void init_matrix_tensor(size_type n, size_type m)
282  { set_to_original(); t.adjust_sizes(n, m); }
283 
284  void init_identity_matrix_tensor(size_type n) {
285  init_matrix_tensor(n, n);
286  auto itw = t.begin();
287  for (size_type i = 0; i < n; ++i)
288  for (size_type j = 0; j < n; ++j)
289  *itw++ = (i == j) ? scalar_type(1) : scalar_type(0);
290  }
291 
292  void init_third_order_tensor(size_type n, size_type m, size_type l)
293  { set_to_original(); t.adjust_sizes(n, m, l); }
294 
295  void init_fourth_order_tensor(size_type n, size_type m,
296  size_type l, size_type k)
297  { set_to_original(); t.adjust_sizes(n, m, l, k); }
298 
299  const bgeot::multi_index &sizes() const { return t.sizes(); }
300 
301  assembly_tensor()
302  : is_copied(false), sparsity_(0), qdim_(0), tensor_copied(0) {}
303  };
304 
305  struct ga_tree_node;
306  typedef ga_tree_node *pga_tree_node;
307 
308 
309  struct ga_tree_node {
310  GA_NODE_TYPE node_type;
311  GA_TOKEN_TYPE op_type;
312  assembly_tensor t;
313  size_type test_function_type; // -1 = undetermined
314  // 0 = no test function,
315  // 1 = first order, 2 = second order,
316  // 3 = both with always first order in first
317  std::string name_test1, name_test2; // variable names corresponding to test
318  // functions when test_function_type > 0.
319  std::string interpolate_name_test1, interpolate_name_test2; // name
320  // of interpolation transformation if any
321  size_type qdim1, qdim2; // Qdims when test_function_type > 0.
322  size_type nbc1, nbc2, nbc3; // For X (nbc1=coordinate number),
323  // macros (nbc1=param number, nbc2,nbc3 type))
324  // and C_MATRIX (nbc1=order).
325  size_type pos; // Position of the first character in string
326  pstring expr; // Original string, for error messages.
327  std::string name; // variable/constant/function/operator name
328  std::string interpolate_name; // For Interpolate : name of transformation
329  std::string interpolate_name_der; // For Interpolate derivative:
330  // name of transformation
331  std::string elementary_name; // For Elementary_transformation :
332  // name of transformation
333  std::string elementary_target;// For Elementary_transformation :
334  // target variable (for its mesh_fem)
335  size_type der1, der2; // For functions and nonlinear operators,
336  // optional derivative or second derivative.
337  bool symmetric_op;
338  pga_tree_node parent; // Parent node
339  std::vector<pga_tree_node> children; // Children nodes
340  scalar_type hash_value; // Hash value to identify nodes.
341  bool marked; // For specific use of some algorithms
342 
343  inline const base_tensor &tensor() const { return t.tensor(); }
344  inline base_tensor &tensor() { return t.tensor(); }
345  int sparsity() const { return t.sparsity(); }
346 
347  inline size_type nb_test_functions() const {
348  if (test_function_type == size_type(-1)) return 0;
349  return test_function_type - (test_function_type >= 2 ? 1 : 0);
350  }
351 
352  inline size_type tensor_order() const
353  { return t.sizes().size() - nb_test_functions(); }
354 
355  inline size_type tensor_test_size() const {
356  size_type st = nb_test_functions();
357  return (st >= 1 ? t.sizes()[0] : 1) * (st == 2 ? t.sizes()[1] : 1);
358  }
359 
360  inline size_type tensor_proper_size() const
361  { return t.org_tensor().size() / tensor_test_size(); }
362 
363  inline size_type tensor_proper_size(size_type i) const
364  { return t.sizes()[nb_test_functions()+i]; }
365 
366 
367  void mult_test(const pga_tree_node n0, const pga_tree_node n1);
368 
369  bool tensor_is_zero() {
370  if (node_type == GA_NODE_ZERO) return true;
371  if (node_type != GA_NODE_CONSTANT) return false;
372  for (size_type i = 0; i < tensor().size(); ++i)
373  if (tensor()[i] != scalar_type(0)) return false;
374  return true;
375  }
376 
377  inline bool is_constant() {
378  return (node_type == GA_NODE_CONSTANT ||
379  (node_type == GA_NODE_ZERO && test_function_type == 0));
380  }
381 
382  inline void init_scalar_tensor(scalar_type v)
383  { t.init_scalar_tensor(v); test_function_type = 0; }
384 
385  inline void init_vector_tensor(size_type d)
386  { t.init_vector_tensor(d); test_function_type = 0; }
387 
388  inline void init_matrix_tensor(size_type n, size_type m)
389  { t.init_matrix_tensor(n, m); test_function_type = 0; }
390 
391  inline void init_identity_matrix_tensor(size_type n)
392  { t.init_identity_matrix_tensor(n); test_function_type = 0; }
393 
394  inline void init_third_order_tensor(size_type n, size_type m, size_type l)
395  { t.init_third_order_tensor(n, m, l); test_function_type = 0; }
396 
397  inline void init_fourth_order_tensor(size_type n, size_type m,
398  size_type l, size_type k)
399  { t.init_fourth_order_tensor(n, m, l, k); test_function_type = 0; }
400 
401  inline void adopt_child(pga_tree_node new_child)
402  { children.push_back(new_child); children.back()->parent = this; }
403 
404  inline void replace_child(pga_tree_node oldchild,
405  pga_tree_node newchild) {
406  bool found = false;
407  for (pga_tree_node &child : children)
408  if (child == oldchild) { child = newchild; found = true; }
409  GMM_ASSERT1(found, "Internal error");
410  }
411 
412  ga_tree_node()
413  : node_type(GA_NODE_VOID), test_function_type(-1), qdim1(0), qdim2(0),
414  nbc1(0), nbc2(0), nbc3(0), pos(0), expr(0), der1(0), der2(0),
415  symmetric_op(false), hash_value(0) {}
416  ga_tree_node(GA_NODE_TYPE ty, size_type p, pstring expr_)
417  : node_type(ty), test_function_type(-1),
418  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
419  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
420  hash_value(0) {}
421  ga_tree_node(scalar_type v, size_type p, pstring expr_)
422  : node_type(GA_NODE_CONSTANT), test_function_type(-1),
423  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
424  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
425  hash_value(0)
426  { init_scalar_tensor(v); }
427  ga_tree_node(const char *n, size_type l, size_type p, pstring expr_)
428  : node_type(GA_NODE_NAME), test_function_type(-1),
429  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
430  pos(p), expr(expr_), name(n, l), der1(0), der2(0), symmetric_op(false),
431  hash_value(0) {}
432  ga_tree_node(GA_TOKEN_TYPE op, size_type p, pstring expr_)
433  : node_type(GA_NODE_OP), op_type(op), test_function_type(-1),
434  qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0),
435  pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false),
436  hash_value(0) {}
437  };
438 
439  struct ga_tree {
440  pga_tree_node root, current_node;
441  std::string secondary_domain;
442 
443  void add_scalar(scalar_type val, size_type pos, pstring expr);
444  void add_allindices(size_type pos, pstring expr);
445  void add_name(const char *name, size_type length, size_type pos,
446  pstring expr);
447  void add_sub_tree(ga_tree &sub_tree);
448  void add_params(size_type pos, pstring expr);
449  void add_matrix(size_type pos, pstring expr);
450  void add_op(GA_TOKEN_TYPE op_type, size_type pos, pstring expr);
451  void clear_node_rec(pga_tree_node pnode);
452  void clear_node(pga_tree_node pnode);
453  void clear() { clear_node_rec(root); root = current_node = nullptr; }
454  void clear_children(pga_tree_node pnode);
455  void replace_node_by_child(pga_tree_node pnode, size_type i);
456  void copy_node(pga_tree_node pnode, pga_tree_node parent,
457  pga_tree_node &child);
458  void duplicate_with_operation(pga_tree_node pnode, GA_TOKEN_TYPE op_type);
459  void duplicate_with_addition(pga_tree_node pnode)
460  { duplicate_with_operation(pnode, GA_PLUS); }
461  void duplicate_with_substraction(pga_tree_node pnode)
462  { duplicate_with_operation(pnode, GA_MINUS); }
463  void insert_node(pga_tree_node pnode, GA_NODE_TYPE node_type);
464  void add_child(pga_tree_node pnode, GA_NODE_TYPE node_type = GA_NODE_VOID);
465  void swap(ga_tree &tree) {
466  std::swap(root, tree.root);
467  std::swap(current_node, tree.current_node);
468  std::swap(secondary_domain, tree.secondary_domain);
469  }
470 
471  ga_tree() : root(nullptr), current_node(nullptr), secondary_domain() {}
472 
473  ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr),
474  secondary_domain(tree.secondary_domain)
475  { if (tree.root) copy_node(tree.root, nullptr, root); }
476 
477  ga_tree &operator =(const ga_tree &tree) {
478  clear(); secondary_domain = tree.secondary_domain;
479  if (tree.root)
480  copy_node(tree.root,nullptr,root);
481  return *this;
482  }
483 
484  ~ga_tree() { clear(); }
485  };
486 
487  // Test equality or equivalence of two sub trees.
488  // version = 0 : strict equality
489  // 1 : give the same result
490  // 2 : give the same result with transposition of test functions
491  bool sub_tree_are_equal
492  (const pga_tree_node pnode1, const pga_tree_node pnode2,
493  const ga_workspace &workspace, int version);
494 
495  // Transform the expression of a node and its sub-nodes in the equivalent
496  // assembly string sent to ostream str
497  void ga_print_node(const pga_tree_node pnode,
498  std::ostream &str);
499  // The same for the whole tree, the result is a std::string
500  std::string ga_tree_to_string(const ga_tree &tree);
501 
502  // Syntax analysis of an assembly string. Conversion to a tree.
503  // No semantic analysis is done. The tree can be inconsistent.
504  void ga_read_string(const std::string &expr, ga_tree &tree,
505  const ga_macro_dictionary &macro_dict);
506  void ga_read_string_reg(const std::string &expr, ga_tree &tree,
507  ga_macro_dictionary &macro_dict);
508 
509 
510 } /* end of namespace */
511 
512 
513 #endif /* GETFEM_GENERIC_ASSEMBLY_TREE_H__ */
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
dal_singleton.h
A simple singleton implementation.
getfem_copyable_ptr.h
A smart pointer that copies the value it points to on copy operations.
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
getfem_models.h
Model representation in Getfem.
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
bgeot_geotrans_inv.h
Inversion of geometric transformations.
gmm_blas.h
Basic linear algebra functions.