1 // 2 // Copyright (c) 2000-2009 3 // Joerg Walter, Mathias Koch, Gunter Winkler 4 // 5 // Distributed under the Boost Software License, Version 1.0. (See 6 // accompanying file LICENSE_1_0.txt or copy at 7 // http://www.boost.org/LICENSE_1_0.txt) 8 // 9 // The authors gratefully acknowledge the support of 10 // GeNeSys mbH & Co. KG in producing this work. 11 // 12 13 #ifndef _BOOST_UBLAS_FUNCTIONAL_ 14 #define _BOOST_UBLAS_FUNCTIONAL_ 15 16 #include <functional> 17 18 #include <boost/core/ignore_unused.hpp> 19 20 #include <boost/numeric/ublas/traits.hpp> 21 #ifdef BOOST_UBLAS_USE_DUFF_DEVICE 22 #include <boost/numeric/ublas/detail/duff.hpp> 23 #endif 24 #ifdef BOOST_UBLAS_USE_SIMD 25 #include <boost/numeric/ublas/detail/raw.hpp> 26 #else 27 namespace boost { namespace numeric { namespace ublas { namespace raw { 28 }}}} 29 #endif 30 #ifdef BOOST_UBLAS_HAVE_BINDINGS 31 #include <boost/numeric/bindings/traits/std_vector.hpp> 32 #include <boost/numeric/bindings/traits/ublas_vector.hpp> 33 #include <boost/numeric/bindings/traits/ublas_matrix.hpp> 34 #include <boost/numeric/bindings/atlas/cblas.hpp> 35 #endif 36 37 #include <boost/numeric/ublas/detail/definitions.hpp> 38 39 40 41 namespace boost { namespace numeric { namespace ublas { 42 43 // Scalar functors 44 45 // Unary 46 template<class T> 47 struct scalar_unary_functor { 48 typedef T value_type; 49 typedef typename type_traits<T>::const_reference argument_type; 50 typedef typename type_traits<T>::value_type result_type; 51 }; 52 53 template<class T> 54 struct scalar_identity: 55 public scalar_unary_functor<T> { 56 typedef typename scalar_unary_functor<T>::argument_type argument_type; 57 typedef typename scalar_unary_functor<T>::result_type result_type; 58 59 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_identity60 result_type apply (argument_type t) { 61 return t; 62 } 63 }; 64 template<class T> 65 struct scalar_negate: 66 public scalar_unary_functor<T> { 67 typedef typename scalar_unary_functor<T>::argument_type argument_type; 68 typedef typename scalar_unary_functor<T>::result_type result_type; 69 70 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_negate71 result_type apply (argument_type t) { 72 return - t; 73 } 74 }; 75 template<class T> 76 struct scalar_conj: 77 public scalar_unary_functor<T> { 78 typedef typename scalar_unary_functor<T>::value_type value_type; 79 typedef typename scalar_unary_functor<T>::argument_type argument_type; 80 typedef typename scalar_unary_functor<T>::result_type result_type; 81 82 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_conj83 result_type apply (argument_type t) { 84 return type_traits<value_type>::conj (t); 85 } 86 }; 87 88 // Unary returning real 89 template<class T> 90 struct scalar_real_unary_functor { 91 typedef T value_type; 92 typedef typename type_traits<T>::const_reference argument_type; 93 typedef typename type_traits<T>::real_type result_type; 94 }; 95 96 template<class T> 97 struct scalar_real: 98 public scalar_real_unary_functor<T> { 99 typedef typename scalar_real_unary_functor<T>::value_type value_type; 100 typedef typename scalar_real_unary_functor<T>::argument_type argument_type; 101 typedef typename scalar_real_unary_functor<T>::result_type result_type; 102 103 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_real104 result_type apply (argument_type t) { 105 return type_traits<value_type>::real (t); 106 } 107 }; 108 template<class T> 109 struct scalar_imag: 110 public scalar_real_unary_functor<T> { 111 typedef typename scalar_real_unary_functor<T>::value_type value_type; 112 typedef typename scalar_real_unary_functor<T>::argument_type argument_type; 113 typedef typename scalar_real_unary_functor<T>::result_type result_type; 114 115 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_imag116 result_type apply (argument_type t) { 117 return type_traits<value_type>::imag (t); 118 } 119 }; 120 121 // Binary 122 template<class T1, class T2> 123 struct scalar_binary_functor { 124 typedef typename type_traits<T1>::const_reference argument1_type; 125 typedef typename type_traits<T2>::const_reference argument2_type; 126 typedef typename promote_traits<T1, T2>::promote_type result_type; 127 }; 128 129 template<class T1, class T2> 130 struct scalar_plus: 131 public scalar_binary_functor<T1, T2> { 132 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 133 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 134 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 135 136 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_plus137 result_type apply (argument1_type t1, argument2_type t2) { 138 return t1 + t2; 139 } 140 }; 141 template<class T1, class T2> 142 struct scalar_minus: 143 public scalar_binary_functor<T1, T2> { 144 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 145 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 146 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 147 148 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_minus149 result_type apply (argument1_type t1, argument2_type t2) { 150 return t1 - t2; 151 } 152 }; 153 template<class T1, class T2> 154 struct scalar_multiplies: 155 public scalar_binary_functor<T1, T2> { 156 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 157 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 158 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 159 160 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_multiplies161 result_type apply (argument1_type t1, argument2_type t2) { 162 return t1 * t2; 163 } 164 }; 165 template<class T1, class T2> 166 struct scalar_divides: 167 public scalar_binary_functor<T1, T2> { 168 typedef typename scalar_binary_functor<T1, T2>::argument1_type argument1_type; 169 typedef typename scalar_binary_functor<T1, T2>::argument2_type argument2_type; 170 typedef typename scalar_binary_functor<T1, T2>::result_type result_type; 171 172 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_divides173 result_type apply (argument1_type t1, argument2_type t2) { 174 return t1 / t2; 175 } 176 }; 177 178 template<class T1, class T2> 179 struct scalar_binary_assign_functor { 180 // ISSUE Remove reference to avoid reference to reference problems 181 typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type; 182 typedef typename type_traits<T2>::const_reference argument2_type; 183 }; 184 185 struct assign_tag {}; 186 struct computed_assign_tag {}; 187 188 template<class T1, class T2> 189 struct scalar_assign: 190 public scalar_binary_assign_functor<T1, T2> { 191 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 192 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 193 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 194 static const bool computed ; 195 #else 196 static const bool computed = false ; 197 #endif 198 199 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_assign200 void apply (argument1_type t1, argument2_type t2) { 201 t1 = t2; 202 } 203 204 template<class U1, class U2> 205 struct rebind { 206 typedef scalar_assign<U1, U2> other; 207 }; 208 }; 209 210 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 211 template<class T1, class T2> 212 const bool scalar_assign<T1,T2>::computed = false; 213 #endif 214 215 template<class T1, class T2> 216 struct scalar_plus_assign: 217 public scalar_binary_assign_functor<T1, T2> { 218 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 219 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 220 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 221 static const bool computed ; 222 #else 223 static const bool computed = true ; 224 #endif 225 226 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_plus_assign227 void apply (argument1_type t1, argument2_type t2) { 228 t1 += t2; 229 } 230 231 template<class U1, class U2> 232 struct rebind { 233 typedef scalar_plus_assign<U1, U2> other; 234 }; 235 }; 236 237 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 238 template<class T1, class T2> 239 const bool scalar_plus_assign<T1,T2>::computed = true; 240 #endif 241 242 template<class T1, class T2> 243 struct scalar_minus_assign: 244 public scalar_binary_assign_functor<T1, T2> { 245 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 246 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 247 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 248 static const bool computed ; 249 #else 250 static const bool computed = true ; 251 #endif 252 253 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_minus_assign254 void apply (argument1_type t1, argument2_type t2) { 255 t1 -= t2; 256 } 257 258 template<class U1, class U2> 259 struct rebind { 260 typedef scalar_minus_assign<U1, U2> other; 261 }; 262 }; 263 264 #if BOOST_WORKAROUND( __IBMCPP__, <=600 ) 265 template<class T1, class T2> 266 const bool scalar_minus_assign<T1,T2>::computed = true; 267 #endif 268 269 template<class T1, class T2> 270 struct scalar_multiplies_assign: 271 public scalar_binary_assign_functor<T1, T2> { 272 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 273 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 274 static const bool computed = true; 275 276 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_multiplies_assign277 void apply (argument1_type t1, argument2_type t2) { 278 t1 *= t2; 279 } 280 281 template<class U1, class U2> 282 struct rebind { 283 typedef scalar_multiplies_assign<U1, U2> other; 284 }; 285 }; 286 template<class T1, class T2> 287 struct scalar_divides_assign: 288 public scalar_binary_assign_functor<T1, T2> { 289 typedef typename scalar_binary_assign_functor<T1, T2>::argument1_type argument1_type; 290 typedef typename scalar_binary_assign_functor<T1, T2>::argument2_type argument2_type; 291 static const bool computed ; 292 293 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_divides_assign294 void apply (argument1_type t1, argument2_type t2) { 295 t1 /= t2; 296 } 297 298 template<class U1, class U2> 299 struct rebind { 300 typedef scalar_divides_assign<U1, U2> other; 301 }; 302 }; 303 template<class T1, class T2> 304 const bool scalar_divides_assign<T1,T2>::computed = true; 305 306 template<class T1, class T2> 307 struct scalar_binary_swap_functor { 308 typedef typename type_traits<typename boost::remove_reference<T1>::type>::reference argument1_type; 309 typedef typename type_traits<typename boost::remove_reference<T2>::type>::reference argument2_type; 310 }; 311 312 template<class T1, class T2> 313 struct scalar_swap: 314 public scalar_binary_swap_functor<T1, T2> { 315 typedef typename scalar_binary_swap_functor<T1, T2>::argument1_type argument1_type; 316 typedef typename scalar_binary_swap_functor<T1, T2>::argument2_type argument2_type; 317 318 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::scalar_swap319 void apply (argument1_type t1, argument2_type t2) { 320 std::swap (t1, t2); 321 } 322 323 template<class U1, class U2> 324 struct rebind { 325 typedef scalar_swap<U1, U2> other; 326 }; 327 }; 328 329 // Vector functors 330 331 // Unary returning scalar 332 template<class V> 333 struct vector_scalar_unary_functor { 334 typedef typename V::value_type value_type; 335 typedef typename V::value_type result_type; 336 }; 337 338 template<class V> 339 struct vector_sum: 340 public vector_scalar_unary_functor<V> { 341 typedef typename vector_scalar_unary_functor<V>::value_type value_type; 342 typedef typename vector_scalar_unary_functor<V>::result_type result_type; 343 344 template<class E> 345 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_sum346 result_type apply (const vector_expression<E> &e) { 347 result_type t = result_type (0); 348 typedef typename E::size_type vector_size_type; 349 vector_size_type size (e ().size ()); 350 for (vector_size_type i = 0; i < size; ++ i) 351 t += e () (i); 352 return t; 353 } 354 // Dense case 355 template<class D, class I> 356 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_sum357 result_type apply (D size, I it) { 358 result_type t = result_type (0); 359 while (-- size >= 0) 360 t += *it, ++ it; 361 return t; 362 } 363 // Sparse case 364 template<class I> 365 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_sum366 result_type apply (I it, const I &it_end) { 367 result_type t = result_type (0); 368 while (it != it_end) 369 t += *it, ++ it; 370 return t; 371 } 372 }; 373 374 // Unary returning real scalar 375 template<class V> 376 struct vector_scalar_real_unary_functor { 377 typedef typename V::value_type value_type; 378 typedef typename type_traits<value_type>::real_type real_type; 379 typedef real_type result_type; 380 }; 381 382 template<class V> 383 struct vector_norm_1: 384 public vector_scalar_real_unary_functor<V> { 385 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 386 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 387 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 388 389 template<class E> 390 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_1391 result_type apply (const vector_expression<E> &e) { 392 real_type t = real_type (); 393 typedef typename E::size_type vector_size_type; 394 vector_size_type size (e ().size ()); 395 for (vector_size_type i = 0; i < size; ++ i) { 396 real_type u (type_traits<value_type>::type_abs (e () (i))); 397 t += u; 398 } 399 return t; 400 } 401 // Dense case 402 template<class D, class I> 403 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_1404 result_type apply (D size, I it) { 405 real_type t = real_type (); 406 while (-- size >= 0) { 407 real_type u (type_traits<value_type>::norm_1 (*it)); 408 t += u; 409 ++ it; 410 } 411 return t; 412 } 413 // Sparse case 414 template<class I> 415 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_1416 result_type apply (I it, const I &it_end) { 417 real_type t = real_type (); 418 while (it != it_end) { 419 real_type u (type_traits<value_type>::norm_1 (*it)); 420 t += u; 421 ++ it; 422 } 423 return t; 424 } 425 }; 426 template<class V> 427 struct vector_norm_2: 428 public vector_scalar_real_unary_functor<V> { 429 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 430 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 431 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 432 433 template<class E> 434 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_2435 result_type apply (const vector_expression<E> &e) { 436 typedef typename E::size_type vector_size_type; 437 vector_size_type size (e ().size ()); 438 #ifndef BOOST_UBLAS_SCALED_NORM 439 real_type t = real_type (); 440 for (vector_size_type i = 0; i < size; ++ i) { 441 real_type u (type_traits<value_type>::norm_2 (e () (i))); 442 t += u * u; 443 } 444 return static_cast<result_type>(type_traits<real_type>::type_sqrt (t)); 445 #else 446 real_type scale = real_type (); 447 real_type sum_squares (1); 448 for (vector_size_type i = 0; i < size; ++ i) { 449 real_type u (type_traits<value_type>::norm_2 (e () (i))); 450 if ( real_type () /* zero */ == u ) continue; 451 if (scale < u) { 452 real_type v (scale / u); 453 sum_squares = sum_squares * v * v + real_type (1); 454 scale = u; 455 } else { 456 real_type v (u / scale); 457 sum_squares += v * v; 458 } 459 } 460 return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares)); 461 #endif 462 } 463 // Dense case 464 template<class D, class I> 465 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_2466 result_type apply (D size, I it) { 467 #ifndef BOOST_UBLAS_SCALED_NORM 468 real_type t = real_type (); 469 while (-- size >= 0) { 470 real_type u (type_traits<value_type>::norm_2 (*it)); 471 t += u * u; 472 ++ it; 473 } 474 return static_cast<result_type>(type_traits<real_type>::type_sqrt (t)); 475 #else 476 real_type scale = real_type (); 477 real_type sum_squares (1); 478 while (-- size >= 0) { 479 real_type u (type_traits<value_type>::norm_2 (*it)); 480 if (scale < u) { 481 real_type v (scale / u); 482 sum_squares = sum_squares * v * v + real_type (1); 483 scale = u; 484 } else { 485 real_type v (u / scale); 486 sum_squares += v * v; 487 } 488 ++ it; 489 } 490 return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares)); 491 #endif 492 } 493 // Sparse case 494 template<class I> 495 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_2496 result_type apply (I it, const I &it_end) { 497 #ifndef BOOST_UBLAS_SCALED_NORM 498 real_type t = real_type (); 499 while (it != it_end) { 500 real_type u (type_traits<value_type>::norm_2 (*it)); 501 t += u * u; 502 ++ it; 503 } 504 return static_cast<result_type>(type_traits<real_type>::type_sqrt (t)); 505 #else 506 real_type scale = real_type (); 507 real_type sum_squares (1); 508 while (it != it_end) { 509 real_type u (type_traits<value_type>::norm_2 (*it)); 510 if (scale < u) { 511 real_type v (scale / u); 512 sum_squares = sum_squares * v * v + real_type (1); 513 scale = u; 514 } else { 515 real_type v (u / scale); 516 sum_squares += v * v; 517 } 518 ++ it; 519 } 520 return static_cast<result_type>(scale * type_traits<real_type>::type_sqrt (sum_squares)); 521 #endif 522 } 523 }; 524 525 template<class V> 526 struct vector_norm_2_square : 527 public vector_scalar_real_unary_functor<V> { 528 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 529 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 530 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 531 532 template<class E> 533 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_2_square534 result_type apply (const vector_expression<E> &e) { 535 real_type t = real_type (); 536 typedef typename E::size_type vector_size_type; 537 vector_size_type size (e ().size ()); 538 for (vector_size_type i = 0; i < size; ++ i) { 539 real_type u (type_traits<value_type>::norm_2 (e () (i))); 540 t += u * u; 541 } 542 return t; 543 } 544 // Dense case 545 template<class D, class I> 546 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_2_square547 result_type apply (D size, I it) { 548 real_type t = real_type (); 549 while (-- size >= 0) { 550 real_type u (type_traits<value_type>::norm_2 (*it)); 551 t += u * u; 552 ++ it; 553 } 554 return t; 555 } 556 // Sparse case 557 template<class I> 558 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_2_square559 result_type apply (I it, const I &it_end) { 560 real_type t = real_type (); 561 while (it != it_end) { 562 real_type u (type_traits<value_type>::norm_2 (*it)); 563 t += u * u; 564 ++ it; 565 } 566 return t; 567 } 568 }; 569 570 template<class V> 571 struct vector_norm_inf: 572 public vector_scalar_real_unary_functor<V> { 573 typedef typename vector_scalar_real_unary_functor<V>::value_type value_type; 574 typedef typename vector_scalar_real_unary_functor<V>::real_type real_type; 575 typedef typename vector_scalar_real_unary_functor<V>::result_type result_type; 576 577 template<class E> 578 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_inf579 result_type apply (const vector_expression<E> &e) { 580 real_type t = real_type (); 581 typedef typename E::size_type vector_size_type; 582 vector_size_type size (e ().size ()); 583 for (vector_size_type i = 0; i < size; ++ i) { 584 real_type u (type_traits<value_type>::norm_inf (e () (i))); 585 if (u > t) 586 t = u; 587 } 588 return t; 589 } 590 // Dense case 591 template<class D, class I> 592 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_inf593 result_type apply (D size, I it) { 594 real_type t = real_type (); 595 while (-- size >= 0) { 596 real_type u (type_traits<value_type>::norm_inf (*it)); 597 if (u > t) 598 t = u; 599 ++ it; 600 } 601 return t; 602 } 603 // Sparse case 604 template<class I> 605 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_norm_inf606 result_type apply (I it, const I &it_end) { 607 real_type t = real_type (); 608 while (it != it_end) { 609 real_type u (type_traits<value_type>::norm_inf (*it)); 610 if (u > t) 611 t = u; 612 ++ it; 613 } 614 return t; 615 } 616 }; 617 618 // Unary returning index 619 template<class V> 620 struct vector_scalar_index_unary_functor { 621 typedef typename V::value_type value_type; 622 typedef typename type_traits<value_type>::real_type real_type; 623 typedef typename V::size_type result_type; 624 }; 625 626 template<class V> 627 struct vector_index_norm_inf: 628 public vector_scalar_index_unary_functor<V> { 629 typedef typename vector_scalar_index_unary_functor<V>::value_type value_type; 630 typedef typename vector_scalar_index_unary_functor<V>::real_type real_type; 631 typedef typename vector_scalar_index_unary_functor<V>::result_type result_type; 632 633 template<class E> 634 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_index_norm_inf635 result_type apply (const vector_expression<E> &e) { 636 // ISSUE For CBLAS compatibility return 0 index in empty case 637 result_type i_norm_inf (0); 638 real_type t = real_type (); 639 typedef typename E::size_type vector_size_type; 640 vector_size_type size (e ().size ()); 641 for (vector_size_type i = 0; i < size; ++ i) { 642 real_type u (type_traits<value_type>::norm_inf (e () (i))); 643 if (u > t) { 644 i_norm_inf = i; 645 t = u; 646 } 647 } 648 return i_norm_inf; 649 } 650 // Dense case 651 template<class D, class I> 652 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_index_norm_inf653 result_type apply (D size, I it) { 654 // ISSUE For CBLAS compatibility return 0 index in empty case 655 result_type i_norm_inf (0); 656 real_type t = real_type (); 657 while (-- size >= 0) { 658 real_type u (type_traits<value_type>::norm_inf (*it)); 659 if (u > t) { 660 i_norm_inf = it.index (); 661 t = u; 662 } 663 ++ it; 664 } 665 return i_norm_inf; 666 } 667 // Sparse case 668 template<class I> 669 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_index_norm_inf670 result_type apply (I it, const I &it_end) { 671 // ISSUE For CBLAS compatibility return 0 index in empty case 672 result_type i_norm_inf (0); 673 real_type t = real_type (); 674 while (it != it_end) { 675 real_type u (type_traits<value_type>::norm_inf (*it)); 676 if (u > t) { 677 i_norm_inf = it.index (); 678 t = u; 679 } 680 ++ it; 681 } 682 return i_norm_inf; 683 } 684 }; 685 686 // Binary returning scalar 687 template<class V1, class V2, class TV> 688 struct vector_scalar_binary_functor { 689 typedef TV value_type; 690 typedef TV result_type; 691 }; 692 693 template<class V1, class V2, class TV> 694 struct vector_inner_prod: 695 public vector_scalar_binary_functor<V1, V2, TV> { 696 typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type; 697 typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type; 698 699 template<class C1, class C2> 700 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_inner_prod701 result_type apply (const vector_container<C1> &c1, 702 const vector_container<C2> &c2) { 703 #ifdef BOOST_UBLAS_USE_SIMD 704 using namespace raw; 705 typedef typename C1::size_type vector_size_type; 706 vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ())); 707 const typename V1::value_type *data1 = data_const (c1 ()); 708 const typename V1::value_type *data2 = data_const (c2 ()); 709 vector_size_type s1 = stride (c1 ()); 710 vector_size_type s2 = stride (c2 ()); 711 result_type t = result_type (0); 712 if (s1 == 1 && s2 == 1) { 713 for (vector_size_type i = 0; i < size; ++ i) 714 t += data1 [i] * data2 [i]; 715 } else if (s2 == 1) { 716 for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1) 717 t += data1 [i1] * data2 [i]; 718 } else if (s1 == 1) { 719 for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2) 720 t += data1 [i] * data2 [i2]; 721 } else { 722 for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2) 723 t += data1 [i1] * data2 [i2]; 724 } 725 return t; 726 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 727 return boost::numeric::bindings::atlas::dot (c1 (), c2 ()); 728 #else 729 return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2)); 730 #endif 731 } 732 template<class E1, class E2> 733 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_inner_prod734 result_type apply (const vector_expression<E1> &e1, 735 const vector_expression<E2> &e2) { 736 typedef typename E1::size_type vector_size_type; 737 vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ())); 738 result_type t = result_type (0); 739 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 740 for (vector_size_type i = 0; i < size; ++ i) 741 t += e1 () (i) * e2 () (i); 742 #else 743 vector_size_type i (0); 744 DD (size, 4, r, (t += e1 () (i) * e2 () (i), ++ i)); 745 #endif 746 return t; 747 } 748 // Dense case 749 template<class D, class I1, class I2> 750 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_inner_prod751 result_type apply (D size, I1 it1, I2 it2) { 752 result_type t = result_type (0); 753 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 754 while (-- size >= 0) 755 t += *it1 * *it2, ++ it1, ++ it2; 756 #else 757 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 758 #endif 759 return t; 760 } 761 // Packed case 762 template<class I1, class I2> 763 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_inner_prod764 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { 765 result_type t = result_type (0); 766 typedef typename I1::difference_type vector_difference_type; 767 vector_difference_type it1_size (it1_end - it1); 768 vector_difference_type it2_size (it2_end - it2); 769 vector_difference_type diff (0); 770 if (it1_size > 0 && it2_size > 0) 771 diff = it2.index () - it1.index (); 772 if (diff != 0) { 773 vector_difference_type size = (std::min) (diff, it1_size); 774 if (size > 0) { 775 it1 += size; 776 it1_size -= size; 777 diff -= size; 778 } 779 size = (std::min) (- diff, it2_size); 780 if (size > 0) { 781 it2 += size; 782 it2_size -= size; 783 diff += size; 784 } 785 } 786 vector_difference_type size ((std::min) (it1_size, it2_size)); 787 while (-- size >= 0) 788 t += *it1 * *it2, ++ it1, ++ it2; 789 return t; 790 } 791 // Sparse case 792 template<class I1, class I2> 793 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::vector_inner_prod794 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { 795 result_type t = result_type (0); 796 if (it1 != it1_end && it2 != it2_end) { 797 for (;;) { 798 if (it1.index () == it2.index ()) { 799 t += *it1 * *it2, ++ it1, ++ it2; 800 if (it1 == it1_end || it2 == it2_end) 801 break; 802 } else if (it1.index () < it2.index ()) { 803 increment (it1, it1_end, it2.index () - it1.index ()); 804 if (it1 == it1_end) 805 break; 806 } else if (it1.index () > it2.index ()) { 807 increment (it2, it2_end, it1.index () - it2.index ()); 808 if (it2 == it2_end) 809 break; 810 } 811 } 812 } 813 return t; 814 } 815 }; 816 817 // Matrix functors 818 819 // Binary returning vector 820 template<class M1, class M2, class TV> 821 struct matrix_vector_binary_functor { 822 typedef typename M1::size_type size_type; 823 typedef typename M1::difference_type difference_type; 824 typedef TV value_type; 825 typedef TV result_type; 826 }; 827 828 template<class M1, class M2, class TV> 829 struct matrix_vector_prod1: 830 public matrix_vector_binary_functor<M1, M2, TV> { 831 typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type; 832 typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type; 833 typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type; 834 typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type; 835 836 template<class C1, class C2> 837 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1838 result_type apply (const matrix_container<C1> &c1, 839 const vector_container<C2> &c2, 840 size_type i) { 841 #ifdef BOOST_UBLAS_USE_SIMD 842 using namespace raw; 843 size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().size ()); 844 const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); 845 const typename M2::value_type *data2 = data_const (c2 ()); 846 size_type s1 = stride2 (c1 ()); 847 size_type s2 = stride (c2 ()); 848 result_type t = result_type (0); 849 if (s1 == 1 && s2 == 1) { 850 for (size_type j = 0; j < size; ++ j) 851 t += data1 [j] * data2 [j]; 852 } else if (s2 == 1) { 853 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) 854 t += data1 [j1] * data2 [j]; 855 } else if (s1 == 1) { 856 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) 857 t += data1 [j] * data2 [j2]; 858 } else { 859 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) 860 t += data1 [j1] * data2 [j2]; 861 } 862 return t; 863 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 864 return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ()); 865 #else 866 return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2, i)); 867 #endif 868 } 869 template<class E1, class E2> 870 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1871 result_type apply (const matrix_expression<E1> &e1, 872 const vector_expression<E2> &e2, 873 size_type i) { 874 size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size ()); 875 result_type t = result_type (0); 876 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 877 for (size_type j = 0; j < size; ++ j) 878 t += e1 () (i, j) * e2 () (j); 879 #else 880 size_type j (0); 881 DD (size, 4, r, (t += e1 () (i, j) * e2 () (j), ++ j)); 882 #endif 883 return t; 884 } 885 // Dense case 886 template<class I1, class I2> 887 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1888 result_type apply (difference_type size, I1 it1, I2 it2) { 889 result_type t = result_type (0); 890 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 891 while (-- size >= 0) 892 t += *it1 * *it2, ++ it1, ++ it2; 893 #else 894 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 895 #endif 896 return t; 897 } 898 // Packed case 899 template<class I1, class I2> 900 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1901 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { 902 result_type t = result_type (0); 903 difference_type it1_size (it1_end - it1); 904 difference_type it2_size (it2_end - it2); 905 difference_type diff (0); 906 if (it1_size > 0 && it2_size > 0) 907 diff = it2.index () - it1.index2 (); 908 if (diff != 0) { 909 difference_type size = (std::min) (diff, it1_size); 910 if (size > 0) { 911 it1 += size; 912 it1_size -= size; 913 diff -= size; 914 } 915 size = (std::min) (- diff, it2_size); 916 if (size > 0) { 917 it2 += size; 918 it2_size -= size; 919 diff += size; 920 } 921 } 922 difference_type size ((std::min) (it1_size, it2_size)); 923 while (-- size >= 0) 924 t += *it1 * *it2, ++ it1, ++ it2; 925 return t; 926 } 927 // Sparse case 928 template<class I1, class I2> 929 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1930 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 931 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { 932 result_type t = result_type (0); 933 if (it1 != it1_end && it2 != it2_end) { 934 size_type it1_index = it1.index2 (), it2_index = it2.index (); 935 for (;;) { 936 difference_type compare = it1_index - it2_index; 937 if (compare == 0) { 938 t += *it1 * *it2, ++ it1, ++ it2; 939 if (it1 != it1_end && it2 != it2_end) { 940 it1_index = it1.index2 (); 941 it2_index = it2.index (); 942 } else 943 break; 944 } else if (compare < 0) { 945 increment (it1, it1_end, - compare); 946 if (it1 != it1_end) 947 it1_index = it1.index2 (); 948 else 949 break; 950 } else if (compare > 0) { 951 increment (it2, it2_end, compare); 952 if (it2 != it2_end) 953 it2_index = it2.index (); 954 else 955 break; 956 } 957 } 958 } 959 return t; 960 } 961 // Sparse packed case 962 template<class I1, class I2> 963 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1964 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, 965 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { 966 result_type t = result_type (0); 967 while (it1 != it1_end) { 968 t += *it1 * it2 () (it1.index2 ()); 969 ++ it1; 970 } 971 return t; 972 } 973 // Packed sparse case 974 template<class I1, class I2> 975 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1976 result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, 977 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { 978 result_type t = result_type (0); 979 while (it2 != it2_end) { 980 t += it1 () (it1.index1 (), it2.index ()) * *it2; 981 ++ it2; 982 } 983 return t; 984 } 985 // Another dispatcher 986 template<class I1, class I2> 987 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod1988 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 989 sparse_bidirectional_iterator_tag) { 990 typedef typename I1::iterator_category iterator1_category; 991 typedef typename I2::iterator_category iterator2_category; 992 return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); 993 } 994 }; 995 996 template<class M1, class M2, class TV> 997 struct matrix_vector_prod2: 998 public matrix_vector_binary_functor<M1, M2, TV> { 999 typedef typename matrix_vector_binary_functor<M1, M2, TV>::size_type size_type; 1000 typedef typename matrix_vector_binary_functor<M1, M2, TV>::difference_type difference_type; 1001 typedef typename matrix_vector_binary_functor<M1, M2, TV>::value_type value_type; 1002 typedef typename matrix_vector_binary_functor<M1, M2, TV>::result_type result_type; 1003 1004 template<class C1, class C2> 1005 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21006 result_type apply (const vector_container<C1> &c1, 1007 const matrix_container<C2> &c2, 1008 size_type i) { 1009 #ifdef BOOST_UBLAS_USE_SIMD 1010 using namespace raw; 1011 size_type size = BOOST_UBLAS_SAME (c1 ().size (), c2 ().size1 ()); 1012 const typename M1::value_type *data1 = data_const (c1 ()); 1013 const typename M2::value_type *data2 = data_const (c2 ()) + i * stride2 (c2 ()); 1014 size_type s1 = stride (c1 ()); 1015 size_type s2 = stride1 (c2 ()); 1016 result_type t = result_type (0); 1017 if (s1 == 1 && s2 == 1) { 1018 for (size_type j = 0; j < size; ++ j) 1019 t += data1 [j] * data2 [j]; 1020 } else if (s2 == 1) { 1021 for (size_type j = 0, j1 = 0; j < size; ++ j, j1 += s1) 1022 t += data1 [j1] * data2 [j]; 1023 } else if (s1 == 1) { 1024 for (size_type j = 0, j2 = 0; j < size; ++ j, j2 += s2) 1025 t += data1 [j] * data2 [j2]; 1026 } else { 1027 for (size_type j = 0, j1 = 0, j2 = 0; j < size; ++ j, j1 += s1, j2 += s2) 1028 t += data1 [j1] * data2 [j2]; 1029 } 1030 return t; 1031 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 1032 return boost::numeric::bindings::atlas::dot (c1 (), c2 ().column (i)); 1033 #else 1034 return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); 1035 #endif 1036 } 1037 template<class E1, class E2> 1038 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21039 result_type apply (const vector_expression<E1> &e1, 1040 const matrix_expression<E2> &e2, 1041 size_type i) { 1042 size_type size = BOOST_UBLAS_SAME (e1 ().size (), e2 ().size1 ()); 1043 result_type t = result_type (0); 1044 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 1045 for (size_type j = 0; j < size; ++ j) 1046 t += e1 () (j) * e2 () (j, i); 1047 #else 1048 size_type j (0); 1049 DD (size, 4, r, (t += e1 () (j) * e2 () (j, i), ++ j)); 1050 #endif 1051 return t; 1052 } 1053 // Dense case 1054 template<class I1, class I2> 1055 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21056 result_type apply (difference_type size, I1 it1, I2 it2) { 1057 result_type t = result_type (0); 1058 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 1059 while (-- size >= 0) 1060 t += *it1 * *it2, ++ it1, ++ it2; 1061 #else 1062 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 1063 #endif 1064 return t; 1065 } 1066 // Packed case 1067 template<class I1, class I2> 1068 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21069 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) { 1070 result_type t = result_type (0); 1071 difference_type it1_size (it1_end - it1); 1072 difference_type it2_size (it2_end - it2); 1073 difference_type diff (0); 1074 if (it1_size > 0 && it2_size > 0) 1075 diff = it2.index1 () - it1.index (); 1076 if (diff != 0) { 1077 difference_type size = (std::min) (diff, it1_size); 1078 if (size > 0) { 1079 it1 += size; 1080 it1_size -= size; 1081 diff -= size; 1082 } 1083 size = (std::min) (- diff, it2_size); 1084 if (size > 0) { 1085 it2 += size; 1086 it2_size -= size; 1087 diff += size; 1088 } 1089 } 1090 difference_type size ((std::min) (it1_size, it2_size)); 1091 while (-- size >= 0) 1092 t += *it1 * *it2, ++ it1, ++ it2; 1093 return t; 1094 } 1095 // Sparse case 1096 template<class I1, class I2> 1097 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21098 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 1099 sparse_bidirectional_iterator_tag, sparse_bidirectional_iterator_tag) { 1100 result_type t = result_type (0); 1101 if (it1 != it1_end && it2 != it2_end) { 1102 size_type it1_index = it1.index (), it2_index = it2.index1 (); 1103 for (;;) { 1104 difference_type compare = it1_index - it2_index; 1105 if (compare == 0) { 1106 t += *it1 * *it2, ++ it1, ++ it2; 1107 if (it1 != it1_end && it2 != it2_end) { 1108 it1_index = it1.index (); 1109 it2_index = it2.index1 (); 1110 } else 1111 break; 1112 } else if (compare < 0) { 1113 increment (it1, it1_end, - compare); 1114 if (it1 != it1_end) 1115 it1_index = it1.index (); 1116 else 1117 break; 1118 } else if (compare > 0) { 1119 increment (it2, it2_end, compare); 1120 if (it2 != it2_end) 1121 it2_index = it2.index1 (); 1122 else 1123 break; 1124 } 1125 } 1126 } 1127 return t; 1128 } 1129 // Packed sparse case 1130 template<class I1, class I2> 1131 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21132 result_type apply (I1 it1, const I1 &/* it1_end */, I2 it2, const I2 &it2_end, 1133 packed_random_access_iterator_tag, sparse_bidirectional_iterator_tag) { 1134 result_type t = result_type (0); 1135 while (it2 != it2_end) { 1136 t += it1 () (it2.index1 ()) * *it2; 1137 ++ it2; 1138 } 1139 return t; 1140 } 1141 // Sparse packed case 1142 template<class I1, class I2> 1143 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21144 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &/* it2_end */, 1145 sparse_bidirectional_iterator_tag, packed_random_access_iterator_tag) { 1146 result_type t = result_type (0); 1147 while (it1 != it1_end) { 1148 t += *it1 * it2 () (it1.index (), it2.index2 ()); 1149 ++ it1; 1150 } 1151 return t; 1152 } 1153 // Another dispatcher 1154 template<class I1, class I2> 1155 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_vector_prod21156 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, 1157 sparse_bidirectional_iterator_tag) { 1158 typedef typename I1::iterator_category iterator1_category; 1159 typedef typename I2::iterator_category iterator2_category; 1160 return apply (it1, it1_end, it2, it2_end, iterator1_category (), iterator2_category ()); 1161 } 1162 }; 1163 1164 // Binary returning matrix 1165 template<class M1, class M2, class TV> 1166 struct matrix_matrix_binary_functor { 1167 typedef typename M1::size_type size_type; 1168 typedef typename M1::difference_type difference_type; 1169 typedef TV value_type; 1170 typedef TV result_type; 1171 }; 1172 1173 template<class M1, class M2, class TV> 1174 struct matrix_matrix_prod: 1175 public matrix_matrix_binary_functor<M1, M2, TV> { 1176 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::size_type size_type; 1177 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::difference_type difference_type; 1178 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::value_type value_type; 1179 typedef typename matrix_matrix_binary_functor<M1, M2, TV>::result_type result_type; 1180 1181 template<class C1, class C2> 1182 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_matrix_prod1183 result_type apply (const matrix_container<C1> &c1, 1184 const matrix_container<C2> &c2, 1185 size_type i, size_type j) { 1186 #ifdef BOOST_UBLAS_USE_SIMD 1187 using namespace raw; 1188 size_type size = BOOST_UBLAS_SAME (c1 ().size2 (), c2 ().sizc1 ()); 1189 const typename M1::value_type *data1 = data_const (c1 ()) + i * stride1 (c1 ()); 1190 const typename M2::value_type *data2 = data_const (c2 ()) + j * stride2 (c2 ()); 1191 size_type s1 = stride2 (c1 ()); 1192 size_type s2 = stride1 (c2 ()); 1193 result_type t = result_type (0); 1194 if (s1 == 1 && s2 == 1) { 1195 for (size_type k = 0; k < size; ++ k) 1196 t += data1 [k] * data2 [k]; 1197 } else if (s2 == 1) { 1198 for (size_type k = 0, k1 = 0; k < size; ++ k, k1 += s1) 1199 t += data1 [k1] * data2 [k]; 1200 } else if (s1 == 1) { 1201 for (size_type k = 0, k2 = 0; k < size; ++ k, k2 += s2) 1202 t += data1 [k] * data2 [k2]; 1203 } else { 1204 for (size_type k = 0, k1 = 0, k2 = 0; k < size; ++ k, k1 += s1, k2 += s2) 1205 t += data1 [k1] * data2 [k2]; 1206 } 1207 return t; 1208 #elif defined(BOOST_UBLAS_HAVE_BINDINGS) 1209 return boost::numeric::bindings::atlas::dot (c1 ().row (i), c2 ().column (j)); 1210 #else 1211 boost::ignore_unused(j); 1212 return apply (static_cast<const matrix_expression<C1> > (c1), static_cast<const matrix_expression<C2> > (c2, i)); 1213 #endif 1214 } 1215 template<class E1, class E2> 1216 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_matrix_prod1217 result_type apply (const matrix_expression<E1> &e1, 1218 const matrix_expression<E2> &e2, 1219 size_type i, size_type j) { 1220 size_type size = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); 1221 result_type t = result_type (0); 1222 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 1223 for (size_type k = 0; k < size; ++ k) 1224 t += e1 () (i, k) * e2 () (k, j); 1225 #else 1226 size_type k (0); 1227 DD (size, 4, r, (t += e1 () (i, k) * e2 () (k, j), ++ k)); 1228 #endif 1229 return t; 1230 } 1231 // Dense case 1232 template<class I1, class I2> 1233 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_matrix_prod1234 result_type apply (difference_type size, I1 it1, I2 it2) { 1235 result_type t = result_type (0); 1236 #ifndef BOOST_UBLAS_USE_DUFF_DEVICE 1237 while (-- size >= 0) 1238 t += *it1 * *it2, ++ it1, ++ it2; 1239 #else 1240 DD (size, 4, r, (t += *it1 * *it2, ++ it1, ++ it2)); 1241 #endif 1242 return t; 1243 } 1244 // Packed case 1245 template<class I1, class I2> 1246 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_matrix_prod1247 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, packed_random_access_iterator_tag) { 1248 result_type t = result_type (0); 1249 difference_type it1_size (it1_end - it1); 1250 difference_type it2_size (it2_end - it2); 1251 difference_type diff (0); 1252 if (it1_size > 0 && it2_size > 0) 1253 diff = it2.index1 () - it1.index2 (); 1254 if (diff != 0) { 1255 difference_type size = (std::min) (diff, it1_size); 1256 if (size > 0) { 1257 it1 += size; 1258 it1_size -= size; 1259 diff -= size; 1260 } 1261 size = (std::min) (- diff, it2_size); 1262 if (size > 0) { 1263 it2 += size; 1264 it2_size -= size; 1265 diff += size; 1266 } 1267 } 1268 difference_type size ((std::min) (it1_size, it2_size)); 1269 while (-- size >= 0) 1270 t += *it1 * *it2, ++ it1, ++ it2; 1271 return t; 1272 } 1273 // Sparse case 1274 template<class I1, class I2> 1275 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_matrix_prod1276 result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) { 1277 result_type t = result_type (0); 1278 if (it1 != it1_end && it2 != it2_end) { 1279 size_type it1_index = it1.index2 (), it2_index = it2.index1 (); 1280 for (;;) { 1281 difference_type compare = difference_type (it1_index - it2_index); 1282 if (compare == 0) { 1283 t += *it1 * *it2, ++ it1, ++ it2; 1284 if (it1 != it1_end && it2 != it2_end) { 1285 it1_index = it1.index2 (); 1286 it2_index = it2.index1 (); 1287 } else 1288 break; 1289 } else if (compare < 0) { 1290 increment (it1, it1_end, - compare); 1291 if (it1 != it1_end) 1292 it1_index = it1.index2 (); 1293 else 1294 break; 1295 } else if (compare > 0) { 1296 increment (it2, it2_end, compare); 1297 if (it2 != it2_end) 1298 it2_index = it2.index1 (); 1299 else 1300 break; 1301 } 1302 } 1303 } 1304 return t; 1305 } 1306 }; 1307 1308 // Unary returning scalar norm 1309 template<class M> 1310 struct matrix_scalar_real_unary_functor { 1311 typedef typename M::value_type value_type; 1312 typedef typename type_traits<value_type>::real_type real_type; 1313 typedef real_type result_type; 1314 }; 1315 1316 template<class M> 1317 struct matrix_norm_1: 1318 public matrix_scalar_real_unary_functor<M> { 1319 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; 1320 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; 1321 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; 1322 1323 template<class E> 1324 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_norm_11325 result_type apply (const matrix_expression<E> &e) { 1326 real_type t = real_type (); 1327 typedef typename E::size_type matrix_size_type; 1328 matrix_size_type size2 (e ().size2 ()); 1329 for (matrix_size_type j = 0; j < size2; ++ j) { 1330 real_type u = real_type (); 1331 matrix_size_type size1 (e ().size1 ()); 1332 for (matrix_size_type i = 0; i < size1; ++ i) { 1333 real_type v (type_traits<value_type>::norm_1 (e () (i, j))); 1334 u += v; 1335 } 1336 if (u > t) 1337 t = u; 1338 } 1339 return t; 1340 } 1341 }; 1342 1343 template<class M> 1344 struct matrix_norm_frobenius: 1345 public matrix_scalar_real_unary_functor<M> { 1346 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; 1347 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; 1348 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; 1349 1350 template<class E> 1351 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_norm_frobenius1352 result_type apply (const matrix_expression<E> &e) { 1353 real_type t = real_type (); 1354 typedef typename E::size_type matrix_size_type; 1355 matrix_size_type size1 (e ().size1 ()); 1356 for (matrix_size_type i = 0; i < size1; ++ i) { 1357 matrix_size_type size2 (e ().size2 ()); 1358 for (matrix_size_type j = 0; j < size2; ++ j) { 1359 real_type u (type_traits<value_type>::norm_2 (e () (i, j))); 1360 t += u * u; 1361 } 1362 } 1363 return type_traits<real_type>::type_sqrt (t); 1364 } 1365 }; 1366 1367 template<class M> 1368 struct matrix_norm_inf: 1369 public matrix_scalar_real_unary_functor<M> { 1370 typedef typename matrix_scalar_real_unary_functor<M>::value_type value_type; 1371 typedef typename matrix_scalar_real_unary_functor<M>::real_type real_type; 1372 typedef typename matrix_scalar_real_unary_functor<M>::result_type result_type; 1373 1374 template<class E> 1375 static BOOST_UBLAS_INLINE applyboost::numeric::ublas::matrix_norm_inf1376 result_type apply (const matrix_expression<E> &e) { 1377 real_type t = real_type (); 1378 typedef typename E::size_type matrix_size_type; 1379 matrix_size_type size1 (e ().size1 ()); 1380 for (matrix_size_type i = 0; i < size1; ++ i) { 1381 real_type u = real_type (); 1382 matrix_size_type size2 (e ().size2 ()); 1383 for (matrix_size_type j = 0; j < size2; ++ j) { 1384 real_type v (type_traits<value_type>::norm_inf (e () (i, j))); 1385 u += v; 1386 } 1387 if (u > t) 1388 t = u; 1389 } 1390 return t; 1391 } 1392 }; 1393 1394 // forward declaration 1395 template <class Z, class D> struct basic_column_major; 1396 1397 // This functor defines storage layout and it's properties 1398 // matrix (i,j) -> storage [i * size_i + j] 1399 template <class Z, class D> 1400 struct basic_row_major { 1401 typedef Z size_type; 1402 typedef D difference_type; 1403 typedef row_major_tag orientation_category; 1404 typedef basic_column_major<Z,D> transposed_layout; 1405 1406 static 1407 BOOST_UBLAS_INLINE storage_sizeboost::numeric::ublas::basic_row_major1408 size_type storage_size (size_type size_i, size_type size_j) { 1409 // Guard against size_type overflow 1410 BOOST_UBLAS_CHECK (size_j == 0 || size_i <= (std::numeric_limits<size_type>::max) () / size_j, bad_size ()); 1411 return size_i * size_j; 1412 } 1413 1414 // Indexing conversion to storage element 1415 static 1416 BOOST_UBLAS_INLINE elementboost::numeric::ublas::basic_row_major1417 size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { 1418 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 1419 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 1420 detail::ignore_unused_variable_warning(size_i); 1421 // Guard against size_type overflow 1422 BOOST_UBLAS_CHECK (i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ()); 1423 return i * size_j + j; 1424 } 1425 static 1426 BOOST_UBLAS_INLINE addressboost::numeric::ublas::basic_row_major1427 size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { 1428 BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); 1429 BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); 1430 // Guard against size_type overflow - address may be size_j past end of storage 1431 BOOST_UBLAS_CHECK (size_j == 0 || i <= ((std::numeric_limits<size_type>::max) () - j) / size_j, bad_index ()); 1432 detail::ignore_unused_variable_warning(size_i); 1433 return i * size_j + j; 1434 } 1435 1436 // Storage element to index conversion 1437 static 1438 BOOST_UBLAS_INLINE distance_iboost::numeric::ublas::basic_row_major1439 difference_type distance_i (difference_type k, size_type /* size_i */, size_type size_j) { 1440 return size_j != 0 ? k / size_j : 0; 1441 } 1442 static 1443 BOOST_UBLAS_INLINE distance_jboost::numeric::ublas::basic_row_major1444 difference_type distance_j (difference_type k, size_type /* size_i */, size_type /* size_j */) { 1445 return k; 1446 } 1447 static 1448 BOOST_UBLAS_INLINE index_iboost::numeric::ublas::basic_row_major1449 size_type index_i (difference_type k, size_type /* size_i */, size_type size_j) { 1450 return size_j != 0 ? k / size_j : 0; 1451 } 1452 static 1453 BOOST_UBLAS_INLINE index_jboost::numeric::ublas::basic_row_major1454 size_type index_j (difference_type k, size_type /* size_i */, size_type size_j) { 1455 return size_j != 0 ? k % size_j : 0; 1456 } 1457 static 1458 BOOST_UBLAS_INLINE fast_iboost::numeric::ublas::basic_row_major1459 bool fast_i () { 1460 return false; 1461 } 1462 static 1463 BOOST_UBLAS_INLINE fast_jboost::numeric::ublas::basic_row_major1464 bool fast_j () { 1465 return true; 1466 } 1467 1468 // Iterating storage elements 1469 template<class I> 1470 static 1471 BOOST_UBLAS_INLINE increment_iboost::numeric::ublas::basic_row_major1472 void increment_i (I &it, size_type /* size_i */, size_type size_j) { 1473 it += size_j; 1474 } 1475 template<class I> 1476 static 1477 BOOST_UBLAS_INLINE increment_iboost::numeric::ublas::basic_row_major1478 void increment_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { 1479 it += n * size_j; 1480 } 1481 template<class I> 1482 static 1483 BOOST_UBLAS_INLINE decrement_iboost::numeric::ublas::basic_row_major1484 void decrement_i (I &it, size_type /* size_i */, size_type size_j) { 1485 it -= size_j; 1486 } 1487 template<class I> 1488 static 1489 BOOST_UBLAS_INLINE decrement_iboost::numeric::ublas::basic_row_major1490 void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type size_j) { 1491 it -= n * size_j; 1492 } 1493 template<class I> 1494 static 1495 BOOST_UBLAS_INLINE increment_jboost::numeric::ublas::basic_row_major1496 void increment_j (I &it, size_type /* size_i */, size_type /* size_j */) { 1497 ++ it; 1498 } 1499 template<class I> 1500 static 1501 BOOST_UBLAS_INLINE increment_jboost::numeric::ublas::basic_row_major1502 void increment_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 1503 it += n; 1504 } 1505 template<class I> 1506 static 1507 BOOST_UBLAS_INLINE decrement_jboost::numeric::ublas::basic_row_major1508 void decrement_j (I &it, size_type /* size_i */, size_type /* size_j */) { 1509 -- it; 1510 } 1511 template<class I> 1512 static 1513 BOOST_UBLAS_INLINE decrement_jboost::numeric::ublas::basic_row_major1514 void decrement_j (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 1515 it -= n; 1516 } 1517 1518 // Triangular access 1519 static 1520 BOOST_UBLAS_INLINE triangular_sizeboost::numeric::ublas::basic_row_major1521 size_type triangular_size (size_type size_i, size_type size_j) { 1522 size_type size = (std::max) (size_i, size_j); 1523 // Guard against size_type overflow - simplified 1524 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ()); 1525 return ((size + 1) * size) / 2; 1526 } 1527 static 1528 BOOST_UBLAS_INLINE lower_elementboost::numeric::ublas::basic_row_major1529 size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { 1530 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 1531 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 1532 BOOST_UBLAS_CHECK (i >= j, bad_index ()); 1533 detail::ignore_unused_variable_warning(size_i); 1534 detail::ignore_unused_variable_warning(size_j); 1535 // FIXME size_type overflow 1536 // sigma_i (i + 1) = (i + 1) * i / 2 1537 // i = 0 1 2 3, sigma = 0 1 3 6 1538 return ((i + 1) * i) / 2 + j; 1539 } 1540 static 1541 BOOST_UBLAS_INLINE upper_elementboost::numeric::ublas::basic_row_major1542 size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { 1543 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 1544 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 1545 BOOST_UBLAS_CHECK (i <= j, bad_index ()); 1546 // FIXME size_type overflow 1547 // sigma_i (size - i) = size * i - i * (i - 1) / 2 1548 // i = 0 1 2 3, sigma = 0 4 7 9 1549 return (i * (2 * (std::max) (size_i, size_j) - i + 1)) / 2 + j - i; 1550 } 1551 1552 // Major and minor indices 1553 static 1554 BOOST_UBLAS_INLINE index_Mboost::numeric::ublas::basic_row_major1555 size_type index_M (size_type index1, size_type /* index2 */) { 1556 return index1; 1557 } 1558 static 1559 BOOST_UBLAS_INLINE index_mboost::numeric::ublas::basic_row_major1560 size_type index_m (size_type /* index1 */, size_type index2) { 1561 return index2; 1562 } 1563 static 1564 BOOST_UBLAS_INLINE size_Mboost::numeric::ublas::basic_row_major1565 size_type size_M (size_type size_i, size_type /* size_j */) { 1566 return size_i; 1567 } 1568 static 1569 BOOST_UBLAS_INLINE size_mboost::numeric::ublas::basic_row_major1570 size_type size_m (size_type /* size_i */, size_type size_j) { 1571 return size_j; 1572 } 1573 }; 1574 1575 // This functor defines storage layout and it's properties 1576 // matrix (i,j) -> storage [i + j * size_i] 1577 template <class Z, class D> 1578 struct basic_column_major { 1579 typedef Z size_type; 1580 typedef D difference_type; 1581 typedef column_major_tag orientation_category; 1582 typedef basic_row_major<Z,D> transposed_layout; 1583 1584 static 1585 BOOST_UBLAS_INLINE storage_sizeboost::numeric::ublas::basic_column_major1586 size_type storage_size (size_type size_i, size_type size_j) { 1587 // Guard against size_type overflow 1588 BOOST_UBLAS_CHECK (size_i == 0 || size_j <= (std::numeric_limits<size_type>::max) () / size_i, bad_size ()); 1589 return size_i * size_j; 1590 } 1591 1592 // Indexing conversion to storage element 1593 static 1594 BOOST_UBLAS_INLINE elementboost::numeric::ublas::basic_column_major1595 size_type element (size_type i, size_type size_i, size_type j, size_type size_j) { 1596 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 1597 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 1598 detail::ignore_unused_variable_warning(size_j); 1599 // Guard against size_type overflow 1600 BOOST_UBLAS_CHECK (j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ()); 1601 return i + j * size_i; 1602 } 1603 static 1604 BOOST_UBLAS_INLINE addressboost::numeric::ublas::basic_column_major1605 size_type address (size_type i, size_type size_i, size_type j, size_type size_j) { 1606 BOOST_UBLAS_CHECK (i <= size_i, bad_index ()); 1607 BOOST_UBLAS_CHECK (j <= size_j, bad_index ()); 1608 detail::ignore_unused_variable_warning(size_j); 1609 // Guard against size_type overflow - address may be size_i past end of storage 1610 BOOST_UBLAS_CHECK (size_i == 0 || j <= ((std::numeric_limits<size_type>::max) () - i) / size_i, bad_index ()); 1611 return i + j * size_i; 1612 } 1613 1614 // Storage element to index conversion 1615 static 1616 BOOST_UBLAS_INLINE distance_iboost::numeric::ublas::basic_column_major1617 difference_type distance_i (difference_type k, size_type /* size_i */, size_type /* size_j */) { 1618 return k; 1619 } 1620 static 1621 BOOST_UBLAS_INLINE distance_jboost::numeric::ublas::basic_column_major1622 difference_type distance_j (difference_type k, size_type size_i, size_type /* size_j */) { 1623 return size_i != 0 ? k / size_i : 0; 1624 } 1625 static 1626 BOOST_UBLAS_INLINE index_iboost::numeric::ublas::basic_column_major1627 size_type index_i (difference_type k, size_type size_i, size_type /* size_j */) { 1628 return size_i != 0 ? k % size_i : 0; 1629 } 1630 static 1631 BOOST_UBLAS_INLINE index_jboost::numeric::ublas::basic_column_major1632 size_type index_j (difference_type k, size_type size_i, size_type /* size_j */) { 1633 return size_i != 0 ? k / size_i : 0; 1634 } 1635 static 1636 BOOST_UBLAS_INLINE fast_iboost::numeric::ublas::basic_column_major1637 bool fast_i () { 1638 return true; 1639 } 1640 static 1641 BOOST_UBLAS_INLINE fast_jboost::numeric::ublas::basic_column_major1642 bool fast_j () { 1643 return false; 1644 } 1645 1646 // Iterating 1647 template<class I> 1648 static 1649 BOOST_UBLAS_INLINE increment_iboost::numeric::ublas::basic_column_major1650 void increment_i (I &it, size_type /* size_i */, size_type /* size_j */) { 1651 ++ it; 1652 } 1653 template<class I> 1654 static 1655 BOOST_UBLAS_INLINE increment_iboost::numeric::ublas::basic_column_major1656 void increment_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 1657 it += n; 1658 } 1659 template<class I> 1660 static 1661 BOOST_UBLAS_INLINE decrement_iboost::numeric::ublas::basic_column_major1662 void decrement_i (I &it, size_type /* size_i */, size_type /* size_j */) { 1663 -- it; 1664 } 1665 template<class I> 1666 static 1667 BOOST_UBLAS_INLINE decrement_iboost::numeric::ublas::basic_column_major1668 void decrement_i (I &it, difference_type n, size_type /* size_i */, size_type /* size_j */) { 1669 it -= n; 1670 } 1671 template<class I> 1672 static 1673 BOOST_UBLAS_INLINE increment_jboost::numeric::ublas::basic_column_major1674 void increment_j (I &it, size_type size_i, size_type /* size_j */) { 1675 it += size_i; 1676 } 1677 template<class I> 1678 static 1679 BOOST_UBLAS_INLINE increment_jboost::numeric::ublas::basic_column_major1680 void increment_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { 1681 it += n * size_i; 1682 } 1683 template<class I> 1684 static 1685 BOOST_UBLAS_INLINE decrement_jboost::numeric::ublas::basic_column_major1686 void decrement_j (I &it, size_type size_i, size_type /* size_j */) { 1687 it -= size_i; 1688 } 1689 template<class I> 1690 static 1691 BOOST_UBLAS_INLINE decrement_jboost::numeric::ublas::basic_column_major1692 void decrement_j (I &it, difference_type n, size_type size_i, size_type /* size_j */) { 1693 it -= n* size_i; 1694 } 1695 1696 // Triangular access 1697 static 1698 BOOST_UBLAS_INLINE triangular_sizeboost::numeric::ublas::basic_column_major1699 size_type triangular_size (size_type size_i, size_type size_j) { 1700 size_type size = (std::max) (size_i, size_j); 1701 // Guard against size_type overflow - simplified 1702 BOOST_UBLAS_CHECK (size == 0 || size / 2 < (std::numeric_limits<size_type>::max) () / size /* +1/2 */, bad_size ()); 1703 return ((size + 1) * size) / 2; 1704 } 1705 static 1706 BOOST_UBLAS_INLINE lower_elementboost::numeric::ublas::basic_column_major1707 size_type lower_element (size_type i, size_type size_i, size_type j, size_type size_j) { 1708 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 1709 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 1710 BOOST_UBLAS_CHECK (i >= j, bad_index ()); 1711 // FIXME size_type overflow 1712 // sigma_j (size - j) = size * j - j * (j - 1) / 2 1713 // j = 0 1 2 3, sigma = 0 4 7 9 1714 return i - j + (j * (2 * (std::max) (size_i, size_j) - j + 1)) / 2; 1715 } 1716 static 1717 BOOST_UBLAS_INLINE upper_elementboost::numeric::ublas::basic_column_major1718 size_type upper_element (size_type i, size_type size_i, size_type j, size_type size_j) { 1719 BOOST_UBLAS_CHECK (i < size_i, bad_index ()); 1720 BOOST_UBLAS_CHECK (j < size_j, bad_index ()); 1721 BOOST_UBLAS_CHECK (i <= j, bad_index ()); 1722 boost::ignore_unused(size_i, size_j); 1723 // FIXME size_type overflow 1724 // sigma_j (j + 1) = (j + 1) * j / 2 1725 // j = 0 1 2 3, sigma = 0 1 3 6 1726 return i + ((j + 1) * j) / 2; 1727 } 1728 1729 // Major and minor indices 1730 static 1731 BOOST_UBLAS_INLINE index_Mboost::numeric::ublas::basic_column_major1732 size_type index_M (size_type /* index1 */, size_type index2) { 1733 return index2; 1734 } 1735 static 1736 BOOST_UBLAS_INLINE index_mboost::numeric::ublas::basic_column_major1737 size_type index_m (size_type index1, size_type /* index2 */) { 1738 return index1; 1739 } 1740 static 1741 BOOST_UBLAS_INLINE size_Mboost::numeric::ublas::basic_column_major1742 size_type size_M (size_type /* size_i */, size_type size_j) { 1743 return size_j; 1744 } 1745 static 1746 BOOST_UBLAS_INLINE size_mboost::numeric::ublas::basic_column_major1747 size_type size_m (size_type size_i, size_type /* size_j */) { 1748 return size_i; 1749 } 1750 }; 1751 1752 1753 template <class Z> 1754 struct basic_full { 1755 typedef Z size_type; 1756 1757 template<class L> 1758 static 1759 BOOST_UBLAS_INLINE packed_sizeboost::numeric::ublas::basic_full1760 size_type packed_size (L, size_type size_i, size_type size_j) { 1761 return L::storage_size (size_i, size_j); 1762 } 1763 1764 static 1765 BOOST_UBLAS_INLINE zeroboost::numeric::ublas::basic_full1766 bool zero (size_type /* i */, size_type /* j */) { 1767 return false; 1768 } 1769 static 1770 BOOST_UBLAS_INLINE oneboost::numeric::ublas::basic_full1771 bool one (size_type /* i */, size_type /* j */) { 1772 return false; 1773 } 1774 static 1775 BOOST_UBLAS_INLINE otherboost::numeric::ublas::basic_full1776 bool other (size_type /* i */, size_type /* j */) { 1777 return true; 1778 } 1779 // FIXME: this should not be used at all 1780 static 1781 BOOST_UBLAS_INLINE restrict1boost::numeric::ublas::basic_full1782 size_type restrict1 (size_type i, size_type /* j */) { 1783 return i; 1784 } 1785 static 1786 BOOST_UBLAS_INLINE restrict2boost::numeric::ublas::basic_full1787 size_type restrict2 (size_type /* i */, size_type j) { 1788 return j; 1789 } 1790 static 1791 BOOST_UBLAS_INLINE mutable_restrict1boost::numeric::ublas::basic_full1792 size_type mutable_restrict1 (size_type i, size_type /* j */) { 1793 return i; 1794 } 1795 static 1796 BOOST_UBLAS_INLINE mutable_restrict2boost::numeric::ublas::basic_full1797 size_type mutable_restrict2 (size_type /* i */, size_type j) { 1798 return j; 1799 } 1800 }; 1801 1802 namespace detail { 1803 template < class L > 1804 struct transposed_structure { 1805 typedef typename L::size_type size_type; 1806 1807 template<class LAYOUT> 1808 static 1809 BOOST_UBLAS_INLINE packed_sizeboost::numeric::ublas::detail::transposed_structure1810 size_type packed_size (LAYOUT l, size_type size_i, size_type size_j) { 1811 return L::packed_size(l, size_j, size_i); 1812 } 1813 1814 static 1815 BOOST_UBLAS_INLINE zeroboost::numeric::ublas::detail::transposed_structure1816 bool zero (size_type i, size_type j) { 1817 return L::zero(j, i); 1818 } 1819 static 1820 BOOST_UBLAS_INLINE oneboost::numeric::ublas::detail::transposed_structure1821 bool one (size_type i, size_type j) { 1822 return L::one(j, i); 1823 } 1824 static 1825 BOOST_UBLAS_INLINE otherboost::numeric::ublas::detail::transposed_structure1826 bool other (size_type i, size_type j) { 1827 return L::other(j, i); 1828 } 1829 template<class LAYOUT> 1830 static 1831 BOOST_UBLAS_INLINE elementboost::numeric::ublas::detail::transposed_structure1832 size_type element (LAYOUT /* l */, size_type i, size_type size_i, size_type j, size_type size_j) { 1833 return L::element(typename LAYOUT::transposed_layout(), j, size_j, i, size_i); 1834 } 1835 1836 static 1837 BOOST_UBLAS_INLINE restrict1boost::numeric::ublas::detail::transposed_structure1838 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 1839 return L::restrict2(j, i, size2, size1); 1840 } 1841 static 1842 BOOST_UBLAS_INLINE restrict2boost::numeric::ublas::detail::transposed_structure1843 size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { 1844 return L::restrict1(j, i, size2, size1); 1845 } 1846 static 1847 BOOST_UBLAS_INLINE mutable_restrict1boost::numeric::ublas::detail::transposed_structure1848 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 1849 return L::mutable_restrict2(j, i, size2, size1); 1850 } 1851 static 1852 BOOST_UBLAS_INLINE mutable_restrict2boost::numeric::ublas::detail::transposed_structure1853 size_type mutable_restrict2 (size_type i, size_type j, size_type size1, size_type size2) { 1854 return L::mutable_restrict1(j, i, size2, size1); 1855 } 1856 1857 static 1858 BOOST_UBLAS_INLINE global_restrict1boost::numeric::ublas::detail::transposed_structure1859 size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { 1860 return L::global_restrict2(index2, size2, index1, size1); 1861 } 1862 static 1863 BOOST_UBLAS_INLINE global_restrict2boost::numeric::ublas::detail::transposed_structure1864 size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { 1865 return L::global_restrict1(index2, size2, index1, size1); 1866 } 1867 static 1868 BOOST_UBLAS_INLINE global_mutable_restrict1boost::numeric::ublas::detail::transposed_structure1869 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { 1870 return L::global_mutable_restrict2(index2, size2, index1, size1); 1871 } 1872 static 1873 BOOST_UBLAS_INLINE global_mutable_restrict2boost::numeric::ublas::detail::transposed_structure1874 size_type global_mutable_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { 1875 return L::global_mutable_restrict1(index2, size2, index1, size1); 1876 } 1877 }; 1878 } 1879 1880 template <class Z> 1881 struct basic_lower { 1882 typedef Z size_type; 1883 typedef lower_tag triangular_type; 1884 1885 template<class L> 1886 static 1887 BOOST_UBLAS_INLINE packed_sizeboost::numeric::ublas::basic_lower1888 size_type packed_size (L, size_type size_i, size_type size_j) { 1889 return L::triangular_size (size_i, size_j); 1890 } 1891 1892 static 1893 BOOST_UBLAS_INLINE zeroboost::numeric::ublas::basic_lower1894 bool zero (size_type i, size_type j) { 1895 return j > i; 1896 } 1897 static 1898 BOOST_UBLAS_INLINE oneboost::numeric::ublas::basic_lower1899 bool one (size_type /* i */, size_type /* j */) { 1900 return false; 1901 } 1902 static 1903 BOOST_UBLAS_INLINE otherboost::numeric::ublas::basic_lower1904 bool other (size_type i, size_type j) { 1905 return j <= i; 1906 } 1907 template<class L> 1908 static 1909 BOOST_UBLAS_INLINE elementboost::numeric::ublas::basic_lower1910 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { 1911 return L::lower_element (i, size_i, j, size_j); 1912 } 1913 1914 // return nearest valid index in column j 1915 static 1916 BOOST_UBLAS_INLINE restrict1boost::numeric::ublas::basic_lower1917 size_type restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { 1918 return (std::max)(j, (std::min) (size1, i)); 1919 } 1920 // return nearest valid index in row i 1921 static 1922 BOOST_UBLAS_INLINE restrict2boost::numeric::ublas::basic_lower1923 size_type restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { 1924 return (std::max)(size_type(0), (std::min) (i+1, j)); 1925 } 1926 // return nearest valid mutable index in column j 1927 static 1928 BOOST_UBLAS_INLINE mutable_restrict1boost::numeric::ublas::basic_lower1929 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { 1930 return (std::max)(j, (std::min) (size1, i)); 1931 } 1932 // return nearest valid mutable index in row i 1933 static 1934 BOOST_UBLAS_INLINE mutable_restrict2boost::numeric::ublas::basic_lower1935 size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { 1936 return (std::max)(size_type(0), (std::min) (i+1, j)); 1937 } 1938 1939 // return an index between the first and (1+last) filled row 1940 static 1941 BOOST_UBLAS_INLINE global_restrict1boost::numeric::ublas::basic_lower1942 size_type global_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { 1943 return (std::max)(size_type(0), (std::min)(size1, index1) ); 1944 } 1945 // return an index between the first and (1+last) filled column 1946 static 1947 BOOST_UBLAS_INLINE global_restrict2boost::numeric::ublas::basic_lower1948 size_type global_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { 1949 return (std::max)(size_type(0), (std::min)(size2, index2) ); 1950 } 1951 1952 // return an index between the first and (1+last) filled mutable row 1953 static 1954 BOOST_UBLAS_INLINE global_mutable_restrict1boost::numeric::ublas::basic_lower1955 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { 1956 return (std::max)(size_type(0), (std::min)(size1, index1) ); 1957 } 1958 // return an index between the first and (1+last) filled mutable column 1959 static 1960 BOOST_UBLAS_INLINE global_mutable_restrict2boost::numeric::ublas::basic_lower1961 size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { 1962 return (std::max)(size_type(0), (std::min)(size2, index2) ); 1963 } 1964 }; 1965 1966 // the first row only contains a single 1. Thus it is not stored. 1967 template <class Z> 1968 struct basic_unit_lower : public basic_lower<Z> { 1969 typedef Z size_type; 1970 typedef unit_lower_tag triangular_type; 1971 1972 template<class L> 1973 static 1974 BOOST_UBLAS_INLINE packed_sizeboost::numeric::ublas::basic_unit_lower1975 size_type packed_size (L, size_type size_i, size_type size_j) { 1976 // Zero size strict triangles are bad at this point 1977 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); 1978 return L::triangular_size (size_i - 1, size_j - 1); 1979 } 1980 1981 static 1982 BOOST_UBLAS_INLINE oneboost::numeric::ublas::basic_unit_lower1983 bool one (size_type i, size_type j) { 1984 return j == i; 1985 } 1986 static 1987 BOOST_UBLAS_INLINE otherboost::numeric::ublas::basic_unit_lower1988 bool other (size_type i, size_type j) { 1989 return j < i; 1990 } 1991 template<class L> 1992 static 1993 BOOST_UBLAS_INLINE elementboost::numeric::ublas::basic_unit_lower1994 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { 1995 // Zero size strict triangles are bad at this point 1996 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); 1997 return L::lower_element (i-1, size_i - 1, j, size_j - 1); 1998 } 1999 2000 static 2001 BOOST_UBLAS_INLINE mutable_restrict1boost::numeric::ublas::basic_unit_lower2002 size_type mutable_restrict1 (size_type i, size_type j, size_type size1, size_type /* size2 */) { 2003 return (std::max)(j+1, (std::min) (size1, i)); 2004 } 2005 static 2006 BOOST_UBLAS_INLINE mutable_restrict2boost::numeric::ublas::basic_unit_lower2007 size_type mutable_restrict2 (size_type i, size_type j, size_type /* size1 */, size_type /* size2 */) { 2008 return (std::max)(size_type(0), (std::min) (i, j)); 2009 } 2010 2011 // return an index between the first and (1+last) filled mutable row 2012 static 2013 BOOST_UBLAS_INLINE global_mutable_restrict1boost::numeric::ublas::basic_unit_lower2014 size_type global_mutable_restrict1 (size_type index1, size_type size1, size_type /* index2 */, size_type /* size2 */) { 2015 return (std::max)(size_type(1), (std::min)(size1, index1) ); 2016 } 2017 // return an index between the first and (1+last) filled mutable column 2018 static 2019 BOOST_UBLAS_INLINE global_mutable_restrict2boost::numeric::ublas::basic_unit_lower2020 size_type global_mutable_restrict2 (size_type /* index1 */, size_type /* size1 */, size_type index2, size_type size2) { 2021 BOOST_UBLAS_CHECK( size2 >= 1 , external_logic() ); 2022 return (std::max)(size_type(0), (std::min)(size2-1, index2) ); 2023 } 2024 }; 2025 2026 // the first row only contains no element. Thus it is not stored. 2027 template <class Z> 2028 struct basic_strict_lower : public basic_unit_lower<Z> { 2029 typedef Z size_type; 2030 typedef strict_lower_tag triangular_type; 2031 2032 template<class L> 2033 static 2034 BOOST_UBLAS_INLINE packed_sizeboost::numeric::ublas::basic_strict_lower2035 size_type packed_size (L, size_type size_i, size_type size_j) { 2036 // Zero size strict triangles are bad at this point 2037 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0, bad_index ()); 2038 return L::triangular_size (size_i - 1, size_j - 1); 2039 } 2040 2041 static 2042 BOOST_UBLAS_INLINE zeroboost::numeric::ublas::basic_strict_lower2043 bool zero (size_type i, size_type j) { 2044 return j >= i; 2045 } 2046 static 2047 BOOST_UBLAS_INLINE oneboost::numeric::ublas::basic_strict_lower2048 bool one (size_type /*i*/, size_type /*j*/) { 2049 return false; 2050 } 2051 static 2052 BOOST_UBLAS_INLINE otherboost::numeric::ublas::basic_strict_lower2053 bool other (size_type i, size_type j) { 2054 return j < i; 2055 } 2056 template<class L> 2057 static 2058 BOOST_UBLAS_INLINE elementboost::numeric::ublas::basic_strict_lower2059 size_type element (L, size_type i, size_type size_i, size_type j, size_type size_j) { 2060 // Zero size strict triangles are bad at this point 2061 BOOST_UBLAS_CHECK (size_i != 0 && size_j != 0 && i != 0, bad_index ()); 2062 return L::lower_element (i-1, size_i - 1, j, size_j - 1); 2063 } 2064 2065 static 2066 BOOST_UBLAS_INLINE restrict1boost::numeric::ublas::basic_strict_lower2067 size_type restrict1 (size_type i, size_type j, size_type size1, size_type size2) { 2068 return basic_unit_lower<Z>::mutable_restrict1(i, j, size1, size2); 2069 } 2070 static 2071 BOOST_UBLAS_INLINE restrict2boost::numeric::ublas::basic_strict_lower2072 size_type restrict2 (size_type i, size_type j, size_type size1, size_type size2) { 2073 return basic_unit_lower<Z>::mutable_restrict2(i, j, size1, size2); 2074 } 2075 2076 // return an index between the first and (1+last) filled row 2077 static 2078 BOOST_UBLAS_INLINE global_restrict1boost::numeric::ublas::basic_strict_lower2079 size_type global_restrict1 (size_type index1, size_type size1, size_type index2, size_type size2) { 2080 return basic_unit_lower<Z>::global_mutable_restrict1(index1, size1, index2, size2); 2081 } 2082 // return an index between the first and (1+last) filled column 2083 static 2084 BOOST_UBLAS_INLINE global_restrict2boost::numeric::ublas::basic_strict_lower2085 size_type global_restrict2 (size_type index1, size_type size1, size_type index2, size_type size2) { 2086 return basic_unit_lower<Z>::global_mutable_restrict2(index1, size1, index2, size2); 2087 } 2088 }; 2089 2090 2091 template <class Z> 2092 struct basic_upper : public detail::transposed_structure<basic_lower<Z> > 2093 { 2094 typedef upper_tag triangular_type; 2095 }; 2096 2097 template <class Z> 2098 struct basic_unit_upper : public detail::transposed_structure<basic_unit_lower<Z> > 2099 { 2100 typedef unit_upper_tag triangular_type; 2101 }; 2102 2103 template <class Z> 2104 struct basic_strict_upper : public detail::transposed_structure<basic_strict_lower<Z> > 2105 { 2106 typedef strict_upper_tag triangular_type; 2107 }; 2108 2109 2110 }}} 2111 2112 #endif 2113