SquareMatrix.h
Go to the documentation of this file.00001 // -*- C++ -*- 00002 00008 #if !defined(__SquareMatrix_h__) 00009 #define __SquareMatrix_h__ 00010 00011 #if defined(DEBUG_ads) && !defined(DEBUG_SquareMatrix) 00012 #define DEBUG_SquareMatrix 00013 #endif 00014 00015 #include "TensorTypes.h" 00016 00017 #include "../array/FixedArray.h" 00018 00019 #include "../../third-party/loki/static_check.h" 00020 00021 #include <iosfwd> 00022 #include <algorithm> 00023 #include <functional> 00024 #include <numeric> 00025 00026 #include <cassert> 00027 #include <cmath> 00028 00029 BEGIN_NAMESPACE_ADS 00030 00031 // 00032 // NxN matrices. 00033 // 00034 00036 00040 template<int N, typename T = double> 00041 class SquareMatrix : 00042 public TensorTypes<T> { 00043 private: 00044 00045 // 00046 // Private types. 00047 // 00048 00049 typedef TensorTypes<T> base_type; 00050 00051 public: 00052 00053 // 00054 // public types 00055 // 00056 00058 typedef typename base_type::value_type value_type; 00060 typedef typename base_type::pointer pointer; 00062 typedef typename base_type::const_pointer const_pointer; 00064 typedef typename base_type::iterator iterator; 00066 typedef typename base_type::const_iterator const_iterator; 00068 typedef typename base_type::reference reference; 00070 typedef typename base_type::const_reference const_reference; 00072 typedef typename base_type::size_type size_type; 00074 typedef typename base_type::difference_type difference_type; 00076 typedef typename base_type::index_type index_type; 00077 00078 protected: 00079 00080 // 00081 // Data 00082 // 00083 00085 value_type _elem[N][N]; 00086 00087 public: 00088 00089 // 00090 // Constructors 00091 // 00092 00094 SquareMatrix() { 00095 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00096 } 00097 00099 ~SquareMatrix() { 00100 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00101 } 00102 00104 SquareMatrix(const SquareMatrix& x); 00105 00107 template<typename T2> 00108 SquareMatrix(const SquareMatrix<N,T2>& x); 00109 00111 SquareMatrix(const_pointer x); 00112 00114 SquareMatrix(const value_type x); 00115 00116 // 00117 // Assignment operators 00118 // 00119 00121 SquareMatrix& 00122 operator=(const SquareMatrix& x); 00123 00125 SquareMatrix& 00126 operator=(const value_type x); 00127 00129 template<typename T2> 00130 SquareMatrix& 00131 operator=(const SquareMatrix<N,T2>& x); 00132 00133 // 00134 // Accessors and Manipulators 00135 // 00136 00138 iterator 00139 begin() { 00140 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00141 return &(_elem[0][0]); 00142 } 00143 00145 iterator 00146 end() { 00147 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00148 return begin() + N * N; 00149 } 00150 00152 const_iterator 00153 begin() const { 00154 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00155 return &(_elem[0][0]); 00156 } 00157 00159 const_iterator 00160 end() const { 00161 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00162 return begin() + N * N; 00163 } 00164 00166 pointer 00167 data() { 00168 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00169 return &(_elem[0][0]); 00170 } 00171 00173 const_pointer 00174 data() const { 00175 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00176 return &(_elem[0][0]); 00177 } 00178 00180 size_type 00181 size() const { 00182 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00183 return N * N; 00184 } 00185 00187 value_type 00188 operator()(const index_type i) const { 00189 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00190 #ifdef DEBUG_SquareMatrix 00191 assert(0 <= i && i < size()); 00192 #endif 00193 return *(begin() + i); 00194 } 00195 00197 reference 00198 operator()(const index_type i) { 00199 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00200 #ifdef DEBUG_SquareMatrix 00201 assert(0 <= i && i < size()); 00202 #endif 00203 return *(begin() + i); 00204 } 00205 00207 value_type 00208 operator[](const index_type i) const { 00209 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00210 #ifdef DEBUG_SquareMatrix 00211 assert(0 <= i && i < size()); 00212 #endif 00213 return *(begin() + i); 00214 } 00215 00217 value_type& 00218 operator[](const index_type i) { 00219 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00220 #ifdef DEBUG_SquareMatrix 00221 assert(0 <= i && i < size()); 00222 #endif 00223 return *(begin() + i); 00224 } 00225 00227 value_type 00228 operator()(const index_type i, const index_type j) const { 00229 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00230 #ifdef DEBUG_SquareMatrix 00231 assert(0 <= i && i < N && 0 <= j && j < N); 00232 #endif 00233 return _elem[i][j]; 00234 } 00235 00237 reference 00238 operator()(const index_type i, const index_type j) { 00239 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00240 #ifdef DEBUG_SquareMatrix 00241 assert(0 <= i && i < N && 0 <= j && j < N); 00242 #endif 00243 return _elem[i][j]; 00244 } 00245 00247 void 00248 get(pointer x) const; 00249 00251 void 00252 set(const_pointer x); 00253 00255 00258 void 00259 getMinor(const int i, const int j, SquareMatrix<N-1,T>& minor) const; 00260 00262 void 00263 negate(); 00264 00266 void 00267 transpose(); 00268 00269 // 00270 // Assignment operators with scalar operand. 00271 // 00272 00274 SquareMatrix& 00275 operator+=(const value_type x); 00276 00278 SquareMatrix& 00279 operator-=(const value_type x); 00280 00282 SquareMatrix& 00283 operator*=(const value_type x); 00284 00286 SquareMatrix& 00287 operator/=(const value_type x); 00288 00290 SquareMatrix& 00291 operator%=(const value_type x); 00292 00293 // 00294 // Assignment operators with SquareMatrix operand 00295 // 00296 00298 template<typename T2> 00299 SquareMatrix& 00300 operator+=(const SquareMatrix<N,T2> & x); 00301 00303 template<typename T2> 00304 SquareMatrix& 00305 operator-=(const SquareMatrix<N,T2> & x); 00306 00308 template<typename T2> 00309 SquareMatrix& 00310 operator*=(const SquareMatrix<N,T2> & x); 00311 00312 // 00313 // Unary operators 00314 // 00315 00317 SquareMatrix 00318 operator+() { 00319 return *this; 00320 } 00321 00323 SquareMatrix 00324 operator-(); 00325 00326 }; 00327 00328 // 00329 // Binary operators 00330 // 00331 00333 template<typename T, int N> 00334 inline 00335 SquareMatrix<N,T> 00336 operator+(const SquareMatrix<N,T>& m, const T x) { 00337 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00338 SquareMatrix<N,T> result(m); 00339 result += x; 00340 return result; 00341 } 00342 00344 template<typename T, int N> 00345 inline 00346 SquareMatrix<N,T> 00347 operator+(const T x, const SquareMatrix<N,T>& m) { 00348 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00349 return m + x; 00350 } 00351 00353 template<typename T, int N> 00354 inline 00355 SquareMatrix<N,T> 00356 operator+(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y) { 00357 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00358 SquareMatrix<N,T> result(x); 00359 result += y; 00360 return result; 00361 } 00362 00364 template<typename T, int N> 00365 inline 00366 SquareMatrix<N,T> 00367 operator-(const SquareMatrix<N,T>& m, const T x) { 00368 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00369 SquareMatrix<N,T> result(m); 00370 result -= x; 00371 return result; 00372 } 00373 00375 template<typename T, int N> 00376 inline 00377 SquareMatrix<N,T> 00378 operator-(const T x, const SquareMatrix<N,T>& m) { 00379 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00380 SquareMatrix<N,T> result(x); 00381 result -= m; 00382 return result; 00383 } 00384 00386 template<typename T, int N> 00387 inline 00388 SquareMatrix<N,T> 00389 operator-(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y) { 00390 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00391 SquareMatrix<N,T> result(x); 00392 result -= y; 00393 return result; 00394 } 00395 00397 template<int N, typename T> 00398 inline 00399 SquareMatrix<N,T> 00400 operator*(const SquareMatrix<N,T>& m, const T x) { 00401 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00402 SquareMatrix<N,T> result(m); 00403 result *= x; 00404 return result; 00405 } 00406 00408 template<int N, typename T> 00409 inline 00410 SquareMatrix<N,T> 00411 operator*(const T x, const SquareMatrix<N,T>& m) { 00412 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00413 return m * x; 00414 } 00415 00417 template<int N, typename T> 00418 SquareMatrix<N,T> 00419 operator*(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y); 00420 00422 template<int N, typename T> 00423 inline 00424 SquareMatrix<N,T> 00425 operator/(const SquareMatrix<N,T>& m, const T x) { 00426 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00427 #ifdef DEBUG_SquareMatrix 00428 assert(x != 0); 00429 #endif 00430 SquareMatrix<N,T> result(m); 00431 result /= x; 00432 return result; 00433 } 00434 00436 template<int N, typename T> 00437 SquareMatrix<N,T> 00438 operator/(const T x, const SquareMatrix<N,T>& m); 00439 00440 // 00441 // Math operators. 00442 // 00443 00445 template<int N, typename T> 00446 inline 00447 T 00448 computeSum(const SquareMatrix<N,T>& x) { 00449 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00450 return std::accumulate(x.begin(), x.end(), T(0)); 00451 } 00452 00454 template<int N, typename T> 00455 inline 00456 T 00457 computeProduct(const SquareMatrix<N,T>& x) { 00458 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00459 return std::accumulate(x.begin(), x.end(), T(1), std::multiplies<T>()); 00460 } 00461 00463 template<int N, typename T> 00464 inline 00465 T 00466 computeMinimum(const SquareMatrix<N,T>& x) { 00467 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00468 return *std::min_element(x.begin(), x.end()); 00469 } 00470 00472 template<int N, typename T> 00473 inline 00474 T 00475 computeMaximum(const SquareMatrix<N,T>& x) { 00476 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00477 return *std::max_element(x.begin(), x.end()); 00478 } 00479 00481 template<int N, typename T> 00482 T 00483 computeDeterminant(const SquareMatrix<N,T>& x); 00484 00486 template<int N, typename T> 00487 T 00488 computeTrace(const SquareMatrix<N,T>& x); 00489 00491 template<int N, typename T> 00492 SquareMatrix<N,T> 00493 computeTranspose(const SquareMatrix<N,T>& x); 00494 00496 template<int N, typename T> 00497 SquareMatrix<N,T> 00498 computeInverse(const SquareMatrix<N,T>& x); 00499 00501 template<int N, typename T> 00502 void 00503 computeInverse(const SquareMatrix<N,T>& x, SquareMatrix<N,T>* y); 00504 00506 template<int N, typename T> 00507 void 00508 computeScaledInverse(const SquareMatrix<N,T>& x, SquareMatrix<N,T>* si); 00509 00511 template<int N, typename T> 00512 SquareMatrix<N,T> 00513 computeScaledInverse(const SquareMatrix<N,T>& x); 00514 00516 template<int N, typename T> 00517 T 00518 computeFrobeniusNormSquared(const SquareMatrix<N,T>& x); 00519 00521 template<int N, typename T> 00522 inline 00523 T 00524 computeFrobeniusNorm(const SquareMatrix<N,T>& x) { 00525 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00526 return std::sqrt(computeFrobeniusNormSquared(x)); 00527 } 00528 00530 template<int N, typename T> 00531 inline 00532 T 00533 computeInnerProduct(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y) { 00534 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00535 T result = 0; 00536 const typename SquareMatrix<N,T>::const_iterator finish = x.end(); 00537 for (typename SquareMatrix<N,T>::const_iterator i = x.begin(), 00538 j = y.begin(); i != finish; ++i, ++j) { 00539 result += *i * *j; 00540 } 00541 return result; 00542 } 00543 00544 00546 template<int N, typename T> 00547 inline 00548 void 00549 computeProduct(const SquareMatrix<N,T>& m, const FixedArray<N,T>& v, 00550 FixedArray<N,T>* x) { 00551 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00552 *x = 0; 00553 typename SquareMatrix<N,T>::const_iterator mi = x->begin(); 00554 typename FixedArray<N,T>::const_iterator vi; 00555 const typename FixedArray<N,T>::const_iterator vi_end = v.end(); 00556 // Loop over the rows. 00557 for (int m = 0; m != N; ++m) { 00558 // Loop over the columns. 00559 for (vi = v.begin(); vi != vi_end; ++vi, ++mi) { 00560 (*x)[m] += *mi * *vi; 00561 } 00562 } 00563 } 00564 00565 00566 // 00567 // Equality 00568 // 00569 00571 template<int N, typename T1, typename T2> 00572 inline 00573 bool 00574 operator==(const SquareMatrix<N,T1>& a, const SquareMatrix<N,T2>& b) { 00575 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00576 return std::equal(a.begin(), a.end(), b.begin()); 00577 } 00578 00580 template<int N, typename T1, typename T2> 00581 inline 00582 bool 00583 operator!=(const SquareMatrix<N,T1>& a, const SquareMatrix<N,T2>& b) { 00584 LOKI_STATIC_CHECK(N > 3, N_greater_than_3); 00585 return !(a == b); 00586 } 00587 00588 // 00589 // I/O 00590 // 00591 00593 template<int N, typename T> 00594 std::ostream& 00595 operator<<(std::ostream& out, const SquareMatrix<N,T>& x); 00596 00598 template<int N, typename T> 00599 std::istream& 00600 operator>>(std::istream& in, SquareMatrix<N,T>& x); 00601 00602 00603 00604 00605 //---------------------------------------------------------------------------- 00606 // 1x1 matrices. 00607 //---------------------------------------------------------------------------- 00608 00610 00613 template<typename T> 00614 class SquareMatrix<1,T> : 00615 public TensorTypes<T> { 00616 private: 00617 00618 // 00619 // private types 00620 // 00621 00622 typedef TensorTypes<T> base_type; 00623 00624 public: 00625 00626 // 00627 // public types 00628 // 00629 00631 typedef typename base_type::value_type value_type; 00633 typedef typename base_type::pointer pointer; 00635 typedef typename base_type::const_pointer const_pointer; 00637 typedef typename base_type::iterator iterator; 00639 typedef typename base_type::const_iterator const_iterator; 00641 typedef typename base_type::reference reference; 00643 typedef typename base_type::const_reference const_reference; 00645 typedef typename base_type::size_type size_type; 00647 typedef typename base_type::difference_type difference_type; 00649 typedef typename base_type::index_type index_type; 00650 00651 protected: 00652 00653 // 00654 // Data 00655 // 00656 00658 value_type _elem; 00659 00660 public: 00661 00662 // 00663 // Constructors 00664 // 00665 00667 SquareMatrix() 00668 {} 00669 00671 ~SquareMatrix() 00672 {} 00673 00675 SquareMatrix(const SquareMatrix& x); 00676 00678 template<typename T2> 00679 SquareMatrix(const SquareMatrix<1,T2>& x); 00680 00682 SquareMatrix(const value_type e); 00683 00685 SquareMatrix(const_pointer x); 00686 00687 // 00688 // Assignment operators 00689 // 00690 00692 SquareMatrix& 00693 operator=(const SquareMatrix& x); 00694 00696 SquareMatrix& 00697 operator=(const value_type x); 00698 00700 template<typename T2> 00701 SquareMatrix& 00702 operator=(const SquareMatrix<1,T2>& x); 00703 00704 // 00705 // Accessors and Manipulators 00706 // 00707 00709 iterator 00710 begin() { 00711 return &_elem; 00712 } 00713 00715 iterator 00716 end() { 00717 return &_elem + 1; 00718 } 00719 00721 const_iterator 00722 begin() const { 00723 return &_elem; 00724 } 00725 00727 const_iterator 00728 end() const { 00729 return &_elem + 1; 00730 } 00731 00733 pointer 00734 data() { 00735 return &_elem; 00736 } 00737 00739 const_pointer 00740 data() const { 00741 return &_elem; 00742 } 00743 00745 size_type 00746 size() const { 00747 return 1; 00748 } 00749 00751 value_type 00752 operator()(const index_type i) const { 00753 #ifdef DEBUG_SquareMatrix 00754 assert(i == 0); 00755 #endif 00756 return _elem; 00757 } 00758 00760 reference 00761 operator()(const index_type i) { 00762 #ifdef DEBUG_SquareMatrix 00763 assert(i == 0); 00764 #endif 00765 return _elem; 00766 } 00767 00769 value_type 00770 operator[](const index_type i) const { 00771 #ifdef DEBUG_SquareMatrix 00772 assert(i == 0); 00773 #endif 00774 return _elem; 00775 } 00776 00778 value_type& 00779 operator[](const index_type i) { 00780 #ifdef DEBUG_SquareMatrix 00781 assert(i == 0); 00782 #endif 00783 return _elem; 00784 } 00785 00787 value_type 00788 operator()(const index_type i, const index_type j) const { 00789 #ifdef DEBUG_SquareMatrix 00790 assert(i == 0 && j == 0); 00791 #endif 00792 return _elem; 00793 } 00794 00796 reference 00797 operator()(const index_type i, const index_type j) { 00798 #ifdef DEBUG_SquareMatrix 00799 assert(i == 0 && j == 0); 00800 #endif 00801 return _elem; 00802 } 00803 00805 void 00806 get(pointer x) const; 00807 00809 void 00810 set(const_pointer x); 00811 00813 void 00814 get(reference e) const; 00815 00817 void 00818 set(const value_type e); 00819 00821 void 00822 negate(); 00823 00825 void 00826 transpose(); 00827 00828 // 00829 // Assignment operators with scalar operand. 00830 // 00831 00833 SquareMatrix& 00834 operator+=(const value_type x); 00835 00837 SquareMatrix& 00838 operator-=(const value_type x); 00839 00841 SquareMatrix& 00842 operator*=(const value_type x); 00843 00845 SquareMatrix& 00846 operator/=(const value_type x); 00847 00849 SquareMatrix& 00850 operator%=(const value_type x); 00851 00852 // 00853 // Assignment operators with SquareMatrix operand 00854 // 00855 00857 template<typename T2> 00858 SquareMatrix& 00859 operator+=(const SquareMatrix<1,T2>& x); 00860 00862 template<typename T2> 00863 SquareMatrix& 00864 operator-=(const SquareMatrix<1,T2>& x); 00865 00867 template<typename T2> 00868 SquareMatrix& 00869 operator*=(const SquareMatrix<1,T2>& x); 00870 00871 // 00872 // Unary operators 00873 // 00874 00876 SquareMatrix 00877 operator+() { 00878 return *this; 00879 } 00880 00882 SquareMatrix 00883 operator-(); 00884 00885 }; 00886 00887 // 00888 // Binary operators 00889 // 00890 00892 template<typename T> 00893 inline 00894 SquareMatrix<1,T> 00895 operator+(const SquareMatrix<1,T>& m, const T x) { 00896 return SquareMatrix<1,T>(m[0] + x); 00897 } 00898 00900 template<typename T> 00901 inline 00902 SquareMatrix<1,T> 00903 operator+(const T x, const SquareMatrix<1,T>& m) { 00904 return m + x; 00905 } 00906 00908 template<typename T> 00909 inline 00910 SquareMatrix<1,T> 00911 operator+(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) { 00912 return SquareMatrix<1,T>(x[0] + y[0]); 00913 } 00914 00916 template<typename T> 00917 inline 00918 SquareMatrix<1,T> 00919 operator-(const SquareMatrix<1,T>& m, const T x) { 00920 return SquareMatrix<1,T>(m[0] - x); 00921 } 00922 00924 template<typename T> 00925 inline 00926 SquareMatrix<1,T> 00927 operator-(const T x, const SquareMatrix<1,T>& m) { 00928 return SquareMatrix<1,T>(x - m[0]); 00929 } 00930 00932 template<typename T> 00933 inline 00934 SquareMatrix<1,T> 00935 operator-(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) { 00936 return SquareMatrix<1,T>(x[0] - y[0]); 00937 } 00938 00940 template<typename T> 00941 inline 00942 SquareMatrix<1,T> 00943 operator*(const SquareMatrix<1,T>& m, const T x) { 00944 return SquareMatrix<1,T>(m[0] * x); 00945 } 00946 00948 template<typename T> 00949 inline 00950 SquareMatrix<1,T> 00951 operator*(const T x, const SquareMatrix<1,T>& m) { 00952 return m * x; 00953 } 00954 00956 template<typename T> 00957 inline 00958 SquareMatrix<1,T> 00959 operator*(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) { 00960 return SquareMatrix<1,T>(x[0]*y[0]); 00961 } 00962 00964 template<typename T> 00965 inline 00966 SquareMatrix<1,T> 00967 operator/(const SquareMatrix<1,T>& m, const T x) { 00968 #ifdef DEBUG_SquareMatrix 00969 assert(x != 0); 00970 #endif 00971 return SquareMatrix<1,T>(m[0] / x); 00972 } 00973 00975 template<typename T> 00976 inline 00977 SquareMatrix<1,T> 00978 operator/(const T x, const SquareMatrix<1,T>& m) { 00979 #ifdef DEBUG_SquareMatrix 00980 assert(m[0] != 0); 00981 #endif 00982 return SquareMatrix<1,T>(x / m[0]); 00983 } 00984 00985 // 00986 // Math operators. 00987 // 00988 00990 template<typename T> 00991 inline 00992 T 00993 computeSum(const SquareMatrix<1,T>& x) { 00994 return (x[0]); 00995 } 00996 00998 template<typename T> 00999 inline 01000 T 01001 computeProduct(const SquareMatrix<1,T>& x) { 01002 return (x[0]); 01003 } 01004 01006 template<typename T> 01007 inline 01008 T 01009 computeMinimum(const SquareMatrix<1,T>& x) { 01010 return x[0]; 01011 } 01012 01014 template<typename T> 01015 inline 01016 T 01017 computeMaximum(const SquareMatrix<1,T>& x) { 01018 return x[0]; 01019 } 01020 01022 template<typename T> 01023 inline 01024 T 01025 computeDeterminant(const SquareMatrix<1,T>& x) { 01026 return x[0]; 01027 } 01028 01030 template<typename T> 01031 inline 01032 T 01033 computeTrace(const SquareMatrix<1,T>& x) { 01034 return x[0]; 01035 } 01036 01038 template<typename T> 01039 inline 01040 SquareMatrix<1,T> 01041 computeTranspose(const SquareMatrix<1,T>& x) { 01042 return x; 01043 } 01044 01046 template<typename T> 01047 inline 01048 SquareMatrix<1,T> 01049 computeInverse(const SquareMatrix<1,T>& x) { 01050 SquareMatrix<1,T> y; 01051 computeInverse(x, &y); 01052 return y; 01053 } 01054 01056 template<typename T> 01057 inline 01058 void 01059 computeInverse(const SquareMatrix<1,T>& x, SquareMatrix<1,T>* y) { 01060 #ifdef DEBUG_SquareMatrix 01061 assert(x[0] != 0); 01062 #endif 01063 y[0] = 1 / x[0]; 01064 } 01065 01067 template<typename T> 01068 inline 01069 void 01070 computeScaledInverse(const SquareMatrix<1,T>& x, SquareMatrix<1,T>* si) { 01071 *si = x; 01072 } 01073 01075 template<typename T> 01076 inline 01077 SquareMatrix<1,T> 01078 computeScaledInverse(const SquareMatrix<1,T>& x) { 01079 return SquareMatrix<1,T>(1); 01080 } 01081 01083 template<typename T> 01084 inline 01085 T 01086 computeFrobeniusNorm(const SquareMatrix<1,T>& x) { 01087 return std::abs(x[0]); 01088 } 01089 01091 template<typename T> 01092 inline 01093 T 01094 computeFrobeniusNormSquared(const SquareMatrix<1,T>& x) { 01095 return (x[0]*x[0]); 01096 } 01097 01099 template<typename T> 01100 inline 01101 T 01102 computeInnerProduct(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) { 01103 return x * y; 01104 } 01105 01107 template<typename T> 01108 inline 01109 void 01110 computeProduct(const SquareMatrix<1,T>& m, const FixedArray<1,T>& v, 01111 FixedArray<1,T>* x) { 01112 (*x)[0] = m[0] * v[0]; 01113 } 01114 01115 // 01116 // Equality 01117 // 01118 01120 template<typename T1, typename T2> 01121 inline 01122 bool 01123 operator==(const SquareMatrix<1,T1>& a, const SquareMatrix<1,T2>& b) { 01124 return (a[0] == b[0]); 01125 } 01126 01128 template<typename T1, typename T2> 01129 inline 01130 bool 01131 operator!=(const SquareMatrix<1,T1>& a, const SquareMatrix<1,T2>& b) { 01132 return !(a == b); 01133 } 01134 01135 // 01136 // I/O 01137 // 01138 01140 template<typename T> 01141 std::ostream& 01142 operator<<(std::ostream& out, const SquareMatrix<1,T>& x); 01143 01145 template<typename T> 01146 std::istream& 01147 operator>>(std::istream& in, SquareMatrix<1,T>& x); 01148 01149 01150 01151 //---------------------------------------------------------------------------- 01152 // 2x2 matrices. 01153 //---------------------------------------------------------------------------- 01154 01156 01159 template<typename T> 01160 class SquareMatrix<2,T> : 01161 public TensorTypes<T> { 01162 private: 01163 01164 // 01165 // private types 01166 // 01167 01168 typedef TensorTypes<T> base_type; 01169 01170 public: 01171 01172 // 01173 // public types 01174 // 01175 01177 typedef typename base_type::value_type value_type; 01179 typedef typename base_type::pointer pointer; 01181 typedef typename base_type::const_pointer const_pointer; 01183 typedef typename base_type::iterator iterator; 01185 typedef typename base_type::const_iterator const_iterator; 01187 typedef typename base_type::reference reference; 01189 typedef typename base_type::const_reference const_reference; 01191 typedef typename base_type::size_type size_type; 01193 typedef typename base_type::difference_type difference_type; 01195 typedef typename base_type::index_type index_type; 01196 01197 protected: 01198 01199 // 01200 // Data 01201 // 01202 01204 value_type _elem[4]; 01205 01206 public: 01207 01208 // 01209 // Constructors 01210 // 01211 01213 SquareMatrix() 01214 {} 01215 01217 ~SquareMatrix() 01218 {} 01219 01221 SquareMatrix(const SquareMatrix& x); 01222 01224 template<typename T2> 01225 SquareMatrix(const SquareMatrix<2,T2>& x); 01226 01228 SquareMatrix(const value_type e00, const value_type e01, 01229 const value_type e10, const value_type e11); 01230 01232 SquareMatrix(const_pointer x); 01233 01235 SquareMatrix(const value_type x); 01236 01237 // 01238 // Assignment operators 01239 // 01240 01242 SquareMatrix& 01243 operator=(const SquareMatrix& x); 01244 01246 SquareMatrix& 01247 operator=(const value_type x); 01248 01250 template<typename T2> 01251 SquareMatrix& 01252 operator=(const SquareMatrix<2,T2>& x); 01253 01254 // 01255 // Accessors and Manipulators 01256 // 01257 01259 iterator 01260 begin() { 01261 return _elem; 01262 } 01263 01265 iterator 01266 end() { 01267 return _elem + 4; 01268 } 01269 01271 const_iterator 01272 begin() const { 01273 return _elem; 01274 } 01275 01277 const_iterator 01278 end() const { 01279 return _elem + 4; 01280 } 01281 01283 pointer 01284 data() { 01285 return _elem; 01286 } 01287 01289 const_pointer 01290 data() const { 01291 return _elem; 01292 } 01293 01295 size_type 01296 size() const { 01297 return 4; 01298 } 01299 01301 value_type 01302 operator()(const index_type i) const { 01303 #ifdef DEBUG_SquareMatrix 01304 assert(0 <= i && i < 4); 01305 #endif 01306 return _elem[i]; 01307 } 01308 01310 reference 01311 operator()(const index_type i) { 01312 #ifdef DEBUG_SquareMatrix 01313 assert(0 <= i && i < 4); 01314 #endif 01315 return _elem[i]; 01316 } 01317 01319 value_type 01320 operator[](const index_type i) const { 01321 #ifdef DEBUG_SquareMatrix 01322 assert(0 <= i && i < 4); 01323 #endif 01324 return _elem[i]; 01325 } 01326 01328 value_type& 01329 operator[](const index_type i) { 01330 #ifdef DEBUG_SquareMatrix 01331 assert(0 <= i && i < 4); 01332 #endif 01333 return _elem[i]; 01334 } 01335 01337 value_type 01338 operator()(const index_type i, const index_type j) const { 01339 #ifdef DEBUG_SquareMatrix 01340 assert(0 <= i && i < 2 && 0 <= j && j < 2); 01341 #endif 01342 return _elem[i * 2 + j]; 01343 } 01344 01346 reference 01347 operator()(const index_type i, const index_type j) { 01348 #ifdef DEBUG_SquareMatrix 01349 assert(0 <= i && i < 2 && 0 <= j && j < 2); 01350 #endif 01351 return _elem[i * 2 + j]; 01352 } 01353 01355 void 01356 get(pointer x) const; 01357 01359 void 01360 set(const_pointer x); 01361 01363 void 01364 get(reference e00, reference e01, 01365 reference e10, reference e11) const; 01366 01368 void 01369 set(const value_type e00, const value_type e01, 01370 const value_type e10, const value_type e11); 01371 01373 void 01374 negate(); 01375 01377 void 01378 transpose(); 01379 01380 // 01381 // Assignment operators with scalar operand. 01382 // 01383 01385 SquareMatrix& 01386 operator+=(const value_type x); 01387 01389 SquareMatrix& 01390 operator-=(const value_type x); 01391 01393 SquareMatrix& 01394 operator*=(const value_type x); 01395 01397 SquareMatrix& 01398 operator/=(const value_type x); 01399 01401 SquareMatrix& 01402 operator%=(const value_type x); 01403 01404 // 01405 // Assignment operators with SquareMatrix operand 01406 // 01407 01409 template<typename T2> 01410 SquareMatrix& 01411 operator+=(const SquareMatrix<2,T2>& x); 01412 01414 template<typename T2> 01415 SquareMatrix& 01416 operator-=(const SquareMatrix<2,T2>& x); 01417 01419 template<typename T2> 01420 SquareMatrix& 01421 operator*=(const SquareMatrix<2,T2>& x); 01422 01423 // 01424 // Unary operators 01425 // 01426 01428 SquareMatrix 01429 operator+() { 01430 return *this; 01431 } 01432 01434 SquareMatrix 01435 operator-(); 01436 01437 }; 01438 01439 // 01440 // Binary operators 01441 // 01442 01444 template<typename T> 01445 inline 01446 SquareMatrix<2,T> 01447 operator+(const SquareMatrix<2,T>& m, const T x) { 01448 return SquareMatrix<2,T>(m[0] + x, m[1] + x, 01449 m[2] + x, m[3] + x); 01450 } 01451 01453 template<typename T> 01454 inline 01455 SquareMatrix<2,T> 01456 operator+(const T x, const SquareMatrix<2,T>& m) { 01457 return m + x; 01458 } 01459 01461 template<typename T> 01462 inline 01463 SquareMatrix<2,T> 01464 operator+(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) { 01465 return SquareMatrix<2,T>(x[0] + y[0], x[1] + y[1], 01466 x[2] + y[2], x[3] + y[3]); 01467 } 01468 01470 template<typename T> 01471 inline 01472 SquareMatrix<2,T> 01473 operator-(const SquareMatrix<2,T>& m, const T x) { 01474 return SquareMatrix<2,T>(m[0] - x, m[1] - x, 01475 m[2] - x, m[3] - x); 01476 } 01477 01479 template<typename T> 01480 inline 01481 SquareMatrix<2,T> 01482 operator-(const T x, const SquareMatrix<2,T>& m) { 01483 return SquareMatrix<2,T>(x - m[0], x - m[1], 01484 x - m[2], x - m[3]); 01485 } 01486 01488 template<typename T> 01489 inline 01490 SquareMatrix<2,T> 01491 operator-(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) { 01492 return SquareMatrix<2,T>(x[0] - y[0], x[1] - y[1], 01493 x[2] - y[2], x[3] - y[3]); 01494 } 01495 01497 template<typename T> 01498 inline 01499 SquareMatrix<2,T> 01500 operator*(const SquareMatrix<2,T>& m, const T x) { 01501 return SquareMatrix<2,T>(m[0] * x, m[1] * x, 01502 m[2] * x, m[3] * x); 01503 } 01504 01506 template<typename T> 01507 inline 01508 SquareMatrix<2,T> 01509 operator*(const T x, const SquareMatrix<2,T>& m) { 01510 return m * x; 01511 } 01512 01514 template<typename T> 01515 inline 01516 SquareMatrix<2,T> 01517 operator*(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) { 01518 return SquareMatrix<2,T>(x[0]*y[0]+x[1]*y[2], 01519 x[0]*y[1]+x[1]*y[3], 01520 x[2]*y[0]+x[3]*y[2], 01521 x[2]*y[1]+x[3]*y[3]); 01522 } 01523 01525 template<typename T> 01526 inline 01527 SquareMatrix<2,T> 01528 operator/(const SquareMatrix<2,T>& m, const T x) { 01529 #ifdef DEBUG_SquareMatrix 01530 assert(x != 0); 01531 #endif 01532 return SquareMatrix<2,T>(m[0] / x, m[1] / x, 01533 m[2] / x, m[3] / x); 01534 } 01535 01537 template<typename T> 01538 inline 01539 SquareMatrix<2,T> 01540 operator/(const T x, const SquareMatrix<2,T>& m) { 01541 #ifdef DEBUG_SquareMatrix 01542 assert(m[0] != 0 && m[1] != 0 && 01543 m[2] != 0 && m[3] != 0); 01544 #endif 01545 return SquareMatrix<2,T>(x / m[0], x / m[1], 01546 x / m[2], x / m[3]); 01547 } 01548 01549 // 01550 // Math operators. 01551 // 01552 01554 template<typename T> 01555 inline 01556 T 01557 computeSum(const SquareMatrix<2,T>& x) { 01558 return (x[0] + x[1] + 01559 x[2] + x[3]); 01560 } 01561 01563 template<typename T> 01564 inline 01565 T 01566 computeProduct(const SquareMatrix<2,T>& x) { 01567 return (x[0] * x[1] * 01568 x[2] * x[3]); 01569 } 01570 01572 template<typename T> 01573 inline 01574 T 01575 computeMinimum(const SquareMatrix<2,T>& x) { 01576 return *std::min_element(x.begin(), x.end()); 01577 } 01578 01580 template<typename T> 01581 inline 01582 T 01583 computeMaximum(const SquareMatrix<2,T>& x) { 01584 return *std::max_element(x.begin(), x.end()); 01585 } 01586 01588 template<typename T> 01589 inline 01590 T 01591 computeDeterminant(const SquareMatrix<2,T>& x) { 01592 return (x[0] * x[3] - x[1] * x[2]); 01593 } 01594 01596 template<typename T> 01597 inline 01598 T 01599 computeTrace(const SquareMatrix<2,T>& x) { 01600 return (x[0] + x[3]); 01601 } 01602 01604 template<typename T> 01605 inline 01606 SquareMatrix<2,T> 01607 computeTranspose(const SquareMatrix<2,T>& x) { 01608 return SquareMatrix<2,T>(x[0], x[2], 01609 x[1], x[3]); 01610 } 01611 01613 template<typename T> 01614 inline 01615 SquareMatrix<2,T> 01616 computeInverse(const SquareMatrix<2,T>& x) { 01617 SquareMatrix<2,T> y; 01618 computeInverse(x, &y); 01619 return y; 01620 } 01621 01623 template<typename T> 01624 inline 01625 SquareMatrix<2,T> 01626 computeInverse(const SquareMatrix<2,T>& x, const T det) { 01627 SquareMatrix<2,T> y; 01628 computeInverse(x, det, &y); 01629 return y; 01630 } 01631 01633 template<typename T> 01634 inline 01635 void 01636 computeInverse(const SquareMatrix<2,T>& x, SquareMatrix<2,T>* y) { 01637 const T det = computeDeterminant(x); 01638 computeInverse(x, det, y); 01639 } 01640 01642 template<typename T> 01643 inline 01644 void 01645 computeInverse(const SquareMatrix<2,T>& x, const T det, SquareMatrix<2,T>* y) { 01646 #ifdef DEBUG_SquareMatrix 01647 assert(det != 0); 01648 #endif 01649 y->set(x[3] / det, - x[1] / det, 01650 - x[2] / det, x[0] / det); 01651 } 01652 01654 template<typename T> 01655 inline 01656 void 01657 computeScaledInverse(const SquareMatrix<2,T>& x, SquareMatrix<2,T>* si) { 01658 (*si)[0] = x[3]; (*si)[1] = -x[1]; 01659 (*si)[2] = -x[2]; (*si)[3] = x[0]; 01660 } 01661 01663 template<typename T> 01664 inline 01665 SquareMatrix<2,T> 01666 computeScaledInverse(const SquareMatrix<2,T>& x) { 01667 return SquareMatrix<2,T>(x[3], - x[1], 01668 - x[2], x[0]); 01669 } 01670 01672 template<typename T> 01673 inline 01674 T 01675 computeFrobeniusNorm(const SquareMatrix<2,T>& x) { 01676 return std::sqrt(computeFrobeniusNormSquared(x)); 01677 } 01678 01680 template<typename T> 01681 inline 01682 T 01683 computeFrobeniusNormSquared(const SquareMatrix<2,T>& x) { 01684 return (x[0]*x[0] + x[1]*x[1] + 01685 x[2]*x[2] + x[3]*x[3]); 01686 } 01687 01689 template<typename T> 01690 inline 01691 T 01692 computeInnerProduct(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) { 01693 return x[0] * y[0] + x[1] * y[1] + x[2] * y[2] + x[3] * y[3]; 01694 } 01695 01697 template<typename T> 01698 inline 01699 void 01700 computeProduct(const SquareMatrix<2,T>& m, const FixedArray<2,T>& v, 01701 FixedArray<2,T>* x) { 01702 (*x)[0] = m[0] * v[0] + m[1] * v[1]; 01703 (*x)[1] = m[2] * v[0] + m[3] * v[1]; 01704 } 01705 01706 // 01707 // Equality 01708 // 01709 01711 template<typename T1, typename T2> 01712 inline 01713 bool 01714 operator==(const SquareMatrix<2,T1>& a, const SquareMatrix<2,T2>& b) { 01715 return (a[0] == b[0] && a[1] == b[1] && 01716 a[2] == b[2] && a[3] == b[3]); 01717 } 01718 01720 template<typename T1, typename T2> 01721 inline 01722 bool 01723 operator!=(const SquareMatrix<2,T1>& a, const SquareMatrix<2,T2>& b) { 01724 return !(a == b); 01725 } 01726 01727 // 01728 // I/O 01729 // 01730 01732 template<typename T> 01733 std::ostream& 01734 operator<<(std::ostream& out, const SquareMatrix<2,T>& x); 01735 01737 template<typename T> 01738 std::istream& 01739 operator>>(std::istream& in, SquareMatrix<2,T>& x); 01740 01741 01742 01743 //---------------------------------------------------------------------------- 01744 // 3x3 matrices. 01745 //---------------------------------------------------------------------------- 01746 01748 01751 template<typename T> 01752 class SquareMatrix<3,T> : 01753 public TensorTypes<T> { 01754 private: 01755 01756 // 01757 // private types 01758 // 01759 01760 typedef TensorTypes<T> base_type; 01761 01762 public: 01763 01764 // 01765 // public types 01766 // 01767 01769 typedef typename base_type::value_type value_type; 01771 typedef typename base_type::pointer pointer; 01773 typedef typename base_type::const_pointer const_pointer; 01775 typedef typename base_type::iterator iterator; 01777 typedef typename base_type::const_iterator const_iterator; 01779 typedef typename base_type::reference reference; 01781 typedef typename base_type::const_reference const_reference; 01783 typedef typename base_type::size_type size_type; 01785 typedef typename base_type::difference_type difference_type; 01787 typedef typename base_type::index_type index_type; 01788 01789 protected: 01790 01791 // 01792 // Data 01793 // 01794 01796 value_type _elem[9]; 01797 01798 public: 01799 01800 // 01801 // Constructors 01802 // 01803 01805 SquareMatrix() 01806 {} 01807 01809 ~SquareMatrix() 01810 {} 01811 01813 SquareMatrix(const SquareMatrix& x); 01814 01816 template<typename T2> 01817 SquareMatrix(const SquareMatrix<3,T2>& x); 01818 01820 SquareMatrix(const value_type e00, const value_type e01, const value_type e02, 01821 const value_type e10, const value_type e11, const value_type e12, 01822 const value_type e20, const value_type e21, const value_type e22); 01823 01825 SquareMatrix(const_pointer x); 01826 01828 SquareMatrix(const value_type x); 01829 01830 // 01831 // Assignment operators 01832 // 01833 01835 SquareMatrix& 01836 operator=(const SquareMatrix& x); 01837 01839 SquareMatrix& 01840 operator=(const value_type x); 01841 01843 template<typename T2> 01844 SquareMatrix& 01845 operator=(const SquareMatrix<3,T2>& x); 01846 01847 // 01848 // Accessors and Manipulators 01849 // 01850 01852 iterator 01853 begin() { 01854 return _elem; 01855 } 01856 01858 iterator 01859 end() { 01860 return _elem + 9; 01861 } 01862 01864 const_iterator 01865 begin() const { 01866 return _elem; 01867 } 01868 01870 const_iterator 01871 end() const { 01872 return _elem + 9; 01873 } 01874 01876 pointer 01877 data() { 01878 return _elem; 01879 } 01880 01882 const_pointer 01883 data() const { 01884 return _elem; 01885 } 01886 01888 size_type 01889 size() const { 01890 return 9; 01891 } 01892 01894 value_type 01895 operator()(const index_type i) const { 01896 #ifdef DEBUG_SquareMatrix 01897 assert(0 <= i && i < 9); 01898 #endif 01899 return _elem[i]; 01900 } 01901 01903 reference 01904 operator()(const index_type i) { 01905 #ifdef DEBUG_SquareMatrix 01906 assert(0 <= i && i < 9); 01907 #endif 01908 return _elem[i]; 01909 } 01910 01912 value_type 01913 operator[](const index_type i) const { 01914 #ifdef DEBUG_SquareMatrix 01915 assert(0 <= i && i < 9); 01916 #endif 01917 return _elem[i]; 01918 } 01919 01921 value_type& 01922 operator[](const index_type i) { 01923 #ifdef DEBUG_SquareMatrix 01924 assert(0 <= i && i < 9); 01925 #endif 01926 return _elem[i]; 01927 } 01928 01930 value_type 01931 operator()(const index_type i, const index_type j) const { 01932 #ifdef DEBUG_SquareMatrix 01933 assert(0 <= i && i < 3 && 0 <= j && j < 3); 01934 #endif 01935 return _elem[i * 3 + j]; 01936 } 01937 01939 reference 01940 operator()(const index_type i, const index_type j) { 01941 #ifdef DEBUG_SquareMatrix 01942 assert(0 <= i && i < 3 && 0 <= j && j < 3); 01943 #endif 01944 return _elem[i * 3 + j]; 01945 } 01946 01948 void 01949 get(pointer x) const; 01950 01952 void 01953 set(const_pointer x); 01954 01956 void 01957 get(reference e00, reference e01, reference e02, 01958 reference e10, reference e11, reference e12, 01959 reference e20, reference e21, reference e22) const; 01960 01962 void 01963 set(const value_type e00, const value_type e01, const value_type e02, 01964 const value_type e10, const value_type e11, const value_type e12, 01965 const value_type e20, const value_type e21, const value_type e22); 01966 01968 void 01969 negate(); 01970 01972 void 01973 transpose(); 01974 01975 // 01976 // Assignment operators with scalar operand. 01977 // 01978 01980 SquareMatrix& 01981 operator+=(const value_type x); 01982 01984 SquareMatrix& 01985 operator-=(const value_type x); 01986 01988 SquareMatrix& 01989 operator*=(const value_type x); 01990 01992 SquareMatrix& 01993 operator/=(const value_type x); 01994 01996 SquareMatrix& 01997 operator%=(const value_type x); 01998 01999 // 02000 // Assignment operators with SquareMatrix operand 02001 // 02002 02004 template<typename T2> 02005 SquareMatrix& 02006 operator+=(const SquareMatrix<3,T2>& x); 02007 02009 template<typename T2> 02010 SquareMatrix& 02011 operator-=(const SquareMatrix<3,T2>& x); 02012 02014 template<typename T2> 02015 SquareMatrix& 02016 operator*=(const SquareMatrix<3,T2>& x); 02017 02018 // 02019 // Unary operators 02020 // 02021 02023 SquareMatrix 02024 operator+() { 02025 return *this; 02026 } 02027 02029 SquareMatrix 02030 operator-(); 02031 02032 }; 02033 02034 // 02035 // Binary operators 02036 // 02037 02039 template<typename T> 02040 inline 02041 SquareMatrix<3,T> 02042 operator+(const SquareMatrix<3,T>& m, const T x) { 02043 return SquareMatrix<3,T>(m[0] + x, m[1] + x, m[2] + x, 02044 m[3] + x, m[4] + x, m[5] + x, 02045 m[6] + x, m[7] + x, m[8] + x); 02046 } 02047 02049 template<typename T> 02050 inline 02051 SquareMatrix<3,T> 02052 operator+(const T x, const SquareMatrix<3,T>& m) { 02053 return m + x; 02054 } 02055 02057 template<typename T> 02058 inline 02059 SquareMatrix<3,T> 02060 operator+(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) { 02061 return SquareMatrix<3,T>(x[0] + y[0], x[1] + y[1], x[2] + y[2], 02062 x[3] + y[3], x[4] + y[4], x[5] + y[5], 02063 x[6] + y[6], x[7] + y[7], x[8] + y[8]); 02064 } 02065 02067 template<typename T> 02068 inline 02069 SquareMatrix<3,T> 02070 operator-(const SquareMatrix<3,T>& m, const T x) { 02071 return SquareMatrix<3,T>(m[0] - x, m[1] - x, m[2] - x, 02072 m[3] - x, m[4] - x, m[5] - x, 02073 m[6] - x, m[7] - x, m[8] - x); 02074 } 02075 02077 template<typename T> 02078 inline 02079 SquareMatrix<3,T> 02080 operator-(const T x, const SquareMatrix<3,T>& m) { 02081 return SquareMatrix<3,T>(x - m[0], x - m[1], x - m[2], 02082 x - m[3], x - m[4], x - m[5], 02083 x - m[6], x - m[7], x - m[8]); 02084 } 02085 02087 template<typename T> 02088 inline 02089 SquareMatrix<3,T> 02090 operator-(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) { 02091 return SquareMatrix<3,T>(x[0] - y[0], x[1] - y[1], x[2] - y[2], 02092 x[3] - y[3], x[4] - y[4], x[5] - y[5], 02093 x[6] - y[6], x[7] - y[7], x[8] - y[8]); 02094 } 02095 02097 template<typename T> 02098 inline 02099 SquareMatrix<3,T> 02100 operator*(const SquareMatrix<3,T>& m, const T x) { 02101 return SquareMatrix<3,T>(m[0] * x, m[1] * x, m[2] * x, 02102 m[3] * x, m[4] * x, m[5] * x, 02103 m[6] * x, m[7] * x, m[8] * x); 02104 } 02105 02107 template<typename T> 02108 inline 02109 SquareMatrix<3,T> 02110 operator*(const T x, const SquareMatrix<3,T>& m) { 02111 return m * x; 02112 } 02113 02115 template<typename T> 02116 inline 02117 SquareMatrix<3,T> 02118 operator*(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) { 02119 return SquareMatrix<3,T> 02120 (x[0]*y[0]+x[1]*y[3]+x[2]*y[6], 02121 x[0]*y[1]+x[1]*y[4]+x[2]*y[7], 02122 x[0]*y[2]+x[1]*y[5]+x[2]*y[8], 02123 x[3]*y[0]+x[4]*y[3]+x[5]*y[6], 02124 x[3]*y[1]+x[4]*y[4]+x[5]*y[7], 02125 x[3]*y[2]+x[4]*y[5]+x[5]*y[8], 02126 x[6]*y[0]+x[7]*y[3]+x[8]*y[6], 02127 x[6]*y[1]+x[7]*y[4]+x[8]*y[7], 02128 x[6]*y[2]+x[7]*y[5]+x[8]*y[8]); 02129 } 02130 02132 template<typename T> 02133 inline 02134 SquareMatrix<3,T> 02135 operator/(const SquareMatrix<3,T>& m, const T x) { 02136 #ifdef DEBUG_SquareMatrix 02137 assert(x != 0); 02138 #endif 02139 return SquareMatrix<3,T>(m[0] / x, m[1] / x, m[2] / x, 02140 m[3] / x, m[4] / x, m[5] / x, 02141 m[6] / x, m[7] / x, m[8] / x); 02142 } 02143 02145 template<typename T> 02146 inline 02147 SquareMatrix<3,T> 02148 operator/(const T x, const SquareMatrix<3,T>& m) { 02149 #ifdef DEBUG_SquareMatrix 02150 assert(m[0] != 0 && m[1] != 0 && m[2] != 0 && 02151 m[3] != 0 && m[4] != 0 && m[5] != 0 && 02152 m[6] != 0 && m[7] != 0 && m[8] != 0); 02153 #endif 02154 return SquareMatrix<3,T>(x / m[0], x / m[1], x / m[2], 02155 x / m[3], x / m[4], x / m[5], 02156 x / m[6], x / m[7], x / m[8]); 02157 } 02158 02159 // 02160 // Math operators. 02161 // 02162 02164 template<typename T> 02165 inline 02166 T 02167 computeSum(const SquareMatrix<3,T>& x) { 02168 return (x[0] + x[1] + x[2] + 02169 x[3] + x[4] + x[5] + 02170 x[6] + x[7] + x[8]); 02171 } 02172 02174 template<typename T> 02175 inline 02176 T 02177 computeProduct(const SquareMatrix<3,T>& x) { 02178 return (x[0] * x[1] * x[2] * 02179 x[3] * x[4] * x[5] * 02180 x[6] * x[7] * x[8]); 02181 } 02182 02184 template<typename T> 02185 inline 02186 T 02187 computeMinimum(const SquareMatrix<3,T>& x) { 02188 return *std::min_element(x.begin(), x.end()); 02189 } 02190 02192 template<typename T> 02193 inline 02194 T 02195 computeMaximum(const SquareMatrix<3,T>& x) { 02196 return *std::max_element(x.begin(), x.end()); 02197 } 02198 02200 template<typename T> 02201 inline 02202 T 02203 computeDeterminant(const SquareMatrix<3,T>& x) { 02204 return (x[0]*x[4]*x[8] + x[1]*x[5]*x[6] + x[2]*x[3]*x[7] - 02205 x[2]*x[4]*x[6] - x[1]*x[3]*x[8] - x[0]*x[5]*x[7]); 02206 } 02207 02209 template<typename T> 02210 inline 02211 T 02212 computeTrace(const SquareMatrix<3,T>& x) { 02213 return (x[0] + x[4] + x[8]); 02214 } 02215 02217 template<typename T> 02218 inline 02219 SquareMatrix<3,T> 02220 computeTranspose(const SquareMatrix<3,T>& x) { 02221 return SquareMatrix<3,T>(x[0], x[3], x[6], 02222 x[1], x[4], x[7], 02223 x[2], x[5], x[8]); 02224 } 02225 02227 template<typename T> 02228 inline 02229 SquareMatrix<3,T> 02230 computeInverse(const SquareMatrix<3,T>& x) { 02231 SquareMatrix<3,T> y; 02232 computeInverse(x, &y); 02233 return y; 02234 } 02235 02237 template<typename T> 02238 inline 02239 SquareMatrix<3,T> 02240 computeInverse(const SquareMatrix<3,T>& x, const T det) { 02241 SquareMatrix<3,T> y; 02242 computeInverse(x, det, &y); 02243 return y; 02244 } 02245 02247 template<typename T> 02248 inline 02249 void 02250 computeInverse(const SquareMatrix<3,T>& x, SquareMatrix<3,T>* y) { 02251 const T det = computeDeterminant(x); 02252 computeInverse(x, det, y); 02253 } 02254 02256 template<typename T> 02257 inline 02258 void 02259 computeInverse(const SquareMatrix<3,T>& x, const T det, SquareMatrix<3,T>* y) { 02260 #ifdef DEBUG_SquareMatrix 02261 assert(det != 0); 02262 #endif 02263 y->set((x[4] * x[8] - x[5] * x[7]) / det, 02264 (x[2] * x[7] - x[1] * x[8]) / det, 02265 (x[1] * x[5] - x[2] * x[4]) / det, 02266 (x[5] * x[6] - x[3] * x[8]) / det, 02267 (x[0] * x[8] - x[2] * x[6]) / det, 02268 (x[2] * x[3] - x[0] * x[5]) / det, 02269 (x[3] * x[7] - x[4] * x[6]) / det, 02270 (x[1] * x[6] - x[0] * x[7]) / det, 02271 (x[0] * x[4] - x[1] * x[3]) / det); 02272 } 02273 02275 template<typename T> 02276 inline 02277 void 02278 computeScaledInverse(const SquareMatrix<3,T>& x, SquareMatrix<3,T>* si) { 02279 (*si)[0] = x[4] * x[8] - x[5] * x[7]; 02280 (*si)[1] = x[2] * x[7] - x[1] * x[8]; 02281 (*si)[2] = x[1] * x[5] - x[2] * x[4]; 02282 (*si)[3] = x[5] * x[6] - x[3] * x[8]; 02283 (*si)[4] = x[0] * x[8] - x[2] * x[6]; 02284 (*si)[5] = x[2] * x[3] - x[0] * x[5]; 02285 (*si)[6] = x[3] * x[7] - x[4] * x[6]; 02286 (*si)[7] = x[1] * x[6] - x[0] * x[7]; 02287 (*si)[8] = x[0] * x[4] - x[1] * x[3]; 02288 } 02289 02291 template<typename T> 02292 inline 02293 SquareMatrix<3,T> 02294 computeScaledInverse(const SquareMatrix<3,T>& x) { 02295 return SquareMatrix<3,T>(x[4] * x[8] - x[5] * x[7], 02296 x[2] * x[7] - x[1] * x[8], 02297 x[1] * x[5] - x[2] * x[4], 02298 x[5] * x[6] - x[3] * x[8], 02299 x[0] * x[8] - x[2] * x[6], 02300 x[2] * x[3] - x[0] * x[5], 02301 x[3] * x[7] - x[4] * x[6], 02302 x[1] * x[6] - x[0] * x[7], 02303 x[0] * x[4] - x[1] * x[3]); 02304 } 02305 02307 template<typename T> 02308 inline 02309 T 02310 computeFrobeniusNorm(const SquareMatrix<3,T>& x) { 02311 return std::sqrt(computeFrobeniusNormSquared(x)); 02312 } 02313 02315 template<typename T> 02316 inline 02317 T 02318 computeFrobeniusNormSquared(const SquareMatrix<3,T>& x) { 02319 return (x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + 02320 x[3]*x[3] + x[4]*x[4] + x[5]*x[5] + 02321 x[6]*x[6] + x[7]*x[7] + x[8]*x[8]); 02322 } 02323 02325 template<typename T> 02326 inline 02327 T 02328 computeInnerProduct(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) { 02329 return x[0] * y[0] + x[1] * y[1] + x[2] * y[2] + x[3] * y[3] + x[4] * y[4] 02330 + x[5] * y[5] + x[6] * y[6] + x[7] * y[7] + x[8] * y[8]; 02331 } 02332 02334 template<typename T> 02335 inline 02336 void 02337 computeProduct(const SquareMatrix<3,T>& m, const FixedArray<3,T>& v, 02338 FixedArray<3,T>* x) { 02339 (*x)[0] = m[0] * v[0] + m[1] * v[1] + m[2] * v[2]; 02340 (*x)[1] = m[3] * v[0] + m[4] * v[1] + m[5] * v[2]; 02341 (*x)[2] = m[6] * v[0] + m[7] * v[1] + m[8] * v[2]; 02342 } 02343 02344 // 02345 // Equality 02346 // 02347 02349 template<typename T1, typename T2> 02350 inline 02351 bool 02352 operator==(const SquareMatrix<3,T1>& a, const SquareMatrix<3,T2>& b) { 02353 return (a[0] == b[0] && a[1] == b[1] && a[2] == b[2] && 02354 a[3] == b[3] && a[4] == b[4] && a[5] == b[5] && 02355 a[6] == b[6] && a[7] == b[7] && a[8] == b[8]); 02356 } 02357 02359 template<typename T1, typename T2> 02360 inline 02361 bool 02362 operator!=(const SquareMatrix<3,T1>& a, const SquareMatrix<3,T2>& b) { 02363 return !(a == b); 02364 } 02365 02366 // 02367 // I/O 02368 // 02369 02371 template<typename T> 02372 std::ostream& 02373 operator<<(std::ostream& out, const SquareMatrix<3,T>& x); 02374 02376 template<typename T> 02377 std::istream& 02378 operator>>(std::istream& in, SquareMatrix<3,T>& x); 02379 02380 END_NAMESPACE_ADS 02381 02382 #define __SquareMatrix_ipp__ 02383 #include "SquareMatrix.ipp" 02384 #undef __SquareMatrix_ipp__ 02385 02386 #endif
Generated on Thu Jun 30 02:14:50 2016 for Algorithms and Data Structures Package by
1.6.3