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