GetFEM  5.5
gmm_def.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-2026 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
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 gmm_def.h
32  @author Yves Renard <Yves.Renard@insa-lyon.fr>
33  @date October 13, 2002.
34  @brief Basic definitions and tools of GMM.
35 */
36 #ifndef GMM_DEF_H__
37 #define GMM_DEF_H__
38 
39 #include "gmm_ref.h"
40 #include <complex>
41 
42 #ifndef M_PI
43 # define M_E 2.7182818284590452354 /* e */
44 # define M_LOG2E 1.4426950408889634074 /* 1/ln(2) */
45 # define M_LOG10E 0.43429448190325182765 /* 1/ln(10) */
46 # define M_LN2 0.69314718055994530942 /* ln(2) */
47 # define M_LN10 2.30258509299404568402 /* ln(10) */
48 # define M_PI 3.14159265358979323846 /* pi */
49 # define M_PI_2 1.57079632679489661923 /* pi/2 */
50 # define M_PI_4 0.78539816339744830962 /* pi/4 */
51 # define M_1_PI 0.31830988618379067154 /* 1/pi */
52 # define M_2_PI 0.63661977236758134308 /* 2/pi */
53 # define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */
54 # define M_SQRT2 1.41421356237309504880 /* sqrt(2) */
55 # define M_SQRT1_2 0.70710678118654752440 /* sqrt(2)/2 */
56 #endif
57 
58 #ifndef M_PIl
59 # define M_PIl 3.1415926535897932384626433832795029L /* pi */
60 # define M_PI_2l 1.5707963267948966192313216916397514L /* pi/2 */
61 # define M_PI_4l 0.7853981633974483096156608458198757L /* pi/4 */
62 # define M_1_PIl 0.3183098861837906715377675267450287L /* 1/pi */
63 # define M_2_PIl 0.6366197723675813430755350534900574L /* 2/pi */
64 # define M_2_SQRTPIl 1.1283791670955125738961589031215452L /* 2/sqrt(pi) */
65 #endif
66 
67 namespace gmm {
68 
69  typedef size_t size_type;
70 
71  /* ******************************************************************** */
72  /* Specifier types */
73  /* ******************************************************************** */
74  /* not perfectly null, required by aCC 3.33 */
75  struct abstract_null_type {
76  abstract_null_type(int=0) {}
77  template <typename A,typename B,typename C> void operator()(A,B,C) {}
78  }; // specify an information lake.
79 
80  struct linalg_true {};
81  struct linalg_false {};
82 
83  template <typename V, typename W> struct linalg_and
84  { typedef linalg_false bool_type; };
85  template <> struct linalg_and<linalg_true, linalg_true>
86  { typedef linalg_true bool_type; };
87  template <typename V, typename W> struct linalg_or
88  { typedef linalg_true bool_type; };
89  template <> struct linalg_and<linalg_false, linalg_false>
90  { typedef linalg_false bool_type; };
91 
92  struct linalg_const {}; // A reference is either linalg_const,
93  struct linalg_modifiable {}; // linalg_modifiable or linalg_false.
94 
95  struct abstract_vector {}; // The object is a vector
96  struct abstract_matrix {}; // The object is a matrix
97 
98  struct abstract_sparse {}; // sparse matrix or vector
99  struct abstract_skyline {}; // 'sky-line' matrix or vector
100  struct abstract_dense {}; // dense matrix or vector
101  struct abstract_indirect {}; // matrix given by the product with a vector
102 
103  struct row_major {}; // matrix with a row access.
104  struct col_major {}; // matrix with a column access
105  struct row_and_col {}; // both accesses but row preference
106  struct col_and_row {}; // both accesses but column preference
107 
108  template <typename T> struct transposed_type;
109  template<> struct transposed_type<row_major> {typedef col_major t_type;};
110  template<> struct transposed_type<col_major> {typedef row_major t_type;};
111  template<> struct transposed_type<row_and_col> {typedef col_and_row t_type;};
112  template<> struct transposed_type<col_and_row> {typedef row_and_col t_type;};
113 
114  template <typename T> struct principal_orientation_type
115  { typedef abstract_null_type potype; };
116  template<> struct principal_orientation_type<row_major>
117  { typedef row_major potype; };
118  template<> struct principal_orientation_type<col_major>
119  { typedef col_major potype; };
120  template<> struct principal_orientation_type<row_and_col>
121  { typedef row_major potype; };
122  template<> struct principal_orientation_type<col_and_row>
123  { typedef col_major potype; };
124 
125  // template <typename V> struct linalg_traits;
126  template <typename V> struct linalg_traits {
127  typedef abstract_null_type this_type;
128  typedef abstract_null_type linalg_type;
129  typedef abstract_null_type value_type;
130  typedef abstract_null_type is_reference;
131  typedef abstract_null_type& reference;
132  typedef abstract_null_type* iterator;
133  typedef const abstract_null_type* const_iterator;
134  typedef abstract_null_type index_sorted;
135  typedef abstract_null_type storage_type;
136  typedef abstract_null_type origin_type;
137  typedef abstract_null_type const_sub_row_type;
138  typedef abstract_null_type sub_row_type;
139  typedef abstract_null_type const_row_iterator;
140  typedef abstract_null_type row_iterator;
141  typedef abstract_null_type const_sub_col_type;
142  typedef abstract_null_type sub_col_type;
143  typedef abstract_null_type const_col_iterator;
144  typedef abstract_null_type col_iterator;
145  typedef abstract_null_type sub_orientation;
146  };
147 
148  template <typename PT, typename V> struct vect_ref_type;
149  template <typename P, typename V> struct vect_ref_type<P *, V> {
150  typedef typename linalg_traits<V>::reference access_type;
151  typedef typename linalg_traits<V>::iterator iterator;
152  };
153  template <typename P, typename V> struct vect_ref_type<const P *, V> {
154  typedef typename linalg_traits<V>::value_type access_type;
155  typedef typename linalg_traits<V>::const_iterator iterator;
156  };
157 
158  template <typename PT> struct const_pointer;
159  template <typename P> struct const_pointer<P *>
160  { typedef const P* pointer; };
161  template <typename P> struct const_pointer<const P *>
162  { typedef const P* pointer; };
163 
164  template <typename PT> struct modifiable_pointer;
165  template <typename P> struct modifiable_pointer<P *>
166  { typedef P* pointer; };
167  template <typename P> struct modifiable_pointer<const P *>
168  { typedef P* pointer; };
169 
170  template <typename R> struct const_reference;
171  template <typename R> struct const_reference<R &>
172  { typedef const R &reference; };
173  template <typename R> struct const_reference<const R &>
174  { typedef const R &reference; };
175 
176 
177  inline bool is_sparse(abstract_sparse) { return true; }
178  inline bool is_sparse(abstract_dense) { return false; }
179  inline bool is_sparse(abstract_skyline) { return true; }
180  inline bool is_sparse(abstract_indirect) { return false; }
181 
182  template <typename L> inline bool is_sparse(const L &)
183  { return is_sparse(typename linalg_traits<L>::storage_type()); }
184 
185  inline bool is_row_matrix_(row_major) { return true; }
186  inline bool is_row_matrix_(col_major) { return false; }
187  inline bool is_row_matrix_(row_and_col) { return true; }
188  inline bool is_row_matrix_(col_and_row) { return true; }
189 
190  template <typename L> inline bool is_row_matrix(const L &)
191  { return is_row_matrix_(typename linalg_traits<L>::sub_orientation()); }
192 
193  inline bool is_col_matrix_(row_major) { return false; }
194  inline bool is_col_matrix_(col_major) { return true; }
195  inline bool is_col_matrix_(row_and_col) { return true; }
196  inline bool is_col_matrix_(col_and_row) { return true; }
197 
198  template <typename L> inline bool is_col_matrix(const L &)
199  { return is_col_matrix_(typename linalg_traits<L>::sub_orientation()); }
200 
201  inline bool is_col_matrix(row_major) { return false; }
202  inline bool is_col_matrix(col_major) { return true; }
203  inline bool is_row_matrix(row_major) { return true; }
204  inline bool is_row_matrix(col_major) { return false; }
205 
206  template <typename L> inline bool is_const_reference(L) { return false; }
207  inline bool is_const_reference(linalg_const) { return true; }
208 
209 
210  template <typename T> struct is_gmm_interfaced_ {
211  typedef linalg_true result;
212  };
213 
214  template<> struct is_gmm_interfaced_<abstract_null_type> {
215  typedef linalg_false result;
216  };
217 
218  template <typename T> struct is_gmm_interfaced {
219  typedef typename is_gmm_interfaced_<typename gmm::linalg_traits<T>::this_type >::result result;
220  };
221 
222  /* ******************************************************************** */
223  /* Original type from a pointer or a reference. */
224  /* ******************************************************************** */
225 
226  template <typename V> struct org_type { typedef V t; };
227  template <typename V> struct org_type<V *> { typedef V t; };
228  template <typename V> struct org_type<const V *> { typedef V t; };
229  template <typename V> struct org_type<V &> { typedef V t; };
230  template <typename V> struct org_type<const V &> { typedef V t; };
231 
232  /* ******************************************************************** */
233  /* Types to deal with const object representing a modifiable reference */
234  /* ******************************************************************** */
235 
236  template <typename PT, typename R> struct mref_type_
237  { typedef abstract_null_type return_type; };
238  template <typename L, typename R> struct mref_type_<L *, R>
239  { typedef typename org_type<L>::t & return_type; };
240  template <typename L, typename R> struct mref_type_<const L *, R>
241  { typedef const typename org_type<L>::t & return_type; };
242  template <typename L> struct mref_type_<L *, linalg_const>
243  { typedef const typename org_type<L>::t & return_type; };
244  template <typename L> struct mref_type_<const L *, linalg_const>
245  { typedef const typename org_type<L>::t & return_type; };
246  template <typename L> struct mref_type_<const L *, linalg_modifiable>
247  { typedef typename org_type<L>::t & return_type; };
248  template <typename L> struct mref_type_<L *, linalg_modifiable>
249  { typedef typename org_type<L>::t & return_type; };
250 
251  template <typename PT> struct mref_type {
252  typedef typename std::iterator_traits<PT>::value_type L;
253  typedef typename mref_type_<PT,
254  typename linalg_traits<L>::is_reference>::return_type return_type;
255  };
256 
257  template <typename L> typename mref_type<const L *>::return_type
258  linalg_cast(const L &l)
259  { return const_cast<typename mref_type<const L *>::return_type>(l); }
260 
261  template <typename L> typename mref_type<L *>::return_type linalg_cast(L &l)
262  { return const_cast<typename mref_type<L *>::return_type>(l); }
263 
264  template <typename L, typename R> struct cref_type_
265  { typedef abstract_null_type return_type; };
266  template <typename L> struct cref_type_<L, linalg_modifiable>
267  { typedef typename org_type<L>::t & return_type; };
268  template <typename L> struct cref_type {
269  typedef typename cref_type_<L,
270  typename linalg_traits<L>::is_reference>::return_type return_type;
271  };
272 
273  template <typename L> typename cref_type<L>::return_type
274  linalg_const_cast(const L &l)
275  { return const_cast<typename cref_type<L>::return_type>(l); }
276 
277 
278  // To be used to select between a reference or a const refercence for
279  // the return type of a function
280  // select_return<C1, C2, L *> return C1 if L is a const reference,
281  // C2 otherwise.
282  // select_return<C1, C2, const L *> return C2 if L is a modifiable reference
283  // C1 otherwise.
284  template <typename C1, typename C2, typename REF> struct select_return_ {
285  typedef abstract_null_type return_type;
286  };
287  template <typename C1, typename C2, typename L>
288  struct select_return_<C1, C2, const L &> { typedef C1 return_type; };
289  template <typename C1, typename C2, typename L>
290  struct select_return_<C1, C2, L &> { typedef C2 return_type; };
291  template <typename C1, typename C2, typename PT> struct select_return {
292  typedef typename std::iterator_traits<PT>::value_type L;
293  typedef typename select_return_<C1, C2,
294  typename mref_type<PT>::return_type>::return_type return_type;
295  };
296 
297 
298  // To be used to select between a reference or a const refercence inside
299  // a structure or a linagl_traits
300  // select_ref<C1, C2, L *> return C1 if L is a const reference,
301  // C2 otherwise.
302  // select_ref<C1, C2, const L *> return C2 in any case.
303  template <typename C1, typename C2, typename REF>
304  struct select_ref_ { typedef abstract_null_type ref_type; };
305 
306  template <typename C1, typename C2, typename L>
307  struct select_ref_<C1, C2, const L &> { typedef C1 ref_type; };
308 
309  template <typename C1, typename C2, typename L>
310  struct select_ref_<C1, C2, L &> { typedef C2 ref_type; };
311 
312  template <typename C1, typename C2, typename PT>
313  struct select_ref {
314  typedef typename std::iterator_traits<PT>::value_type L;
315  typedef typename select_ref_<C1, C2,
316  typename mref_type<PT>::return_type>::ref_type ref_type;
317  };
318 
319  template <typename C1, typename C2, typename L>
320  struct select_ref<C1, C2, const L *> { typedef C1 ref_type; };
321 
322 
323  template<typename R> struct is_a_reference_
324  { typedef linalg_true reference; };
325  template<> struct is_a_reference_<linalg_false>
326  { typedef linalg_false reference; };
327 
328  template<typename L> struct is_a_reference {
329  typedef typename is_a_reference_<typename linalg_traits<L>::is_reference>
330  ::reference reference;
331  };
332 
333 
334  template <typename L> inline bool is_original_linalg(const L &)
335  { return is_original_linalg(typename is_a_reference<L>::reference()); }
336  inline bool is_original_linalg(linalg_false) { return true; }
337  inline bool is_original_linalg(linalg_true) { return false; }
338 
339 
340  template <typename PT> struct which_reference
341  { typedef abstract_null_type is_reference; };
342  template <typename PT> struct which_reference<PT *>
343  { typedef linalg_modifiable is_reference; };
344  template <typename PT> struct which_reference<const PT *>
345  { typedef linalg_const is_reference; };
346 
347 
348  template <typename C1, typename C2, typename R> struct select_orientation_
349  { typedef abstract_null_type return_type; };
350  template <typename C1, typename C2>
351  struct select_orientation_<C1, C2, row_major>
352  { typedef C1 return_type; };
353  template <typename C1, typename C2>
354  struct select_orientation_<C1, C2, col_major>
355  { typedef C2 return_type; };
356  template <typename C1, typename C2, typename L> struct select_orientation {
357  typedef typename select_orientation_<C1, C2,
358  typename principal_orientation_type<typename
359  linalg_traits<L>::sub_orientation>::potype>::return_type return_type;
360  };
361 
362  /* ******************************************************************** */
363  /* Operations on scalars */
364  /* ******************************************************************** */
365 
366  template <typename T> inline T sqr(T a) { return T(a * a); }
367  template <typename T> inline T abs(T a) { return (a < T(0)) ? T(-a) : a; }
368  template <typename T> inline T abs(std::complex<T> a)
369  { T x = a.real(), y = a.imag(); return T(::sqrt(x*x+y*y)); }
370  template <typename T> inline T abs_sqr(T a) { return T(a*a); }
371  template <typename T> inline T abs_sqr(std::complex<T> a)
372  { return gmm::sqr(a.real()) + gmm::sqr(a.imag()); }
373  template <typename T> inline T pos(T a) { return (a < T(0)) ? T(0) : a; }
374  template <typename T> inline T neg(T a) { return (a < T(0)) ? T(-a) : T(0); }
375  template <typename T> inline T sgn(T a) { return (a < T(0)) ? T(-1) : T(1); }
376  template <typename T> inline T Heaviside(T a) { return (a < T(0)) ? T(0) : T(1); }
377  inline double random() { return double(rand())/(RAND_MAX+0.5); }
378  template <typename T> inline T random(T)
379  { return T(rand()*2.0)/(T(RAND_MAX)+T(1)/T(2)) - T(1); }
380  template <typename T> inline std::complex<T> random(std::complex<T>)
381  { return std::complex<T>(gmm::random(T()), gmm::random(T())); }
382  template <typename T> inline T irandom(T max)
383  { return T(gmm::random() * double(max)); }
384  template <typename T> inline T conj(T a) { return a; }
385  template <typename T> inline std::complex<T> conj(std::complex<T> a)
386  { return std::conj(a); }
387  template <typename T> inline T real(T a) { return a; }
388  template <typename T> inline T real(std::complex<T> a) { return a.real(); }
389  template <typename T> inline T imag(T ) { return T(0); }
390  template <typename T> inline T imag(std::complex<T> a) { return a.imag(); }
391  template <typename T> inline T sqrt(T a) { return T(::sqrt(a)); }
392  template <typename T> inline std::complex<T> sqrt(std::complex<T> a) {
393  T x = a.real(), y = a.imag();
394  if (x == T(0)) {
395  T t = T(::sqrt(gmm::abs(y) / T(2)));
396  return std::complex<T>(t, y < T(0) ? -t : t);
397  }
398  T t = T(::sqrt(T(2) * (gmm::abs(a) + gmm::abs(x)))), u = t / T(2);
399  return x > T(0) ? std::complex<T>(u, y / t)
400  : std::complex<T>(gmm::abs(y) / t, y < T(0) ? -u : u);
401  }
402  using std::swap;
403 
404 
405  template <typename T> struct number_traits {
406  typedef T magnitude_type;
407  };
408 
409  template <typename T> struct number_traits<std::complex<T> > {
410  typedef T magnitude_type;
411  };
412 
413  template <typename T> inline T conj_product(T a, T b) { return a * b; }
414  template <typename T> inline
415  std::complex<T> conj_product(std::complex<T> a, std::complex<T> b)
416  { return std::conj(a) * b; } // to be optimized ?
417 
418  template <typename T> inline bool is_complex(T) { return false; }
419  template <typename T> inline bool is_complex(std::complex<T> )
420  { return true; }
421 
422 # define magnitude_of_linalg(M) typename number_traits<typename \
423  linalg_traits<M>::value_type>::magnitude_type
424 
425  /* ******************************************************************** */
426  /* types promotion */
427  /* ******************************************************************** */
428 
429  /* should be completed for more specific cases <unsigned int, float> etc */
430  template <typename T1, typename T2, bool c>
431  struct strongest_numeric_type_aux {
432  typedef T1 T;
433  };
434  template <typename T1, typename T2>
435  struct strongest_numeric_type_aux<T1,T2,false> {
436  typedef T2 T;
437  };
438 
439  template <typename T1, typename T2>
440  struct strongest_numeric_type {
441  typedef typename
442  strongest_numeric_type_aux<T1,T2,(sizeof(T1)>sizeof(T2))>::T T;
443  };
444  template <typename T1, typename T2>
445  struct strongest_numeric_type<T1,std::complex<T2> > {
446  typedef typename number_traits<T1>::magnitude_type R1;
447  typedef std::complex<typename strongest_numeric_type<R1,T2>::T > T;
448  };
449  template <typename T1, typename T2>
450  struct strongest_numeric_type<std::complex<T1>,T2 > {
451  typedef typename number_traits<T2>::magnitude_type R2;
452  typedef std::complex<typename strongest_numeric_type<T1,R2>::T > T;
453  };
454  template <typename T1, typename T2>
455  struct strongest_numeric_type<std::complex<T1>,std::complex<T2> > {
456  typedef std::complex<typename strongest_numeric_type<T1,T2>::T > T;
457  };
458 
459  template<> struct strongest_numeric_type<int,float> { typedef float T; };
460  template<> struct strongest_numeric_type<float,int> { typedef float T; };
461  template<> struct strongest_numeric_type<long,float> { typedef float T; };
462  template<> struct strongest_numeric_type<float,long> { typedef float T; };
463  template<> struct strongest_numeric_type<long,double> { typedef double T; };
464  template<> struct strongest_numeric_type<double,long> { typedef double T; };
465 
466  template <typename V1, typename V2>
467  struct strongest_value_type {
468  typedef typename
469  strongest_numeric_type<typename linalg_traits<V1>::value_type,
470  typename linalg_traits<V2>::value_type>::T
471  value_type;
472  };
473  template <typename V1, typename V2, typename V3>
474  struct strongest_value_type3 {
475  typedef typename
476  strongest_value_type<V1, typename
477  strongest_value_type<V2,V3>::value_type>::value_type
478  value_type;
479  };
480 
481 
482 
483  /* ******************************************************************** */
484  /* Basic vectors used */
485  /* ******************************************************************** */
486 
487  template<typename T> struct dense_vector_type
488  { typedef std::vector<T> vector_type; };
489 
490  template <typename T> class wsvector;
491  template <typename T> class rsvector;
492  template <typename T> class dsvector;
493  template<typename T> struct sparse_vector_type
494  { typedef wsvector<T> vector_type; };
495 
496  template <typename T> class slvector;
497  template <typename T> class dense_matrix;
498  template <typename VECT> class row_matrix;
499  template <typename VECT> class col_matrix;
500 
501 
502  /* ******************************************************************** */
503  /* Selects a temporary vector type */
504  /* V if V is a valid vector type, */
505  /* wsvector if V is a reference on a sparse vector, */
506  /* std::vector if V is a reference on a dense vector. */
507  /* ******************************************************************** */
508 
509 
510  template <typename R, typename S, typename L, typename V>
511  struct temporary_vector_ {
512  typedef abstract_null_type vector_type;
513  };
514  template <typename V, typename L>
515  struct temporary_vector_<linalg_true, abstract_sparse, L, V>
516  { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; };
517  template <typename V, typename L>
518  struct temporary_vector_<linalg_true, abstract_skyline, L, V>
519  { typedef slvector<typename linalg_traits<V>::value_type> vector_type; };
520  template <typename V, typename L>
521  struct temporary_vector_<linalg_true, abstract_dense, L, V>
522  { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; };
523  template <typename S, typename V>
524  struct temporary_vector_<linalg_false, S, abstract_vector, V>
525  { typedef V vector_type; };
526  template <typename V>
527  struct temporary_vector_<linalg_false, abstract_dense, abstract_matrix, V>
528  { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; };
529  template <typename V>
530  struct temporary_vector_<linalg_false, abstract_sparse, abstract_matrix, V>
531  { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; };
532 
533  template <typename V> struct temporary_vector {
534  typedef typename temporary_vector_<typename is_a_reference<V>::reference,
535  typename linalg_traits<V>::storage_type,
536  typename linalg_traits<V>::linalg_type,
537  V>::vector_type vector_type;
538  };
539 
540  /* ******************************************************************** */
541  /* Selects a temporary matrix type */
542  /* M if M is a valid matrix type, */
543  /* row_matrix<wsvector> if M is a reference on a sparse matrix, */
544  /* dense_matrix if M is a reference on a dense matrix. */
545  /* ******************************************************************** */
546 
547 
548  template <typename R, typename S, typename L, typename V>
549  struct temporary_matrix_ { typedef abstract_null_type matrix_type; };
550  template <typename V, typename L>
551  struct temporary_matrix_<linalg_true, abstract_sparse, L, V> {
552  typedef typename linalg_traits<V>::value_type T;
553  typedef row_matrix<wsvector<T> > matrix_type;
554  };
555  template <typename V, typename L>
556  struct temporary_matrix_<linalg_true, abstract_skyline, L, V> {
557  typedef typename linalg_traits<V>::value_type T;
558  typedef row_matrix<slvector<T> > matrix_type;
559  };
560  template <typename V, typename L>
561  struct temporary_matrix_<linalg_true, abstract_dense, L, V>
562  { typedef dense_matrix<typename linalg_traits<V>::value_type> matrix_type; };
563  template <typename S, typename V>
564  struct temporary_matrix_<linalg_false, S, abstract_matrix, V>
565  { typedef V matrix_type; };
566 
567  template <typename V> struct temporary_matrix {
568  typedef typename temporary_matrix_<typename is_a_reference<V>::reference,
569  typename linalg_traits<V>::storage_type,
570  typename linalg_traits<V>::linalg_type,
571  V>::matrix_type matrix_type;
572  };
573 
574 
575  template <typename S, typename L, typename V>
576  struct temporary_col_matrix_ { typedef abstract_null_type matrix_type; };
577  template <typename V, typename L>
578  struct temporary_col_matrix_<abstract_sparse, L, V> {
579  typedef typename linalg_traits<V>::value_type T;
580  typedef col_matrix<wsvector<T> > matrix_type;
581  };
582  template <typename V, typename L>
583  struct temporary_col_matrix_<abstract_skyline, L, V> {
584  typedef typename linalg_traits<V>::value_type T;
585  typedef col_matrix<slvector<T> > matrix_type;
586  };
587  template <typename V, typename L>
588  struct temporary_col_matrix_<abstract_dense, L, V>
589  { typedef dense_matrix<typename linalg_traits<V>::value_type> matrix_type; };
590 
591  template <typename V> struct temporary_col_matrix {
592  typedef typename temporary_col_matrix_<
593  typename linalg_traits<V>::storage_type,
594  typename linalg_traits<V>::linalg_type,
595  V>::matrix_type matrix_type;
596  };
597 
598 
599 
600 
601  template <typename S, typename L, typename V>
602  struct temporary_row_matrix_ { typedef abstract_null_type matrix_type; };
603  template <typename V, typename L>
604  struct temporary_row_matrix_<abstract_sparse, L, V> {
605  typedef typename linalg_traits<V>::value_type T;
606  typedef row_matrix<wsvector<T> > matrix_type;
607  };
608  template <typename V, typename L>
609  struct temporary_row_matrix_<abstract_skyline, L, V> {
610  typedef typename linalg_traits<V>::value_type T;
611  typedef row_matrix<slvector<T> > matrix_type;
612  };
613  template <typename V, typename L>
614  struct temporary_row_matrix_<abstract_dense, L, V>
615  { typedef dense_matrix<typename linalg_traits<V>::value_type> matrix_type; };
616 
617  template <typename V> struct temporary_row_matrix {
618  typedef typename temporary_row_matrix_<
619  typename linalg_traits<V>::storage_type,
620  typename linalg_traits<V>::linalg_type,
621  V>::matrix_type matrix_type;
622  };
623 
624 
625 
626  /* ******************************************************************** */
627  /* Selects a temporary dense vector type */
628  /* V if V is a valid dense vector type, */
629  /* std::vector if V is a reference or another type of vector */
630  /* ******************************************************************** */
631 
632  template <typename R, typename S, typename V>
633  struct temporary_dense_vector_ { typedef abstract_null_type vector_type; };
634  template <typename S, typename V>
635  struct temporary_dense_vector_<linalg_true, S, V>
636  { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; };
637  template <typename V>
638  struct temporary_dense_vector_<linalg_false, abstract_sparse, V>
639  { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; };
640  template <typename V>
641  struct temporary_dense_vector_<linalg_false, abstract_skyline, V>
642  { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; };
643  template <typename V>
644  struct temporary_dense_vector_<linalg_false, abstract_dense, V>
645  { typedef V vector_type; };
646 
647  template <typename V> struct temporary_dense_vector {
648  typedef typename temporary_dense_vector_<typename
649  is_a_reference<V>::reference,
650  typename linalg_traits<V>::storage_type, V>::vector_type vector_type;
651  };
652 
653  /* ******************************************************************** */
654  /* Selects a temporary sparse vector type */
655  /* V if V is a valid sparse vector type, */
656  /* wsvector if V is a reference or another type of vector */
657  /* ******************************************************************** */
658 
659  template <typename R, typename S, typename V>
660  struct temporary_sparse_vector_ { typedef abstract_null_type vector_type; };
661  template <typename S, typename V>
662  struct temporary_sparse_vector_<linalg_true, S, V>
663  { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; };
664  template <typename V>
665  struct temporary_sparse_vector_<linalg_false, abstract_sparse, V>
666  { typedef V vector_type; };
667  template <typename V>
668  struct temporary_sparse_vector_<linalg_false, abstract_dense, V>
669  { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; };
670  template <typename V>
671  struct temporary_sparse_vector_<linalg_false, abstract_skyline, V>
672  { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; };
673 
674  template <typename V> struct temporary_sparse_vector {
675  typedef typename temporary_sparse_vector_<typename
676  is_a_reference<V>::reference,
677  typename linalg_traits<V>::storage_type, V>::vector_type vector_type;
678  };
679 
680  /* ******************************************************************** */
681  /* Selects a temporary sky-line vector type */
682  /* V if V is a valid sky-line vector type, */
683  /* slvector if V is a reference or another type of vector */
684  /* ******************************************************************** */
685 
686  template <typename R, typename S, typename V>
687  struct temporary_skyline_vector_
688  { typedef abstract_null_type vector_type; };
689  template <typename S, typename V>
690  struct temporary_skyline_vector_<linalg_true, S, V>
691  { typedef slvector<typename linalg_traits<V>::value_type> vector_type; };
692  template <typename V>
693  struct temporary_skyline_vector_<linalg_false, abstract_skyline, V>
694  { typedef V vector_type; };
695  template <typename V>
696  struct temporary_skyline_vector_<linalg_false, abstract_dense, V>
697  { typedef slvector<typename linalg_traits<V>::value_type> vector_type; };
698  template <typename V>
699  struct temporary_skyline_vector_<linalg_false, abstract_sparse, V>
700  { typedef slvector<typename linalg_traits<V>::value_type> vector_type; };
701 
702  template <typename V> struct temporary_skylines_vector {
703  typedef typename temporary_skyline_vector_<typename
704  is_a_reference<V>::reference,
705  typename linalg_traits<V>::storage_type, V>::vector_type vector_type;
706  };
707 
708  /* ********************************************************************* */
709  /* Definition & Comparison of origins. */
710  /* ********************************************************************* */
711 
712  template <typename L>
713  typename select_return<const typename linalg_traits<L>::origin_type *,
714  typename linalg_traits<L>::origin_type *,
715  L *>::return_type
716  linalg_origin(L &l)
717  { return linalg_traits<L>::origin(linalg_cast(l)); }
718 
719  template <typename L>
720  typename select_return<const typename linalg_traits<L>::origin_type *,
721  typename linalg_traits<L>::origin_type *,
722  const L *>::return_type
723  linalg_origin(const L &l)
724  { return linalg_traits<L>::origin(linalg_cast(l)); }
725 
726  template <typename PT1, typename PT2>
727  bool same_porigin(PT1, PT2) { return false; }
728 
729  template <typename PT>
730  bool same_porigin(PT pt1, PT pt2) { return (pt1 == pt2); }
731 
732  template <typename L1, typename L2>
733  bool same_origin(const L1 &l1, const L2 &l2)
734  { return same_porigin(linalg_origin(l1), linalg_origin(l2)); }
735 
736 
737  /* ******************************************************************** */
738  /* Miscellaneous */
739  /* ******************************************************************** */
740 
741  template <typename V> inline size_type vect_size(const V &v)
742  { return linalg_traits<V>::size(v); }
743 
744  template <typename MAT> inline size_type mat_nrows(const MAT &m)
745  { return linalg_traits<MAT>::nrows(m); }
746 
747  template <typename MAT> inline size_type mat_ncols(const MAT &m)
748  { return linalg_traits<MAT>::ncols(m); }
749 
750 
751  template <typename V> inline
752  typename select_return<typename linalg_traits<V>::const_iterator,
753  typename linalg_traits<V>::iterator,
754  V *>::return_type
755  vect_begin(V &v)
756  { return linalg_traits<V>::begin(linalg_cast(v)); }
757 
758  template <typename V> inline
759  typename select_return<typename linalg_traits<V>::const_iterator,
760  typename linalg_traits<V>::iterator,
761  const V *>::return_type
762  vect_begin(const V &v)
763  { return linalg_traits<V>::begin(linalg_cast(v)); }
764 
765  template <typename V> inline
766  typename linalg_traits<V>::const_iterator
767  vect_const_begin(const V &v)
768  { return linalg_traits<V>::begin(v); }
769 
770  template <typename V> inline
771  typename select_return<typename linalg_traits<V>::const_iterator,
772  typename linalg_traits<V>::iterator,
773  V *>::return_type
774  vect_end(V &v)
775  { return linalg_traits<V>::end(linalg_cast(v)); }
776 
777  template <typename V> inline
778  typename select_return<typename linalg_traits<V>::const_iterator,
779  typename linalg_traits<V>::iterator,
780  const V *>::return_type
781  vect_end(const V &v)
782  { return linalg_traits<V>::end(linalg_cast(v)); }
783 
784  template <typename V> inline
785  typename linalg_traits<V>::const_iterator
786  vect_const_end(const V &v)
787  { return linalg_traits<V>::end(v); }
788 
789  template <typename M> inline
790  typename select_return<typename linalg_traits<M>::const_row_iterator,
791  typename linalg_traits<M>::row_iterator,
792  M *>::return_type
793  mat_row_begin(M &m) { return linalg_traits<M>::row_begin(linalg_cast(m)); }
794 
795  template <typename M> inline
796  typename select_return<typename linalg_traits<M>::const_row_iterator,
797  typename linalg_traits<M>::row_iterator,
798  const M *>::return_type
799  mat_row_begin(const M &m)
800  { return linalg_traits<M>::row_begin(linalg_cast(m)); }
801 
802  template <typename M> inline typename linalg_traits<M>::const_row_iterator
803  mat_row_const_begin(const M &m)
804  { return linalg_traits<M>::row_begin(m); }
805 
806  template <typename M> inline
807  typename select_return<typename linalg_traits<M>::const_row_iterator,
808  typename linalg_traits<M>::row_iterator,
809  M *>::return_type
810  mat_row_end(M &v) {
811  return linalg_traits<M>::row_end(linalg_cast(v));
812  }
813 
814  template <typename M> inline
815  typename select_return<typename linalg_traits<M>::const_row_iterator,
816  typename linalg_traits<M>::row_iterator,
817  const M *>::return_type
818  mat_row_end(const M &v) {
819  return linalg_traits<M>::row_end(linalg_cast(v));
820  }
821 
822  template <typename M> inline
823  typename linalg_traits<M>::const_row_iterator
824  mat_row_const_end(const M &v)
825  { return linalg_traits<M>::row_end(v); }
826 
827  template <typename M> inline
828  typename select_return<typename linalg_traits<M>::const_col_iterator,
829  typename linalg_traits<M>::col_iterator,
830  M *>::return_type
831  mat_col_begin(M &v) {
832  return linalg_traits<M>::col_begin(linalg_cast(v));
833  }
834 
835  template <typename M> inline
836  typename select_return<typename linalg_traits<M>::const_col_iterator,
837  typename linalg_traits<M>::col_iterator,
838  const M *>::return_type
839  mat_col_begin(const M &v) {
840  return linalg_traits<M>::col_begin(linalg_cast(v));
841  }
842 
843  template <typename M> inline
844  typename linalg_traits<M>::const_col_iterator
845  mat_col_const_begin(const M &v)
846  { return linalg_traits<M>::col_begin(v); }
847 
848  template <typename M> inline
849  typename linalg_traits<M>::const_col_iterator
850  mat_col_const_end(const M &v)
851  { return linalg_traits<M>::col_end(v); }
852 
853  template <typename M> inline
854  typename select_return<typename linalg_traits<M>::const_col_iterator,
855  typename linalg_traits<M>::col_iterator,
856  M *>::return_type
857  mat_col_end(M &m)
858  { return linalg_traits<M>::col_end(linalg_cast(m)); }
859 
860  template <typename M> inline
861  typename select_return<typename linalg_traits<M>::const_col_iterator,
862  typename linalg_traits<M>::col_iterator,
863  const M *>::return_type
864  mat_col_end(const M &m)
865  { return linalg_traits<M>::col_end(linalg_cast(m)); }
866 
867  template <typename MAT> inline
868  typename select_return<typename linalg_traits<MAT>::const_sub_row_type,
869  typename linalg_traits<MAT>::sub_row_type,
870  const MAT *>::return_type
871  mat_row(const MAT &m, size_type i)
872  { return linalg_traits<MAT>::row(mat_row_begin(m) + i); }
873 
874  template <typename MAT> inline
875  typename select_return<typename linalg_traits<MAT>::const_sub_row_type,
876  typename linalg_traits<MAT>::sub_row_type,
877  MAT *>::return_type
878  mat_row(MAT &m, size_type i)
879  { return linalg_traits<MAT>::row(mat_row_begin(m) + i); }
880 
881  template <typename MAT> inline
882  typename linalg_traits<MAT>::const_sub_row_type
883  mat_const_row(const MAT &m, size_type i)
884  { return linalg_traits<MAT>::row(mat_row_const_begin(m) + i); }
885 
886  template <typename MAT> inline
887  typename select_return<typename linalg_traits<MAT>::const_sub_col_type,
888  typename linalg_traits<MAT>::sub_col_type,
889  const MAT *>::return_type
890  mat_col(const MAT &m, size_type i)
891  { return linalg_traits<MAT>::col(mat_col_begin(m) + i); }
892 
893 
894  template <typename MAT> inline
895  typename select_return<typename linalg_traits<MAT>::const_sub_col_type,
896  typename linalg_traits<MAT>::sub_col_type,
897  MAT *>::return_type
898  mat_col(MAT &m, size_type i)
899  { return linalg_traits<MAT>::col(mat_col_begin(m) + i); }
900 
901  template <typename MAT> inline
902  typename linalg_traits<MAT>::const_sub_col_type
903  mat_const_col(const MAT &m, size_type i)
904  { return linalg_traits<MAT>::col(mat_col_const_begin(m) + i); }
905 
906  /* ********************************************************************* */
907  /* Set to begin end set to end for iterators on non-const sparse vectors.*/
908  /* ********************************************************************* */
909 
910  template <typename IT, typename ORG, typename VECT> inline
911  void set_to_begin(IT &it, ORG o, VECT *, linalg_false)
912  { it = vect_begin(*o); }
913 
914  template <typename IT, typename ORG, typename VECT> inline
915  void set_to_begin(IT &it, ORG o, const VECT *, linalg_false)
916  { it = vect_const_begin(*o); }
917 
918  template <typename IT, typename ORG, typename VECT> inline
919  void set_to_end(IT &it, ORG o, VECT *, linalg_false)
920  { it = vect_end(*o); }
921 
922  template <typename IT, typename ORG, typename VECT> inline
923  void set_to_end(IT &it, ORG o, const VECT *, linalg_false)
924  { it = vect_const_end(*o); }
925 
926 
927  template <typename IT, typename ORG, typename VECT> inline
928  void set_to_begin(IT &, ORG, VECT *, linalg_const) { }
929 
930  template <typename IT, typename ORG, typename VECT> inline
931  void set_to_begin(IT &, ORG, const VECT *, linalg_const) { }
932 
933  template <typename IT, typename ORG, typename VECT> inline
934  void set_to_end(IT &, ORG, VECT *, linalg_const) { }
935 
936  template <typename IT, typename ORG, typename VECT> inline
937  void set_to_end(IT &, ORG, const VECT *, linalg_const) { }
938 
939  template <typename IT, typename ORG, typename VECT> inline
940  void set_to_begin(IT &, ORG, VECT *v, linalg_modifiable)
941  { GMM_ASSERT3(!is_sparse(*v), "internal_error"); (void)v; }
942 
943  template <typename IT, typename ORG, typename VECT> inline
944  void set_to_begin(IT &, ORG, const VECT *v, linalg_modifiable)
945  { GMM_ASSERT3(!is_sparse(*v), "internal_error"); (void)v; }
946 
947  template <typename IT, typename ORG, typename VECT> inline
948  void set_to_end(IT &, ORG, VECT *v, linalg_modifiable)
949  { GMM_ASSERT3(!is_sparse(*v), "internal_error"); (void)v; }
950 
951  template <typename IT, typename ORG, typename VECT> inline
952  void set_to_end(IT &, ORG, const VECT *v, linalg_modifiable)
953  { GMM_ASSERT3(!is_sparse(*v), "internal_error"); (void)v; }
954 
955  /* ******************************************************************** */
956  /* General index for certain algorithms. */
957  /* ******************************************************************** */
958 
959  template<class IT>
960  size_type index_of_it(const IT &it, size_type, abstract_sparse)
961  { return it.index(); }
962  template<class IT>
963  size_type index_of_it(const IT &it, size_type, abstract_skyline)
964  { return it.index(); }
965  template<class IT>
966  size_type index_of_it(const IT &, size_type k, abstract_dense)
967  { return k; }
968 
969  /* ********************************************************************* */
970  /* Numeric limits. */
971  /* ********************************************************************* */
972 
973  template<typename T> inline T default_tol(T) {
974  using namespace std;
975  static T tol(10);
976  if (tol == T(10)) {
977  if (numeric_limits<T>::is_specialized)
978  tol = numeric_limits<T>::epsilon();
979  else {
980  int i=int(sizeof(T)/4); while(i-- > 0) tol*=T(1E-8);
981  GMM_WARNING1("The numeric type " << typeid(T).name()
982  << " has no numeric_limits defined !!\n"
983  << "Taking " << tol << " as default tolerance");
984  }
985  }
986  return tol;
987  }
988  template<typename T> inline T default_tol(std::complex<T>)
989  { return default_tol(T()); }
990 
991  template<typename T> inline T default_min(T) {
992  using namespace std;
993  static T mi(10);
994  if (mi == T(10)) {
995  if (numeric_limits<T>::is_specialized)
996  mi = std::numeric_limits<T>::min();
997  else {
998  mi = T(0);
999  GMM_WARNING1("The numeric type " << typeid(T).name()
1000  << " has no numeric_limits defined !!\n"
1001  << "Taking 0 as default minimum");
1002  }
1003  }
1004  return mi;
1005  }
1006  template<typename T> inline T default_min(std::complex<T>)
1007  { return default_min(T()); }
1008 
1009  template<typename T> inline T default_max(T) {
1010  using namespace std;
1011  static T mi(10);
1012  if (mi == T(10)) {
1013  if (numeric_limits<T>::is_specialized)
1014  mi = std::numeric_limits<T>::max();
1015  else {
1016  mi = T(1);
1017  GMM_WARNING1("The numeric type " << typeid(T).name()
1018  << " has no numeric_limits defined !!\n"
1019  << "Taking 1 as default maximum !");
1020  }
1021  }
1022  return mi;
1023  }
1024  template<typename T> inline T default_max(std::complex<T>)
1025  { return default_max(T()); }
1026 
1027 
1028  /*
1029  use safe_divide to avoid NaNs when dividing very small complex
1030  numbers, for example
1031  std::complex<float>(1e-23,1e-30)/std::complex<float>(1e-23,1e-30)
1032  */
1033  template<typename T> inline T safe_divide(T a, T b) { return a/b; }
1034  template<typename T> inline std::complex<T>
1035  safe_divide(std::complex<T> a, std::complex<T> b) {
1036  T m = std::max(gmm::abs(b.real()), gmm::abs(b.imag()));
1037  a = std::complex<T>(a.real()/m, a.imag()/m);
1038  b = std::complex<T>(b.real()/m, b.imag()/m);
1039  return a / b;
1040  }
1041 
1042 
1043  /* ******************************************************************** */
1044  /* Write */
1045  /* ******************************************************************** */
1046 
1047  template <typename T> struct cast_char_type { typedef T return_type; };
1048  template <> struct cast_char_type<signed char> { typedef int return_type; };
1049  template <> struct cast_char_type<unsigned char>
1050  { typedef unsigned int return_type; };
1051  template <typename T> inline typename cast_char_type<T>::return_type
1052  cast_char(const T &c) { return typename cast_char_type<T>::return_type(c); }
1053 
1054 
1055  template <typename L> inline void write(std::ostream &o, const L &l)
1056  { write(o, l, typename linalg_traits<L>::linalg_type()); }
1057 
1058  template <typename L> void write(std::ostream &o, const L &l,
1059  abstract_vector) {
1060  o << "vector(" << vect_size(l) << ") [";
1061  write(o, l, typename linalg_traits<L>::storage_type());
1062  o << " ]";
1063  }
1064 
1065  template <typename L> void write(std::ostream &o, const L &l,
1066  abstract_sparse) {
1067  typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
1068  ite = vect_const_end(l);
1069  for (; it != ite; ++it)
1070  o << " (r" << it.index() << ", " << cast_char(*it) << ")";
1071  }
1072 
1073  template <typename L> void write(std::ostream &o, const L &l,
1074  abstract_dense) {
1075  typename linalg_traits<L>::const_iterator it = vect_const_begin(l),
1076  ite = vect_const_end(l);
1077  if (it != ite) o << " " << cast_char(*it++);
1078  for (; it != ite; ++it) o << ", " << cast_char(*it);
1079  }
1080 
1081  template <typename L> void write(std::ostream &o, const L &l,
1082  abstract_skyline) {
1083  typedef typename linalg_traits<L>::const_iterator const_iterator;
1084  const_iterator it = vect_const_begin(l), ite = vect_const_end(l);
1085  if (it != ite) {
1086  o << "<r+" << it.index() << ">";
1087  if (it != ite) o << " " << cast_char(*it++);
1088  for (; it != ite; ++it) { o << ", " << cast_char(*it); }
1089  }
1090  }
1091 
1092  template <typename L> inline void write(std::ostream &o, const L &l,
1093  abstract_matrix) {
1094  write(o, l, typename linalg_traits<L>::sub_orientation());
1095  }
1096 
1097 
1098  template <typename L> void write(std::ostream &o, const L &l,
1099  row_major) {
1100  o << "matrix(" << mat_nrows(l) << ", " << mat_ncols(l) << ")" << endl;
1101  for (size_type i = 0; i < mat_nrows(l); ++i) {
1102  o << "(";
1103  write(o, mat_const_row(l, i), typename linalg_traits<L>::storage_type());
1104  o << " )\n";
1105  }
1106  }
1107 
1108  template <typename L> inline
1109  void write(std::ostream &o, const L &l, row_and_col)
1110  { write(o, l, row_major()); }
1111 
1112  template <typename L> inline
1113  void write(std::ostream &o, const L &l, col_and_row)
1114  { write(o, l, row_major()); }
1115 
1116  template <typename L> void write(std::ostream &o, const L &l, col_major) {
1117  o << "matrix(" << mat_nrows(l) << ", " << mat_ncols(l) << ")" << endl;
1118  for (size_type i = 0; i < mat_nrows(l); ++i) {
1119  o << "(";
1120  if (is_sparse(l)) { // not optimized ...
1121  for (size_type j = 0; j < mat_ncols(l); ++j)
1122  if (l(i,j) != typename linalg_traits<L>::value_type(0))
1123  o << " (r" << j << ", " << l(i,j) << ")";
1124  }
1125  else {
1126  if (mat_ncols(l) != 0) o << ' ' << l(i, 0);
1127  for (size_type j = 1; j < mat_ncols(l); ++j) o << ", " << l(i, j);
1128  }
1129  o << " )\n";
1130  }
1131  }
1132 
1133 }
1134 
1135 #endif // GMM_DEF_H__
gmm::rsvector
sparse vector built upon std::vector.
Definition: gmm_vector.h:995
gmm_ref.h
Provide some simple pseudo-containers.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48