1 // -*- c++ -*- (enables emacs c++ mode) 2 //=========================================================================== 3 // 4 // Copyright (C) 2002-2008 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 2.1 of the License, or 11 // (at your option) any later version. 12 // This program is distributed in the hope that it will be useful, but 13 // WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14 // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15 // License for more details. 16 // You should have received a copy of the GNU Lesser General Public License 17 // along with this program; if not, write to the Free Software Foundation, 18 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. 19 // 20 // As a special exception, you may use this file as part of a free software 21 // library without restriction. Specifically, if other files instantiate 22 // templates or use macros or inline functions from this file, or you compile 23 // this file and link it with other files to produce an executable, this 24 // file does not by itself cause the resulting executable to be covered by 25 // the GNU General Public License. This exception does not however 26 // invalidate any other reasons why the executable file might be covered by 27 // the GNU 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) {} operatorabstract_null_type77 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 /* types to deal with const object representing a modifiable reference */ 224 /* ******************************************************************** */ 225 226 template <typename PT, typename R> struct mref_type_ 227 { typedef abstract_null_type return_type; }; 228 template <typename L, typename R> struct mref_type_<L *, R> 229 { typedef L & return_type; }; 230 template <typename L, typename R> struct mref_type_<const L *, R> 231 { typedef const L & return_type; }; 232 template <typename L> struct mref_type_<L *, linalg_const> 233 { typedef const L & return_type; }; 234 template <typename L> struct mref_type_<const L *, linalg_const> 235 { typedef const L & return_type; }; 236 template <typename L> struct mref_type_<const L *, linalg_modifiable> 237 { typedef L & return_type; }; 238 template <typename L> struct mref_type_<L *, linalg_modifiable> 239 { typedef L & return_type; }; 240 241 template <typename PT> struct mref_type { 242 typedef typename std::iterator_traits<PT>::value_type L; 243 typedef typename mref_type_<PT, 244 typename linalg_traits<L>::is_reference>::return_type return_type; 245 }; 246 247 template <typename L> typename mref_type<const L *>::return_type 248 linalg_cast(const L &l) 249 { return const_cast<typename mref_type<const L *>::return_type>(l); } 250 251 template <typename L> typename mref_type<L *>::return_type linalg_cast(L &l) 252 { return const_cast<typename mref_type<L *>::return_type>(l); } 253 254 template <typename L, typename R> struct cref_type_ 255 { typedef abstract_null_type return_type; }; 256 template <typename L> struct cref_type_<L, linalg_modifiable> 257 { typedef L & return_type; }; 258 template <typename L> struct cref_type { 259 typedef typename cref_type_<L, 260 typename linalg_traits<L>::is_reference>::return_type return_type; 261 }; 262 263 template <typename L> typename cref_type<L>::return_type 264 linalg_const_cast(const L &l) 265 { return const_cast<typename cref_type<L>::return_type>(l); } 266 267 268 // To be used to select between a reference or a const refercence for 269 // the return type of a function 270 // select_return<C1, C2, L *> return C1 if L is a const reference, 271 // C2 otherwise. 272 // select_return<C1, C2, const L *> return C2 if L is a modifiable reference 273 // C1 otherwise. 274 template <typename C1, typename C2, typename REF> struct select_return_ { 275 typedef abstract_null_type return_type; 276 }; 277 template <typename C1, typename C2, typename L> 278 struct select_return_<C1, C2, const L &> { typedef C1 return_type; }; 279 template <typename C1, typename C2, typename L> 280 struct select_return_<C1, C2, L &> { typedef C2 return_type; }; 281 template <typename C1, typename C2, typename PT> struct select_return { 282 typedef typename std::iterator_traits<PT>::value_type L; 283 typedef typename select_return_<C1, C2, 284 typename mref_type<PT>::return_type>::return_type return_type; 285 }; 286 287 288 // To be used to select between a reference or a const refercence inside 289 // a structure or a linagl_traits 290 // select_ref<C1, C2, L *> return C1 if L is a const reference, 291 // C2 otherwise. 292 // select_ref<C1, C2, const L *> return C2 in any case. 293 template <typename C1, typename C2, typename REF> struct select_ref_ 294 { typedef abstract_null_type ref_type; }; 295 template <typename C1, typename C2, typename L> 296 struct select_ref_<C1, C2, const L &> { typedef C1 ref_type; }; 297 template <typename C1, typename C2, typename L> 298 struct select_ref_<C1, C2, L &> { typedef C2 ref_type; }; 299 template <typename C1, typename C2, typename PT> struct select_ref { 300 typedef typename std::iterator_traits<PT>::value_type L; 301 typedef typename select_ref_<C1, C2, 302 typename mref_type<PT>::return_type>::ref_type ref_type; 303 }; 304 template <typename C1, typename C2, typename L> 305 struct select_ref<C1, C2, const L *> 306 { typedef C1 ref_type; }; 307 308 309 template<typename R> struct is_a_reference_ 310 { typedef linalg_true reference; }; 311 template<> struct is_a_reference_<linalg_false> 312 { typedef linalg_false reference; }; 313 314 template<typename L> struct is_a_reference { 315 typedef typename is_a_reference_<typename linalg_traits<L>::is_reference> 316 ::reference reference; 317 }; 318 319 320 template <typename L> inline bool is_original_linalg(const L &) 321 { return is_original_linalg(typename is_a_reference<L>::reference()); } 322 inline bool is_original_linalg(linalg_false) { return true; } 323 inline bool is_original_linalg(linalg_true) { return false; } 324 325 326 template <typename PT> struct which_reference 327 { typedef abstract_null_type is_reference; }; 328 template <typename PT> struct which_reference<PT *> 329 { typedef linalg_modifiable is_reference; }; 330 template <typename PT> struct which_reference<const PT *> 331 { typedef linalg_const is_reference; }; 332 333 334 template <typename C1, typename C2, typename R> struct select_orientation_ 335 { typedef abstract_null_type return_type; }; 336 template <typename C1, typename C2> 337 struct select_orientation_<C1, C2, row_major> 338 { typedef C1 return_type; }; 339 template <typename C1, typename C2> 340 struct select_orientation_<C1, C2, col_major> 341 { typedef C2 return_type; }; 342 template <typename C1, typename C2, typename L> struct select_orientation { 343 typedef typename select_orientation_<C1, C2, 344 typename principal_orientation_type<typename 345 linalg_traits<L>::sub_orientation>::potype>::return_type return_type; 346 }; 347 348 /* ******************************************************************** */ 349 /* Operations on scalars */ 350 /* ******************************************************************** */ 351 352 template <typename T> inline T sqr(T a) { return a * a; } 353 template <typename T> inline T abs(T a) { return (a < T(0)) ? T(-a) : a; } 354 template <typename T> inline T abs(std::complex<T> a) 355 { T x = a.real(), y = a.imag(); return ::sqrt(x*x+y*y); } 356 template <typename T> inline T abs_sqr(T a) { return a*a; } 357 template <typename T> inline T abs_sqr(std::complex<T> a) 358 { return gmm::sqr(a.real()) + gmm::sqr(a.imag()); } 359 template <typename T> inline T pos(T a) { return (a < T(0)) ? T(0) : a; } 360 template <typename T> inline T neg(T a) { return (a < T(0)) ? T(-a) : T(0); } 361 template <typename T> inline T sgn(T a) { return (a < T(0)) ? T(-1) : T(1); } 362 inline double random() { return double(rand())/(RAND_MAX+0.5); } 363 template <typename T> inline T random(T) 364 { return T(rand()*2.0)/(T(RAND_MAX)+T(1)/T(2)) - T(1); } 365 template <typename T> inline std::complex<T> random(std::complex<T>) 366 { return std::complex<T>(gmm::random(T()), gmm::random(T())); } 367 template <typename T> inline T irandom(T max) 368 { return T(gmm::random() * max); } 369 template <typename T> inline T conj(T a) { return a; } 370 template <typename T> inline std::complex<T> conj(std::complex<T> a) 371 { return std::conj(a); } 372 template <typename T> inline T real(T a) { return a; } 373 template <typename T> inline T real(std::complex<T> a) { return a.real(); } 374 template <typename T> inline T imag(T ) { return T(0); } 375 template <typename T> inline T imag(std::complex<T> a) { return a.imag(); } 376 template <typename T> inline T sqrt(T a) { return ::sqrt(a); } 377 template <typename T> inline std::complex<T> sqrt(std::complex<T> a) { 378 T x = a.real(), y = a.imag(); 379 if (x == T(0)) { 380 T t = ::sqrt(gmm::abs(y) / T(2)); 381 return std::complex<T>(t, y < T(0) ? -t : t); 382 } 383 T t = ::sqrt(T(2) * (gmm::abs(a) + gmm::abs(x))), u = t / T(2); 384 return x > T(0) ? std::complex<T>(u, y / t) 385 : std::complex<T>(gmm::abs(y) / t, y < T(0) ? -u : u); 386 } 387 using std::swap; 388 389 390 template <typename T> struct number_traits { 391 typedef T magnitude_type; 392 }; 393 394 template <typename T> struct number_traits<std::complex<T> > { 395 typedef T magnitude_type; 396 }; 397 398 template <typename T> inline T conj_product(T a, T b) { return a * b; } 399 template <typename T> inline 400 std::complex<T> conj_product(std::complex<T> a, std::complex<T> b) 401 { return std::conj(a) * b; } // to be optimized ? 402 403 template <typename T> inline bool is_complex(T) { return false; } 404 template <typename T> inline bool is_complex(std::complex<T> ) 405 { return true; } 406 407 # define magnitude_of_linalg(M) typename number_traits<typename \ 408 linalg_traits<M>::value_type>::magnitude_type 409 410 template<typename T> inline std::complex<T> operator*(const std::complex<T>& a, int b) { 411 return a*T(b); 412 } 413 template<typename T> inline std::complex<T> operator*(int b, const std::complex<T>& a) { 414 return a*T(b); 415 } 416 417 /* ******************************************************************** */ 418 /* types promotion */ 419 /* ******************************************************************** */ 420 421 /* should be completed for more specific cases <unsigned int, float> etc */ 422 template <typename T1, typename T2, bool c> 423 struct strongest_numeric_type_aux { 424 typedef T1 T; 425 }; 426 template <typename T1, typename T2> 427 struct strongest_numeric_type_aux<T1,T2,false> { 428 typedef T2 T; 429 }; 430 431 template <typename T1, typename T2> 432 struct strongest_numeric_type { 433 typedef typename 434 strongest_numeric_type_aux<T1,T2,(sizeof(T1)>sizeof(T2))>::T T; 435 }; 436 template <typename T1, typename T2> 437 struct strongest_numeric_type<T1,std::complex<T2> > { 438 typedef typename number_traits<T1>::magnitude_type R1; 439 typedef std::complex<typename strongest_numeric_type<R1,T2>::T > T; 440 }; 441 template <typename T1, typename T2> 442 struct strongest_numeric_type<std::complex<T1>,T2 > { 443 typedef typename number_traits<T2>::magnitude_type R2; 444 typedef std::complex<typename strongest_numeric_type<T1,R2>::T > T; 445 }; 446 template <typename T1, typename T2> 447 struct strongest_numeric_type<std::complex<T1>,std::complex<T2> > { 448 typedef std::complex<typename strongest_numeric_type<T1,T2>::T > T; 449 }; 450 451 template<> struct strongest_numeric_type<int,float> { typedef float T; }; 452 template<> struct strongest_numeric_type<float,int> { typedef float T; }; 453 template<> struct strongest_numeric_type<long,float> { typedef float T; }; 454 template<> struct strongest_numeric_type<float,long> { typedef float T; }; 455 template<> struct strongest_numeric_type<long,double> { typedef double T; }; 456 template<> struct strongest_numeric_type<double,long> { typedef double T; }; 457 458 template <typename V1, typename V2> 459 struct strongest_value_type { 460 typedef typename 461 strongest_numeric_type<typename linalg_traits<V1>::value_type, 462 typename linalg_traits<V2>::value_type>::T 463 value_type; 464 }; 465 template <typename V1, typename V2, typename V3> 466 struct strongest_value_type3 { 467 typedef typename 468 strongest_value_type<V1, typename 469 strongest_value_type<V2,V3>::value_type>::value_type 470 value_type; 471 }; 472 473 474 475 /* ******************************************************************** */ 476 /* Basic vectors used */ 477 /* ******************************************************************** */ 478 479 template<typename T> struct dense_vector_type 480 { typedef std::vector<T> vector_type; }; 481 482 template <typename T> class wsvector; 483 template <typename T> class rsvector; 484 template<typename T> struct sparse_vector_type 485 { typedef wsvector<T> vector_type; }; 486 487 template <typename T> class slvector; 488 template <typename T> class dense_matrix; 489 template <typename VECT> class row_matrix; 490 template <typename VECT> class col_matrix; 491 492 493 /* ******************************************************************** */ 494 /* Selects a temporary vector type */ 495 /* V if V is a valid vector type, */ 496 /* wsvector if V is a reference on a sparse vector, */ 497 /* std::vector if V is a reference on a dense vector. */ 498 /* ******************************************************************** */ 499 500 501 template <typename R, typename S, typename L, typename V> 502 struct temporary_vector_ { 503 typedef abstract_null_type vector_type; 504 }; 505 template <typename V, typename L> 506 struct temporary_vector_<linalg_true, abstract_sparse, L, V> 507 { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; }; 508 template <typename V, typename L> 509 struct temporary_vector_<linalg_true, abstract_skyline, L, V> 510 { typedef slvector<typename linalg_traits<V>::value_type> vector_type; }; 511 template <typename V, typename L> 512 struct temporary_vector_<linalg_true, abstract_dense, L, V> 513 { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; }; 514 template <typename S, typename V> 515 struct temporary_vector_<linalg_false, S, abstract_vector, V> 516 { typedef V vector_type; }; 517 template <typename V> 518 struct temporary_vector_<linalg_false, abstract_dense, abstract_matrix, V> 519 { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; }; 520 template <typename V> 521 struct temporary_vector_<linalg_false, abstract_sparse, abstract_matrix, V> 522 { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; }; 523 524 template <typename V> struct temporary_vector { 525 typedef typename temporary_vector_<typename is_a_reference<V>::reference, 526 typename linalg_traits<V>::storage_type, 527 typename linalg_traits<V>::linalg_type, 528 V>::vector_type vector_type; 529 }; 530 531 /* ******************************************************************** */ 532 /* Selects a temporary matrix type */ 533 /* M if M is a valid matrix type, */ 534 /* row_matrix<wsvector> if M is a reference on a sparse matrix, */ 535 /* dense_matrix if M is a reference on a dense matrix. */ 536 /* ******************************************************************** */ 537 538 539 template <typename R, typename S, typename L, typename V> 540 struct temporary_matrix_ { typedef abstract_null_type matrix_type; }; 541 template <typename V, typename L> 542 struct temporary_matrix_<linalg_true, abstract_sparse, L, V> { 543 typedef typename linalg_traits<V>::value_type T; 544 typedef row_matrix<wsvector<T> > matrix_type; 545 }; 546 template <typename V, typename L> 547 struct temporary_matrix_<linalg_true, abstract_skyline, L, V> { 548 typedef typename linalg_traits<V>::value_type T; 549 typedef row_matrix<slvector<T> > matrix_type; 550 }; 551 template <typename V, typename L> 552 struct temporary_matrix_<linalg_true, abstract_dense, L, V> 553 { typedef dense_matrix<typename linalg_traits<V>::value_type> matrix_type; }; 554 template <typename S, typename V> 555 struct temporary_matrix_<linalg_false, S, abstract_matrix, V> 556 { typedef V matrix_type; }; 557 558 template <typename V> struct temporary_matrix { 559 typedef typename temporary_matrix_<typename is_a_reference<V>::reference, 560 typename linalg_traits<V>::storage_type, 561 typename linalg_traits<V>::linalg_type, 562 V>::matrix_type matrix_type; 563 }; 564 565 566 template <typename S, typename L, typename V> 567 struct temporary_col_matrix_ { typedef abstract_null_type matrix_type; }; 568 template <typename V, typename L> 569 struct temporary_col_matrix_<abstract_sparse, L, V> { 570 typedef typename linalg_traits<V>::value_type T; 571 typedef col_matrix<wsvector<T> > matrix_type; 572 }; 573 template <typename V, typename L> 574 struct temporary_col_matrix_<abstract_skyline, L, V> { 575 typedef typename linalg_traits<V>::value_type T; 576 typedef col_matrix<slvector<T> > matrix_type; 577 }; 578 template <typename V, typename L> 579 struct temporary_col_matrix_<abstract_dense, L, V> 580 { typedef dense_matrix<typename linalg_traits<V>::value_type> matrix_type; }; 581 582 template <typename V> struct temporary_col_matrix { 583 typedef typename temporary_col_matrix_< 584 typename linalg_traits<V>::storage_type, 585 typename linalg_traits<V>::linalg_type, 586 V>::matrix_type matrix_type; 587 }; 588 589 590 591 592 template <typename S, typename L, typename V> 593 struct temporary_row_matrix_ { typedef abstract_null_type matrix_type; }; 594 template <typename V, typename L> 595 struct temporary_row_matrix_<abstract_sparse, L, V> { 596 typedef typename linalg_traits<V>::value_type T; 597 typedef row_matrix<wsvector<T> > matrix_type; 598 }; 599 template <typename V, typename L> 600 struct temporary_row_matrix_<abstract_skyline, L, V> { 601 typedef typename linalg_traits<V>::value_type T; 602 typedef row_matrix<slvector<T> > matrix_type; 603 }; 604 template <typename V, typename L> 605 struct temporary_row_matrix_<abstract_dense, L, V> 606 { typedef dense_matrix<typename linalg_traits<V>::value_type> matrix_type; }; 607 608 template <typename V> struct temporary_row_matrix { 609 typedef typename temporary_row_matrix_< 610 typename linalg_traits<V>::storage_type, 611 typename linalg_traits<V>::linalg_type, 612 V>::matrix_type matrix_type; 613 }; 614 615 616 617 /* ******************************************************************** */ 618 /* Selects a temporary dense vector type */ 619 /* V if V is a valid dense vector type, */ 620 /* std::vector if V is a reference or another type of vector */ 621 /* ******************************************************************** */ 622 623 template <typename R, typename S, typename V> 624 struct temporary_dense_vector_ { typedef abstract_null_type vector_type; }; 625 template <typename S, typename V> 626 struct temporary_dense_vector_<linalg_true, S, V> 627 { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; }; 628 template <typename V> 629 struct temporary_dense_vector_<linalg_false, abstract_sparse, V> 630 { typedef std::vector<typename linalg_traits<V>::value_type> vector_type; }; 631 template <typename V> 632 struct temporary_dense_vector_<linalg_false, abstract_skyline, 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_dense, V> 636 { typedef V vector_type; }; 637 638 template <typename V> struct temporary_dense_vector { 639 typedef typename temporary_dense_vector_<typename 640 is_a_reference<V>::reference, 641 typename linalg_traits<V>::storage_type, V>::vector_type vector_type; 642 }; 643 644 /* ******************************************************************** */ 645 /* Selects a temporary sparse vector type */ 646 /* V if V is a valid sparse vector type, */ 647 /* wsvector if V is a reference or another type of vector */ 648 /* ******************************************************************** */ 649 650 template <typename R, typename S, typename V> 651 struct temporary_sparse_vector_ { typedef abstract_null_type vector_type; }; 652 template <typename S, typename V> 653 struct temporary_sparse_vector_<linalg_true, S, V> 654 { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; }; 655 template <typename V> 656 struct temporary_sparse_vector_<linalg_false, abstract_sparse, V> 657 { typedef V vector_type; }; 658 template <typename V> 659 struct temporary_sparse_vector_<linalg_false, abstract_dense, V> 660 { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; }; 661 template <typename V> 662 struct temporary_sparse_vector_<linalg_false, abstract_skyline, V> 663 { typedef wsvector<typename linalg_traits<V>::value_type> vector_type; }; 664 665 template <typename V> struct temporary_sparse_vector { 666 typedef typename temporary_sparse_vector_<typename 667 is_a_reference<V>::reference, 668 typename linalg_traits<V>::storage_type, V>::vector_type vector_type; 669 }; 670 671 /* ******************************************************************** */ 672 /* Selects a temporary sky-line vector type */ 673 /* V if V is a valid sky-line vector type, */ 674 /* slvector if V is a reference or another type of vector */ 675 /* ******************************************************************** */ 676 677 template <typename R, typename S, typename V> 678 struct temporary_skyline_vector_ 679 { typedef abstract_null_type vector_type; }; 680 template <typename S, typename V> 681 struct temporary_skyline_vector_<linalg_true, S, V> 682 { typedef slvector<typename linalg_traits<V>::value_type> vector_type; }; 683 template <typename V> 684 struct temporary_skyline_vector_<linalg_false, abstract_skyline, V> 685 { typedef V vector_type; }; 686 template <typename V> 687 struct temporary_skyline_vector_<linalg_false, abstract_dense, V> 688 { typedef slvector<typename linalg_traits<V>::value_type> vector_type; }; 689 template <typename V> 690 struct temporary_skyline_vector_<linalg_false, abstract_sparse, V> 691 { typedef slvector<typename linalg_traits<V>::value_type> vector_type; }; 692 693 template <typename V> struct temporary_skylines_vector { 694 typedef typename temporary_skyline_vector_<typename 695 is_a_reference<V>::reference, 696 typename linalg_traits<V>::storage_type, V>::vector_type vector_type; 697 }; 698 699 /* ********************************************************************* */ 700 /* Definition & Comparison of origins. */ 701 /* ********************************************************************* */ 702 703 template <typename L> 704 typename select_return<const typename linalg_traits<L>::origin_type *, 705 typename linalg_traits<L>::origin_type *, 706 L *>::return_type 707 linalg_origin(L &l) 708 { return linalg_traits<L>::origin(linalg_cast(l)); } 709 710 template <typename L> 711 typename select_return<const typename linalg_traits<L>::origin_type *, 712 typename linalg_traits<L>::origin_type *, 713 const L *>::return_type 714 linalg_origin(const L &l) 715 { return linalg_traits<L>::origin(linalg_cast(l)); } 716 717 template <typename PT1, typename PT2> 718 bool same_porigin(PT1, PT2) { return false; } 719 720 template <typename PT> 721 bool same_porigin(PT pt1, PT pt2) { return (pt1 == pt2); } 722 723 template <typename L1, typename L2> 724 bool same_origin(const L1 &l1, const L2 &l2) 725 { return same_porigin(linalg_origin(l1), linalg_origin(l2)); } 726 727 728 /* ******************************************************************** */ 729 /* Miscellaneous */ 730 /* ******************************************************************** */ 731 732 template <typename V> inline size_type vect_size(const V &v) 733 { return linalg_traits<V>::size(v); } 734 735 template <typename MAT> inline size_type mat_nrows(const MAT &m) 736 { return linalg_traits<MAT>::nrows(m); } 737 738 template <typename MAT> inline size_type mat_ncols(const MAT &m) 739 { return linalg_traits<MAT>::ncols(m); } 740 741 742 template <typename V> inline 743 typename select_return<typename linalg_traits<V>::const_iterator, 744 typename linalg_traits<V>::iterator, V *>::return_type 745 vect_begin(V &v) 746 { return linalg_traits<V>::begin(linalg_cast(v)); } 747 748 template <typename V> inline 749 typename select_return<typename linalg_traits<V>::const_iterator, 750 typename linalg_traits<V>::iterator, const V *>::return_type 751 vect_begin(const V &v) 752 { return linalg_traits<V>::begin(linalg_cast(v)); } 753 754 template <typename V> inline 755 typename linalg_traits<V>::const_iterator 756 vect_const_begin(const V &v) 757 { return linalg_traits<V>::begin(v); } 758 759 template <typename V> inline 760 typename select_return<typename linalg_traits<V>::const_iterator, 761 typename linalg_traits<V>::iterator, V *>::return_type 762 vect_end(V &v) 763 { return linalg_traits<V>::end(linalg_cast(v)); } 764 765 template <typename V> inline 766 typename select_return<typename linalg_traits<V>::const_iterator, 767 typename linalg_traits<V>::iterator, const V *>::return_type 768 vect_end(const V &v) 769 { return linalg_traits<V>::end(linalg_cast(v)); } 770 771 template <typename V> inline 772 typename linalg_traits<V>::const_iterator 773 vect_const_end(const V &v) 774 { return linalg_traits<V>::end(v); } 775 776 template <typename M> inline 777 typename select_return<typename linalg_traits<M>::const_row_iterator, 778 typename linalg_traits<M>::row_iterator, M *>::return_type 779 mat_row_begin(M &m) { return linalg_traits<M>::row_begin(linalg_cast(m)); } 780 781 template <typename M> inline 782 typename select_return<typename linalg_traits<M>::const_row_iterator, 783 typename linalg_traits<M>::row_iterator, const M *>::return_type 784 mat_row_begin(const M &m) 785 { return linalg_traits<M>::row_begin(linalg_cast(m)); } 786 787 template <typename M> inline typename linalg_traits<M>::const_row_iterator 788 mat_row_const_begin(const M &m) 789 { return linalg_traits<M>::row_begin(m); } 790 791 template <typename M> inline 792 typename select_return<typename linalg_traits<M>::const_row_iterator, 793 typename linalg_traits<M>::row_iterator, M *>::return_type 794 mat_row_end(M &v) { 795 return linalg_traits<M>::row_end(linalg_cast(v)); 796 } 797 798 template <typename M> inline 799 typename select_return<typename linalg_traits<M>::const_row_iterator, 800 typename linalg_traits<M>::row_iterator, const M *>::return_type 801 mat_row_end(const M &v) { 802 return linalg_traits<M>::row_end(linalg_cast(v)); 803 } 804 805 template <typename M> inline 806 typename linalg_traits<M>::const_row_iterator 807 mat_row_const_end(const M &v) 808 { return linalg_traits<M>::row_end(v); } 809 810 template <typename M> inline 811 typename select_return<typename linalg_traits<M>::const_col_iterator, 812 typename linalg_traits<M>::col_iterator, M *>::return_type 813 mat_col_begin(M &v) { 814 return linalg_traits<M>::col_begin(linalg_cast(v)); 815 } 816 817 template <typename M> inline 818 typename select_return<typename linalg_traits<M>::const_col_iterator, 819 typename linalg_traits<M>::col_iterator, const M *>::return_type 820 mat_col_begin(const M &v) { 821 return linalg_traits<M>::col_begin(linalg_cast(v)); 822 } 823 824 template <typename M> inline 825 typename linalg_traits<M>::const_col_iterator 826 mat_col_const_begin(const M &v) 827 { return linalg_traits<M>::col_begin(v); } 828 829 template <typename M> inline 830 typename linalg_traits<M>::const_col_iterator 831 mat_col_const_end(const M &v) 832 { return linalg_traits<M>::col_end(v); } 833 834 template <typename M> inline 835 typename select_return<typename linalg_traits<M>::const_col_iterator, 836 typename linalg_traits<M>::col_iterator, 837 M *>::return_type 838 mat_col_end(M &m) 839 { return linalg_traits<M>::col_end(linalg_cast(m)); } 840 841 template <typename M> inline 842 typename select_return<typename linalg_traits<M>::const_col_iterator, 843 typename linalg_traits<M>::col_iterator, 844 const M *>::return_type 845 mat_col_end(const M &m) 846 { return linalg_traits<M>::col_end(linalg_cast(m)); } 847 848 template <typename MAT> inline 849 typename select_return<typename linalg_traits<MAT>::const_sub_row_type, 850 typename linalg_traits<MAT>::sub_row_type, 851 const MAT *>::return_type 852 mat_row(const MAT &m, size_type i) 853 { return linalg_traits<MAT>::row(mat_row_begin(m) + i); } 854 855 template <typename MAT> inline 856 typename select_return<typename linalg_traits<MAT>::const_sub_row_type, 857 typename linalg_traits<MAT>::sub_row_type, 858 MAT *>::return_type 859 mat_row(MAT &m, size_type i) 860 { return linalg_traits<MAT>::row(mat_row_begin(m) + i); } 861 862 template <typename MAT> inline 863 typename linalg_traits<MAT>::const_sub_row_type 864 mat_const_row(const MAT &m, size_type i) 865 { return linalg_traits<MAT>::row(mat_row_const_begin(m) + i); } 866 867 template <typename MAT> inline 868 typename select_return<typename linalg_traits<MAT>::const_sub_col_type, 869 typename linalg_traits<MAT>::sub_col_type, 870 const MAT *>::return_type 871 mat_col(const MAT &m, size_type i) 872 { return linalg_traits<MAT>::col(mat_col_begin(m) + i); } 873 874 875 template <typename MAT> inline 876 typename select_return<typename linalg_traits<MAT>::const_sub_col_type, 877 typename linalg_traits<MAT>::sub_col_type, 878 MAT *>::return_type 879 mat_col(MAT &m, size_type i) 880 { return linalg_traits<MAT>::col(mat_col_begin(m) + i); } 881 882 template <typename MAT> inline 883 typename linalg_traits<MAT>::const_sub_col_type 884 mat_const_col(const MAT &m, size_type i) 885 { return linalg_traits<MAT>::col(mat_col_const_begin(m) + i); } 886 887 /* ********************************************************************* */ 888 /* Set to begin end set to end for iterators on non-const sparse vectors.*/ 889 /* ********************************************************************* */ 890 891 template <typename IT, typename ORG, typename VECT> inline 892 void set_to_begin(IT &it, ORG o, VECT *, linalg_false) 893 { it = vect_begin(*o); } 894 895 template <typename IT, typename ORG, typename VECT> inline 896 void set_to_begin(IT &it, ORG o, const VECT *, linalg_false) 897 { it = vect_const_begin(*o); } 898 899 template <typename IT, typename ORG, typename VECT> inline 900 void set_to_end(IT &it, ORG o, VECT *, linalg_false) 901 { it = vect_end(*o); } 902 903 template <typename IT, typename ORG, typename VECT> inline 904 void set_to_end(IT &it, ORG o, const VECT *, linalg_false) 905 { it = vect_const_end(*o); } 906 907 908 template <typename IT, typename ORG, typename VECT> inline 909 void set_to_begin(IT &, ORG, VECT *, linalg_const) { } 910 911 template <typename IT, typename ORG, typename VECT> inline 912 void set_to_begin(IT &, ORG, const VECT *, linalg_const) { } 913 914 template <typename IT, typename ORG, typename VECT> inline 915 void set_to_end(IT &, ORG, VECT *, linalg_const) { } 916 917 template <typename IT, typename ORG, typename VECT> inline 918 void set_to_end(IT &, ORG, const VECT *, linalg_const) { } 919 920 921 template <typename IT, typename ORG, typename VECT> inline 922 void set_to_begin(IT &, ORG, VECT *v, linalg_modifiable) 923 { GMM_ASSERT3(!is_sparse(*v), "internal_error"); v = 0; } 924 925 template <typename IT, typename ORG, typename VECT> inline 926 void set_to_begin(IT &, ORG, const VECT *v, linalg_modifiable) 927 { GMM_ASSERT3(!is_sparse(*v), "internal_error"); v = 0; } 928 929 template <typename IT, typename ORG, typename VECT> inline 930 void set_to_end(IT &, ORG, VECT *v, linalg_modifiable) 931 { GMM_ASSERT3(!is_sparse(*v), "internal_error"); v = 0; } 932 933 template <typename IT, typename ORG, typename VECT> inline 934 void set_to_end(IT &, ORG, const VECT *v, linalg_modifiable) 935 { GMM_ASSERT3(!is_sparse(*v), "internal_error"); v = 0; } 936 937 /* ******************************************************************** */ 938 /* General index for certain algorithms. */ 939 /* ******************************************************************** */ 940 941 template<class IT> 942 size_type index_of_it(const IT &it, size_type, abstract_sparse) 943 { return it.index(); } 944 template<class IT> 945 size_type index_of_it(const IT &it, size_type, abstract_skyline) 946 { return it.index(); } 947 template<class IT> 948 size_type index_of_it(const IT &, size_type k, abstract_dense) 949 { return k; } 950 951 /* ********************************************************************* */ 952 /* Numeric limits. */ 953 /* ********************************************************************* */ 954 955 template<typename T> inline T default_tol(T) { 956 using namespace std; 957 static T tol(10); 958 if (tol == T(10)) { 959 if (numeric_limits<T>::is_specialized) 960 tol = numeric_limits<T>::epsilon(); 961 else { 962 int i=sizeof(T)/4; while(i-- > 0) tol*=T(1E-8); 963 GMM_WARNING1("The numeric type " << typeid(T).name() 964 << " has no numeric_limits defined !!\n" 965 << "Taking " << tol << " as default tolerance"); 966 } 967 } 968 return tol; 969 } 970 template<typename T> inline T default_tol(std::complex<T>) 971 { return default_tol(T()); } 972 973 template<typename T> inline T default_min(T) { 974 using namespace std; 975 static T mi(10); 976 if (mi == T(10)) { 977 if (numeric_limits<T>::is_specialized) 978 mi = std::numeric_limits<T>::min(); 979 else { 980 mi = T(0); 981 GMM_WARNING1("The numeric type " << typeid(T).name() 982 << " has no numeric_limits defined !!\n" 983 << "Taking 0 as default minimum"); 984 } 985 } 986 return mi; 987 } 988 template<typename T> inline T default_min(std::complex<T>) 989 { return default_min(T()); } 990 991 template<typename T> inline T default_max(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>::max(); 997 else { 998 mi = T(1); 999 GMM_WARNING1("The numeric type " << typeid(T).name() 1000 << " has no numeric_limits defined !!\n" 1001 << "Taking 1 as default maximum !"); 1002 } 1003 } 1004 return mi; 1005 } 1006 template<typename T> inline T default_max(std::complex<T>) 1007 { return default_max(T()); } 1008 1009 1010 /* 1011 use safe_divide to avoid NaNs when dividing very small complex 1012 numbers, for example 1013 std::complex<float>(1e-23,1e-30)/std::complex<float>(1e-23,1e-30) 1014 */ 1015 template<typename T> inline T safe_divide(T a, T b) { return a/b; } 1016 template<typename T> inline std::complex<T> 1017 safe_divide(std::complex<T> a, std::complex<T> b) { 1018 T m = std::max(gmm::abs(b.real()), gmm::abs(b.imag())); 1019 a = std::complex<T>(a.real()/m, a.imag()/m); 1020 b = std::complex<T>(b.real()/m, b.imag()/m); 1021 return a / b; 1022 } 1023 1024 1025 /* ******************************************************************** */ 1026 /* Write */ 1027 /* ******************************************************************** */ 1028 1029 template <typename T> struct cast_char_type { typedef T return_type; }; 1030 template <> struct cast_char_type<signed char> { typedef int return_type; }; 1031 template <> struct cast_char_type<unsigned char> 1032 { typedef unsigned int return_type; }; 1033 template <typename T> inline typename cast_char_type<T>::return_type 1034 cast_char(const T &c) { return typename cast_char_type<T>::return_type(c); } 1035 1036 1037 template <typename L> inline void write(std::ostream &o, const L &l) 1038 { write(o, l, typename linalg_traits<L>::linalg_type()); } 1039 1040 template <typename L> void write(std::ostream &o, const L &l, 1041 abstract_vector) { 1042 o << "vector(" << vect_size(l) << ") ["; 1043 write(o, l, typename linalg_traits<L>::storage_type()); 1044 o << " ]"; 1045 } 1046 1047 template <typename L> void write(std::ostream &o, const L &l, 1048 abstract_sparse) { 1049 typename linalg_traits<L>::const_iterator it = vect_const_begin(l), 1050 ite = vect_const_end(l); 1051 for (; it != ite; ++it) 1052 o << " (r" << it.index() << "," << cast_char(*it) << ")"; 1053 } 1054 1055 template <typename L> void write(std::ostream &o, const L &l, 1056 abstract_dense) { 1057 typename linalg_traits<L>::const_iterator it = vect_const_begin(l), 1058 ite = vect_const_end(l); 1059 if (it != ite) o << " " << cast_char(*it++); 1060 for (; it != ite; ++it) o << ", " << cast_char(*it); 1061 } 1062 1063 template <typename L> void write(std::ostream &o, const L &l, 1064 abstract_skyline) { 1065 typedef typename linalg_traits<L>::const_iterator const_iterator; 1066 const_iterator it = vect_const_begin(l), ite = vect_const_end(l); 1067 if (it != ite) { 1068 o << "<r+" << it.index() << ">"; 1069 if (it != ite) o << " " << cast_char(*it++); 1070 for (; it != ite; ++it) { o << ", " << cast_char(*it); } 1071 } 1072 } 1073 1074 template <typename L> inline void write(std::ostream &o, const L &l, 1075 abstract_matrix) { 1076 write(o, l, typename linalg_traits<L>::sub_orientation()); 1077 } 1078 1079 1080 template <typename L> void write(std::ostream &o, const L &l, 1081 row_major) { 1082 o << "matrix(" << mat_nrows(l) << ", " << mat_ncols(l) << ")" << endl; 1083 for (size_type i = 0; i < mat_nrows(l); ++i) { 1084 o << "("; 1085 write(o, mat_const_row(l, i), typename linalg_traits<L>::storage_type()); 1086 o << " )\n"; 1087 } 1088 } 1089 1090 template <typename L> inline 1091 void write(std::ostream &o, const L &l, row_and_col) 1092 { write(o, l, row_major()); } 1093 1094 template <typename L> inline 1095 void write(std::ostream &o, const L &l, col_and_row) 1096 { write(o, l, row_major()); } 1097 1098 template <typename L> void write(std::ostream &o, const L &l, col_major) { 1099 o << "matrix(" << mat_nrows(l) << ", " << mat_ncols(l) << ")" << endl; 1100 for (size_type i = 0; i < mat_nrows(l); ++i) { 1101 o << "("; 1102 if (is_sparse(l)) { // not optimized ... 1103 for (size_type j = 0; j < mat_ncols(l); ++j) 1104 if (l(i,j) != typename linalg_traits<L>::value_type(0)) 1105 o << " (r" << j << ", " << l(i,j) << ")"; 1106 } 1107 else { 1108 if (mat_ncols(l) != 0) o << ' ' << l(i, 0); 1109 for (size_type j = 1; j < mat_ncols(l); ++j) o << ", " << l(i, j); 1110 } 1111 o << " )\n"; 1112 } 1113 } 1114 1115 } 1116 1117 #endif // GMM_DEF_H__ 1118