MeshKit  1.0
Matrix.hpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 
00033 #ifndef MESHKIT_MATRIX_HPP
00034 #define MESHKIT_MATRIX_HPP
00035 
00036 #include <string>
00037 #include <istream>
00038 #include <ostream>
00039 #include <sstream>
00040 #include <cmath>
00041 
00042 namespace MeshKit {
00043 
00049 template <unsigned R, unsigned C>
00050 class Matrix {
00051 protected:
00052   double m[R*C];
00053 
00054 public:
00055   typedef Matrix<R,C> my_type;
00056 
00057   enum { ROWS = R, COLS = C };
00058 
00060   Matrix()                                         { }
00062   Matrix( double v )                               { diag(v); }
00064   Matrix( const double* v )                        { set(v); }
00066   Matrix( const double v[R][C] )                   { set(v); }
00068   Matrix( const Matrix<R,1>* c )                   { set_columns(c); }
00070   Matrix( const Matrix<1,C>* r )                   { set_rows(r); }
00072   Matrix( const char* s )                          { set(s); }
00074   Matrix( const std::string& s )                   { set(s); }
00079   Matrix( const Matrix<R+1,C+1>& m, unsigned r, unsigned c )
00080                                                    { make_minor(m,r,c); }
00081 
00082   Matrix<R,C>& operator=( double v )               { set(v); return *this; }
00083   Matrix<R,C>& operator=( const double* v )        { set(v); return *this; }
00084   Matrix<R,C>& operator=( const char* s )          { set(s); return *this; }
00085   Matrix<R,C>& operator=( const std::string& s )   { set(s); return *this; }
00086 
00087   double& operator()( unsigned r, unsigned c )        { return m[r*C+c]; }
00088   double  operator()( unsigned r, unsigned c ) const  { return m[r*C+c]; }
00089   double* data()                                      { return m; }
00090   const double* data() const                          { return m; }
00091 
00092   void zero()                                         { set(0.0);  }
00093   void identity()                                     { diag(1.0); }
00094   void set( double v ) {
00095     for (unsigned i = 0; i < R*C; ++i) m[i] = v;
00096   }
00097   void set( const double* v ) {
00098     for (unsigned i = 0; i < R*C; ++i) m[i] = v[i];
00099   }
00100   void set( const double v[R][C] );
00101   void set( const char* s )   { std::istringstream i(s); i >> *this; }
00102   void set( const std::string& s ) { set( s.c_str() ); }
00104   inline void diag( double v );
00106   inline void diag( const double* v );
00108   inline void make_minor( const Matrix<R+1,C+1>& m, unsigned r, unsigned c );
00109 
00111   inline Matrix<R,C>& assign_add_transpose( const Matrix<C,R>& other );
00113   inline Matrix<R,C>& assign_product( double s, const Matrix<R,C>& m );
00115   inline Matrix<R,C>& assign_add_product( double s, const Matrix<R,C>& m );
00117   inline Matrix<R,C>& assign_multiply_elements( const Matrix<R,C>& m );
00118 
00119   inline Matrix<R,C>& operator+=( const Matrix<R,C>& other );
00120   inline Matrix<R,C>& operator-=( const Matrix<R,C>& other );
00121   inline Matrix<R,C>& operator+=( double scalar );
00122   inline Matrix<R,C>& operator-=( double scalar );
00123   inline Matrix<R,C>& operator*=( double scalar );
00124   inline Matrix<R,C>& operator/=( double scalar );
00125 
00126   //      Matrix<1,C>& row( unsigned r )       { return *(Matrix<1,C>*)(m+C*r); }
00127   //const Matrix<1,C>& row( unsigned r ) const { return *(Matrix<1,C>*)(m+C*r); }
00128   Matrix<1,C> row( unsigned r ) const { return Matrix<1,C>(m+C*r); }
00129   Matrix<R,1> column( unsigned c ) const;
00130   Matrix<1,R> column_transpose( unsigned c ) const;
00131 
00132   void set_row( unsigned r, const Matrix<1,C>& v );
00133   void add_row( unsigned r, const Matrix<1,C>& v );
00134   void set_row_transpose( unsigned r, const Matrix<C,1>& v );
00135   void set_rows( const Matrix<1,C>* v );
00136   void set_column( unsigned c, const Matrix<R,1>& v );
00137   void add_column( unsigned c, const Matrix<R,1>& v );
00138   void set_column_transpose( unsigned c, const Matrix<1,R>& v );
00139   void set_columns( const Matrix<R,1>* v );
00140 };
00141 
00142 template <>
00143 class Matrix<1,1> {
00144 protected:
00145   double m;
00146 
00147 public:
00148   typedef Matrix<1,1> my_type;
00149 
00150   enum { ROWS = 1, COLS = 1 };
00151 
00153   Matrix()                                         { }
00155   Matrix( double v )                               : m(v) {}
00157   Matrix( const double* v )                        : m(*v) {}
00159   Matrix( const char* s )                          { set(s); }
00161   Matrix( const std::string& s )                   { set(s); }
00166   Matrix( const Matrix<2,2>& M, unsigned r, unsigned c ) :  m(M(r,c)) {}
00167 
00168   Matrix<1,1>& operator=( double v )                 { m = v; return *this; }
00169   Matrix<1,1>& operator=( const double* v )          { m = *v; return *this; }
00170   Matrix<1,1>& operator=( const char* s )            { set(s); return *this; }
00171   Matrix<1,1>& operator=( const std::string& s )     { set(s); return *this; }
00172 
00173   double& operator()( unsigned, unsigned )           { return m; }
00174   double  operator()( unsigned, unsigned ) const     { return m; }
00175   double* data()                                     { return &m; }
00176   const double* data() const                         { return &m; }
00177 
00178   void zero()                                        { m = 0.0; }
00179   void identity()                                    { m = 1.0; }
00180   void set( double v )        { m = v; }
00181   void set( const double* v ) { m= *v; }
00182   void set( const char* s )   { std::istringstream i(s); i >> m; }
00183   void set( const std::string& s ) { set( s.c_str() ); }
00185   inline void diag( double v ) { m = v; }
00187   inline void diag( const double* v ) { m = *v; }
00189   inline void make_minor( const Matrix<2,2>& M, unsigned r, unsigned c )
00190     { m = M(r,c); }
00191 
00193   inline Matrix<1,1>& assign_add_transpose( const Matrix<1,1>& other ) {
00194     m += other.m; return *this;
00195   }
00197   inline Matrix<1,1>& assign_product( double s, const Matrix<1,1>& other ) {
00198     m = s*other.m; return *this;
00199   }
00201   inline Matrix<1,1>& assign_add_product( double s, const Matrix<1,1>& other ) {
00202     m += s*other.m; return *this;
00203   }
00205   inline Matrix<1,1>& assign_multiply_elements( const Matrix<1,1>& other ) {
00206     m *= other.m; return *this;
00207   }
00208 
00209   operator double () const {
00210     return m;
00211   }
00212 };
00213 
00219 template <unsigned L>
00220 class Vector : public Matrix<L,1>
00221 {
00222 public:
00223   Vector()                                         { }
00224   Vector( double v )                               { Matrix<L,1>::set(v); }
00225   Vector( const double* v )                        { Matrix<L,1>::set(v); }
00226   Vector( const char* s )                          { Matrix<L,1>::set(s); }
00227   Vector( const std::string& s )                   { Matrix<L,1>::set(s); }
00228   Vector( const Matrix<L,1>& m)                    : Matrix<L,1>(m) {}
00229 
00230   double& operator[](unsigned idx)       { return Matrix<L,1>::operator()(idx,0); }
00231   double  operator[](unsigned idx) const { return Matrix<L,1>::operator()(idx,0); }
00232   double& operator()(unsigned idx)       { return Matrix<L,1>::operator()(idx,0); }
00233   double  operator()(unsigned idx) const { return Matrix<L,1>::operator()(idx,0); }
00234 
00235   Vector<L>& operator=( const Matrix<L,1>& m ) {
00236     Matrix<L,1>::operator=(m); return *this;
00237   }
00238 };
00239 
00240 template <unsigned R, unsigned C> inline
00241 void Matrix<R,C>::set( const double v[R][C] ) {
00242   for (unsigned r = 0; r < R; ++r)
00243     for (unsigned c = 0; c < C; ++c)
00244       operator()(r,c) = v[r][c];
00245 }
00246 
00247 template <unsigned R, unsigned C> inline
00248 void Matrix<R,C>::diag( double v ) {
00249   //for (unsigned r = 0; r < R; ++r)
00250   //  for (unsigned c = 0; c < C; ++c)
00251   //    operator()(r,c) = (r == c) ? v : 0.0;
00252 
00253   switch (R) {
00254   default: for (unsigned r = 4; r < R; ++r)
00255       switch (C) {
00256       default: for (unsigned k = 4; k < C; ++k)
00257                  operator()(r,k) = r == k ? v : 0.0;
00258       case 4:    operator()(r,3) = 0.0;
00259       case 3:    operator()(r,2) = 0.0;
00260       case 2:    operator()(r,1) = 0.0;
00261       case 1:    operator()(r,0) = 0.0;
00262       }
00263   case 4:
00264     switch (C) {
00265     default: for (unsigned k = 4; k < C; ++k)
00266                operator()(3,k) = 0.0;
00267     case 4:    operator()(3,3) = v;
00268     case 3:    operator()(3,2) = 0.0;
00269     case 2:    operator()(3,1) = 0.0;
00270     case 1:    operator()(3,0) = 0.0;
00271     }
00272   case 3:
00273     switch (C) {
00274     default: for (unsigned k = 4; k < C; ++k)
00275                operator()(2,k) = 0.0;
00276     case 4:    operator()(2,3) = 0.0;
00277     case 3:    operator()(2,2) = v;
00278     case 2:    operator()(2,1) = 0.0;
00279     case 1:    operator()(2,0) = 0.0;
00280     }
00281   case 2:
00282     switch (C) {
00283     default: for (unsigned k = 4; k < C; ++k)
00284                operator()(1,k) = 0.0;
00285     case 4:    operator()(1,3) = 0.0;
00286     case 3:    operator()(1,2) = 0.0;
00287     case 2:    operator()(1,1) = v;
00288     case 1:    operator()(1,0) = 0.0;
00289     }
00290   case 1:
00291     switch (C) {
00292     default: for (unsigned k = 4; k < C; ++k)
00293                operator()(0,k) = 0.0;
00294     case 4:    operator()(0,3) = 0.0;
00295     case 3:    operator()(0,2) = 0.0;
00296     case 2:    operator()(0,1) = 0.0;
00297     case 1:    operator()(0,0) = v;
00298     }
00299   }
00300 }
00301 
00302 template <unsigned R, unsigned C> inline
00303 void Matrix<R,C>::diag( const double* v ) {
00304   //for (unsigned r = 0; r < R; ++r)
00305   //  for (unsigned c = 0; c < C; ++c)
00306   //    operator()(r,c) = (r == c) ? v[r] : 0.0;
00307 
00308   switch (R) {
00309   default: for (unsigned r = 4; r < R; ++r)
00310       switch (C) {
00311       default: for (unsigned k = 4; k < C; ++k)
00312                  operator()(r,k) = r == k ? v[r] : 0.0;
00313       case 4:    operator()(r,3) = 0.0;
00314       case 3:    operator()(r,2) = 0.0;
00315       case 2:    operator()(r,1) = 0.0;
00316       case 1:    operator()(r,0) = 0.0;
00317       }
00318   case 4:
00319     switch (C) {
00320     default: for (unsigned k = 4; k < C; ++k)
00321                operator()(3,k) = 0.0;
00322     case 4:    operator()(3,3) = v[3];
00323     case 3:    operator()(3,2) = 0.0;
00324     case 2:    operator()(3,1) = 0.0;
00325     case 1:    operator()(3,0) = 0.0;
00326     }
00327   case 3:
00328     switch (C) {
00329     default: for (unsigned k = 4; k < C; ++k)
00330                operator()(2,k) = 0.0;
00331     case 4:    operator()(2,3) = 0.0;
00332     case 3:    operator()(2,2) = v[2];
00333     case 2:    operator()(2,1) = 0.0;
00334     case 1:    operator()(2,0) = 0.0;
00335     }
00336   case 2:
00337     switch (C) {
00338     default: for (unsigned k = 4; k < C; ++k)
00339                operator()(1,k) = 0.0;
00340     case 4:    operator()(1,3) = 0.0;
00341     case 3:    operator()(1,2) = 0.0;
00342     case 2:    operator()(1,1) = v[1];
00343     case 1:    operator()(1,0) = 0.0;
00344     }
00345   case 1:
00346     switch (C) {
00347     default: for (unsigned k = 4; k < C; ++k)
00348                operator()(0,k) = 0.0;
00349     case 4:    operator()(0,3) = 0.0;
00350     case 3:    operator()(0,2) = 0.0;
00351     case 2:    operator()(0,1) = 0.0;
00352     case 1:    operator()(0,0) = v[0];
00353     }
00354   }
00355 }
00356 
00362 template <unsigned R, unsigned C> inline
00363 void Matrix<R,C>::make_minor( const Matrix<R+1,C+1>& M, unsigned r, unsigned c ) {
00364   for (unsigned i = 0; i < r; ++i) {
00365     for (unsigned j = 0; j < c; ++j)
00366       operator()(i,j) = M(i,j);
00367     for (unsigned j = c; j < C; ++j)
00368       operator()(i,j) = M(i,j+1);
00369   }
00370   for (unsigned i = r; i < R; ++i) {
00371     for (unsigned j = 0; j < c; ++j)
00372       operator()(i,j) = M(i+1,j);
00373     for (unsigned j = c; j < C; ++j)
00374       operator()(i,j) = M(i+1,j+1);
00375   }
00376 }
00377 
00378 template <unsigned R, unsigned C> inline
00379 void Matrix<R,C>::set_row( unsigned r, const Matrix<1,C>& v ) {
00380   for (unsigned i = 0; i < C; ++i) operator()(r,i) = v(0,i);
00381 }
00382 template <unsigned R, unsigned C> inline
00383 void Matrix<R,C>::add_row( unsigned r, const Matrix<1,C>& v ) {
00384   for (unsigned i = 0; i < C; ++i) operator()(r,i) += v(0,i);
00385 }
00386 template <unsigned R, unsigned C> inline
00387 void Matrix<R,C>::set_row_transpose( unsigned r, const Matrix<C,1>& v ) {
00388   for (unsigned i = 0; i < C; ++i) operator()(r,i) = v(i,0);
00389 }
00390 template <unsigned R, unsigned C> inline
00391 void Matrix<R,C>::set_rows( const Matrix<1,C>* v ) {
00392   for( unsigned r = 0; r < R; ++r) set_row( r, v[r] );
00393 }
00394 
00395 template <unsigned R, unsigned C> inline
00396 void Matrix<R,C>::set_column( unsigned c, const Matrix<R,1>& v ) {
00397   for (unsigned i = 0; i < R; ++i) operator()(i,c) = v(i,0);
00398 }
00399 template <unsigned R, unsigned C> inline
00400 void Matrix<R,C>::add_column( unsigned c, const Matrix<R,1>& v ) {
00401   for (unsigned i = 0; i < R; ++i) operator()(i,c) += v(i,0);
00402 }
00403 template <unsigned R, unsigned C> inline
00404 void Matrix<R,C>::set_column_transpose( unsigned c, const Matrix<1,R>& v ) {
00405   for (unsigned i = 0; i < R; ++i) operator()(i,c) = v(0,i);
00406 }
00407 template <unsigned R, unsigned C> inline
00408 void Matrix<R,C>::set_columns( const Matrix<R,1>* v ) {
00409   for( unsigned c = 0; c < C; ++c) set_column( c, v[c] );
00410 }
00411 
00412 
00413 template <unsigned R, unsigned C> inline
00414 Matrix<R,1> Matrix<R,C>::column( unsigned c ) const {
00415   Matrix<R,1> result;
00416   for (unsigned i = 0; i < R; ++i)
00417     result(i,0) = operator()(i,c);
00418   return result;
00419 }
00420 
00421 template <unsigned R, unsigned C> inline
00422 Matrix<1,R> Matrix<R,C>::column_transpose( unsigned c ) const {
00423   Matrix<1,R> result;
00424   for (unsigned i = 0; i < R; ++i)
00425     result(0,i) = operator()(i,c);
00426   return result;
00427 }
00428 
00430 template <unsigned R1, unsigned C1, unsigned R2, unsigned C2> inline
00431 void set_region( Matrix<R1,C1>& d, unsigned r, unsigned c, Matrix<R2,C2>& s ) {
00432   const unsigned rmax = r+R2 > R1 ? R1 : r+R2;
00433   const unsigned cmax = c+C2 > C1 ? C1 : c+C2;
00434   for (unsigned i = r; i < rmax; ++i)
00435     for (unsigned j = c; j < cmax; ++j)
00436       d(i,j) = s(i-r,j-c);
00437 }
00438 
00439 
00440 template <unsigned R, unsigned C> inline
00441 Matrix<R,C>& Matrix<R,C>::assign_add_transpose( const Matrix<C,R>& other ) {
00442   for (unsigned r = 0; r < R; ++r)
00443     for (unsigned c = 0; c < C; ++c)
00444       operator()(r,c) += other(c,r);
00445   return *this;
00446 }
00447 
00448 template <unsigned R, unsigned C> inline
00449 Matrix<R,C>& Matrix<R,C>::assign_multiply_elements( const Matrix<R,C>& other ) {
00450   for (unsigned i = 0; i < R*C; ++i)  m[i] *= other.data()[i]; return *this;
00451 }
00452 
00453 template <unsigned R, unsigned C> inline
00454 Matrix<R,C>& Matrix<R,C>::assign_product( double s, const Matrix<R,C>& other ) {
00455   for (unsigned i = 0; i < R*C; ++i)  m[i] = s * other.data()[i]; return *this;
00456 }
00457 
00458 template <unsigned R, unsigned C> inline
00459 Matrix<R,C>& Matrix<R,C>::assign_add_product( double s, const Matrix<R,C>& other ) {
00460   for (unsigned i = 0; i < R*C; ++i) m[i] += s * other.data()[i]; return *this;
00461 }
00462 
00463 template <unsigned R, unsigned C> inline
00464 Matrix<R,C>& Matrix<R,C>::operator+=( const Matrix<R,C>& other ) {
00465   for (unsigned i = 0; i < R*C; ++i)  m[i] += other.data()[i]; return *this;
00466 }
00467 
00468 template <unsigned R, unsigned C> inline
00469 Matrix<R,C>& Matrix<R,C>::operator-=( const Matrix<R,C>& other ) {
00470   for (unsigned i = 0; i < R*C; ++i)  m[i] -= other.data()[i]; return *this;
00471 }
00472 
00473 template <unsigned R, unsigned C> inline
00474 Matrix<R,C>& Matrix<R,C>::operator+=( double s ) {
00475   for (unsigned i = 0; i < R*C; ++i)  m[i] += s; return *this;
00476 }
00477 
00478 template <unsigned R, unsigned C> inline
00479 Matrix<R,C>& Matrix<R,C>::operator-=( double s ) {
00480   for (unsigned i = 0; i < R*C; ++i)  m[i] -= s; return *this;
00481 }
00482 
00483 template <unsigned R, unsigned C> inline
00484 Matrix<R,C>& Matrix<R,C>::operator*=( double s ) {
00485   for (unsigned i = 0; i < R*C; ++i)  m[i] *= s; return *this;
00486 }
00487 
00488 template <unsigned R, unsigned C> inline
00489 Matrix<R,C>& Matrix<R,C>::operator/=( double s ) {
00490   for (unsigned i = 0; i < R*C; ++i)  m[i] /= s; return *this;
00491 }
00492 
00493 
00494 template <unsigned R, unsigned C> inline
00495 Matrix<R,C> operator-( const Matrix<R,C>& m ) {
00496   Matrix<R,C> result;
00497   for (unsigned i = 0; i < R*C; ++i)
00498     result.data()[i] = -m.data()[i];
00499   return result;
00500 }
00501 
00502 template <unsigned R, unsigned C> inline
00503 Matrix<R,C> operator+( const Matrix<R,C>& m, double s ) {
00504   Matrix<R,C> tmp(m); tmp += s; return tmp;
00505 }
00506 
00507 template <unsigned R, unsigned C> inline
00508 Matrix<R,C> operator+( double s, const Matrix<R,C>& m ) {
00509   Matrix<R,C> tmp(m); tmp += s; return tmp;
00510 }
00511 
00512 template <unsigned R, unsigned C> inline
00513 Matrix<R,C> operator+( const Matrix<R,C>& A, const Matrix<R,C>& B ) {
00514   Matrix<R,C> tmp(A); tmp += B; return tmp;
00515 }
00516 
00517 
00518 template <unsigned R, unsigned C> inline
00519 Matrix<R,C> operator-( const Matrix<R,C>& m, double s ) {
00520   Matrix<R,C> tmp(m); tmp -= s; return tmp;
00521 }
00522 
00523 template <unsigned R, unsigned C> inline
00524 Matrix<R,C> operator-( double s, const Matrix<R,C>& m ) {
00525   Matrix<R,C> tmp(m); tmp -= s; return tmp;
00526 }
00527 
00528 template <unsigned R, unsigned C> inline
00529 Matrix<R,C> operator-( const Matrix<R,C>& A, const Matrix<R,C>& B ) {
00530   Matrix<R,C> tmp(A); tmp -= B; return tmp;
00531 }
00532 
00533 
00534 template <unsigned R, unsigned C> inline
00535 Matrix<R,C> operator*( const Matrix<R,C>& m, double s ) {
00536   Matrix<R,C> tmp(m); tmp *= s; return tmp;
00537 }
00538 
00539 template <unsigned R, unsigned C> inline
00540 Matrix<R,C> operator*( double s, const Matrix<R,C>& m ) {
00541   Matrix<R,C> tmp(m); tmp *= s; return tmp;
00542 }
00543 
00544 template <unsigned R, unsigned RC, unsigned C> inline
00545 double multiply_helper_result_val( unsigned r, unsigned c,
00546                                    const Matrix<R,RC>& A,
00547                                    const Matrix<RC,C>& B ) {
00548   double tmp = A(r,0)*B(0,c);
00549   switch (RC) {
00550     default: for (unsigned k = 6; k < RC; ++k)
00551             tmp += A(r,k)*B(k,c);
00552     case 6: tmp += A(r,5)*B(5,c);
00553     case 5: tmp += A(r,4)*B(4,c);
00554     case 4: tmp += A(r,3)*B(3,c);
00555     case 3: tmp += A(r,2)*B(2,c);
00556     case 2: tmp += A(r,1)*B(1,c);
00557     case 1: ;
00558   }
00559   return tmp;
00560 }
00561 
00562 template <unsigned R, unsigned RC, unsigned C> inline
00563 Matrix<R,C> operator*( const Matrix<R,RC>& A, const Matrix<RC,C>& B ) {
00564   //Matrix<R,C> result(0.0);
00565   //for (unsigned i = 0; i < R; ++i)
00566   //  for (unsigned j = 0; j < C; ++j)
00567   //    for (unsigned k = 0; k < RC; ++k)
00568   //      result(i,j) += A(i,k) * B(k,j);
00569 
00570   Matrix<R,C> result;
00571   for (unsigned r = 0; r < R; ++r)
00572     for (unsigned c = 0; c < C; ++c)
00573       result(r,c) = multiply_helper_result_val( r, c, A, B );
00574 
00575   return result;
00576 }
00577 
00578 
00579 template <unsigned R, unsigned C> inline
00580 Matrix<R,C> operator/( const Matrix<R,C>& m, double s ) {
00581   Matrix<R,C> tmp(m); tmp /= s; return tmp;
00582 }
00583 
00584 template <unsigned RC> inline
00585 double cofactor( const Matrix<RC,RC>& m, unsigned r, unsigned c ) {
00586   const double sign[] = { 1.0, -1.0 };
00587   return sign[(r+c)%2] * det( Matrix<RC-1,RC-1>( m, r, c ) );
00588 }
00589 
00590 template <unsigned RC> inline
00591 double det( const Matrix<RC,RC>& m ) {
00592   double result = 0.0;
00593   for (unsigned i = 0; i < RC; ++i)
00594     result += m(0,i) * cofactor<RC>( m, 0, i );
00595   return result;
00596 }
00597 
00598 //inline
00599 //double det( const Matrix<1,1>& m ) {
00600 // return m(0,0);
00601 //}
00602 
00603 inline
00604 double det( const Matrix<2,2>& m ) {
00605   return m(0,0)*m(1,1) - m(0,1)*m(1,0);
00606 }
00607 
00608 inline
00609 double det( const Matrix<3,3>& m ) {
00610   return m(0,0)*(m(1,1)*m(2,2) - m(2,1)*m(1,2))
00611          + m(0,1)*(m(2,0)*m(1,2) - m(1,0)*m(2,2))
00612          + m(0,2)*(m(1,0)*m(2,1) - m(2,0)*m(1,1));
00613 }
00614 
00615 inline
00616 Matrix<2,2> adj( const Matrix<2,2>& m ) {
00617   Matrix<2,2> result;
00618   result(0,0) =  m(1,1);
00619   result(0,1) = -m(0,1);
00620   result(1,0) = -m(1,0);
00621   result(1,1) =  m(0,0);
00622   return result;
00623 }
00624 
00625 inline
00626 Matrix<2,2> transpose_adj( const Matrix<2,2>& m ) {
00627   Matrix<2,2> result;
00628   result(0,0) =  m(1,1);
00629   result(1,0) = -m(0,1);
00630   result(0,1) = -m(1,0);
00631   result(1,1) =  m(0,0);
00632   return result;
00633 }
00634 
00635 inline
00636 Matrix<3,3> adj( const Matrix<3,3>& m ) {
00637   Matrix<3,3> result;
00638 
00639   result(0,0) = m(1,1)*m(2,2) - m(1,2)*m(2,1);
00640   result(0,1) = m(0,2)*m(2,1) - m(0,1)*m(2,2);
00641   result(0,2) = m(0,1)*m(1,2) - m(0,2)*m(1,1);
00642 
00643   result(1,0) = m(1,2)*m(2,0) - m(1,0)*m(2,2);
00644   result(1,1) = m(0,0)*m(2,2) - m(0,2)*m(2,0);
00645   result(1,2) = m(0,2)*m(1,0) - m(0,0)*m(1,2);
00646 
00647   result(2,0) = m(1,0)*m(2,1) - m(1,1)*m(2,0);
00648   result(2,1) = m(0,1)*m(2,0) - m(0,0)*m(2,1);
00649   result(2,2) = m(0,0)*m(1,1) - m(0,1)*m(1,0);
00650 
00651   return result;
00652 }
00653 
00654 inline
00655 Matrix<3,3> transpose_adj( const Matrix<3,3>& m ) {
00656   Matrix<3,3> result;
00657 
00658   result(0,0) = m(1,1)*m(2,2) - m(1,2)*m(2,1);
00659   result(0,1) = m(1,2)*m(2,0) - m(1,0)*m(2,2);
00660   result(0,2) = m(1,0)*m(2,1) - m(1,1)*m(2,0);
00661 
00662   result(1,0) = m(0,2)*m(2,1) - m(0,1)*m(2,2);
00663   result(1,1) = m(0,0)*m(2,2) - m(0,2)*m(2,0);
00664   result(1,2) = m(0,1)*m(2,0) - m(0,0)*m(2,1);
00665 
00666   result(2,0) = m(0,1)*m(1,2) - m(0,2)*m(1,1);
00667   result(2,1) = m(0,2)*m(1,0) - m(0,0)*m(1,2);
00668   result(2,2) = m(0,0)*m(1,1) - m(0,1)*m(1,0);
00669 
00670   return result;
00671 }
00672 
00673 template <unsigned R, unsigned C> inline
00674 Matrix<R,C> inverse( const Matrix<R,C>& m ) {
00675   return adj(m) * (1.0 / det(m));
00676 }
00677 
00678 template <unsigned R, unsigned C> inline
00679 Matrix<R,C> transpose( const Matrix<C,R>& m ) {
00680   Matrix<R,C> result;
00681   for (unsigned r = 0; r < R; ++r)
00682     for (unsigned c = 0; c < C; ++c)
00683       result(r,c) = m(c,r);
00684   return result;
00685 }
00686 /*
00687 template <unsigned R> inline
00688 const Matrix<R,1>& transpose( const Matrix<1,R>& m ) {
00689   return *reinterpret_cast<const Matrix<R,1>*>(&m);
00690 }
00691 
00692 template <unsigned C> inline
00693 const Matrix<1,C>& transpose( const Matrix<C,1>& m ) {
00694   return *reinterpret_cast<const Matrix<1,C>*>(&m);
00695 }
00696 */
00697 template <unsigned RC> inline
00698 double trace( const Matrix<RC,RC>& m ) {
00699   double result = m(0,0);
00700   for (unsigned i = 1; i < RC; ++i)
00701     result += m(i,i);
00702   return result;
00703 }
00704 
00705 template <unsigned R, unsigned C> inline
00706 double sqr_Frobenius( const Matrix<R,C>& m ) {
00707   double result = *m.data() * *m.data();
00708   for (unsigned i = 1; i < R*C; ++i)
00709     result += m.data()[i] * m.data()[i];
00710   return result;
00711 }
00712 
00713 template <unsigned R, unsigned C> inline
00714 double Frobenius( const Matrix<R,C>& m ) {
00715   return std::sqrt( sqr_Frobenius<R,C>( m ) );
00716 }
00717 
00718 template <unsigned R, unsigned C> inline
00719 bool operator==( const Matrix<R,C>& A, const Matrix<R,C>& B ) {
00720   for (unsigned i = 0; i < R*C; ++i)
00721     if (A.data()[i] != B.data()[i])
00722       return false;
00723   return true;
00724 }
00725 
00726 template <unsigned R, unsigned C> inline
00727 bool operator!=( const Matrix<R,C>& A, const Matrix<R,C>& B ) {
00728   return !(A == B);
00729 }
00730 
00731 template <unsigned R, unsigned C> inline
00732 std::ostream& operator<<( std::ostream& str, const Matrix<R,C>& m ) {
00733   str << m.data()[0];
00734   for (unsigned i = 1; i < R*C; ++i)
00735     str << ' ' << m.data()[i];
00736   return str;
00737 }
00738 
00739 template <unsigned R, unsigned C> inline
00740 std::istream& operator>>( std::istream& str, Matrix<R,C>& m ) {
00741   for (unsigned i = 0; i < R*C; ++i)
00742     str >> m.data()[i];
00743   return str;
00744 }
00745 
00746 template <unsigned R> inline
00747 double sqr_length( const Matrix<R,1>& v ) {
00748   return sqr_Frobenius(v);
00749 }
00750 
00751 template <unsigned C> inline
00752 double sqr_length( const Matrix<1,C>& v ) {
00753   return sqr_Frobenius(v);
00754 }
00755 
00756 template <unsigned R> inline
00757 double length( const Matrix<R,1>& v ) {
00758   return Frobenius(v);
00759 }
00760 
00761 template <unsigned C> inline
00762 double length( const Matrix<1,C>& v ) {
00763   return Frobenius(v);
00764 }
00765 
00766 template <unsigned R, unsigned C> inline
00767 double inner_product( const Matrix<R,C>& m1, const Matrix<R,C>& m2 ) {
00768   double result = 0.0;
00769   for (unsigned r = 0; r < R; ++r)
00770     for (unsigned c = 0; c < C; ++c)
00771       result += m1(r,c) * m2(r,c);
00772   return result;
00773 }
00774 
00775 template <unsigned R, unsigned C> inline
00776 Matrix<R,C> outer( const Matrix<R,1>& v1, const Matrix<C,1>& v2 ) {
00777   Matrix<R,C> result;
00778   for (unsigned r = 0; r < R; ++r)
00779     for (unsigned c = 0; c < C; ++c)
00780       result(r,c) = v1(r,0) * v2(c,0);
00781   return result;
00782 }
00783 
00784 template <unsigned R, unsigned C> inline
00785 Matrix<R,C> outer( const Matrix<1,R>& v1, const Matrix<1,C>& v2 ) {
00786   Matrix<R,C> result;
00787   for (unsigned r = 0; r < R; ++r)
00788     for (unsigned c = 0; c < C; ++c)
00789       result(r,c) = v1(0,r) * v2(0,c);
00790   return result;
00791 }
00792 
00793 inline
00794 double vector_product( const Matrix<2,1>& v1, const Matrix<2,1>& v2 ) {
00795   return v1(0,0) * v2(1,0) - v1(1,0) * v2(0,0);
00796 }
00797 
00798 inline
00799 double vector_product( const Matrix<1,2>& v1, const Matrix<1,2>& v2 ) {
00800   return v1(0,0) * v2(0,1) - v1(0,1) * v2(0,0);
00801 }
00802 
00803 inline
00804 Matrix<3,1> vector_product( const Matrix<3,1>& a, const Matrix<3,1>& b ) {
00805   Matrix<3,1> result;
00806   result(0,0) = a(1,0)*b(2,0) - a(2,0)*b(1,0);
00807   result(1,0) = a(2,0)*b(0,0) - a(0,0)*b(2,0);
00808   result(2,0) = a(0,0)*b(1,0) - a(1,0)*b(0,0);
00809   return result;
00810 }
00811 
00812 inline
00813 Matrix<1,3> vector_product( const Matrix<1,3>& a, const Matrix<1,3>& b ) {
00814   Matrix<1,3> result;
00815   result(0,0) = a(0,1)*b(0,2) - a(0,2)*b(0,1);
00816   result(0,1) = a(0,2)*b(0,0) - a(0,0)*b(0,2);
00817   result(0,2) = a(0,0)*b(0,1) - a(0,1)*b(0,0);
00818   return result;
00819 }
00820 
00821 template <unsigned R, unsigned C> inline
00822 double operator%( const Matrix<R,C>& v1, const Matrix<R,C>& v2 ) {
00823   return inner_product( v1, v2 );
00824 }
00825 
00826 inline
00827 double operator*( const Matrix<2,1>& v1, const Matrix<2,1>& v2 ) {
00828   return vector_product( v1, v2 );
00829 }
00830 
00831 inline
00832 double operator*( const Matrix<1,2>& v1, const Matrix<1,2>& v2 ) {
00833   return vector_product( v1, v2 );
00834 }
00835 
00836 inline
00837 Matrix<3,1> operator*( const Matrix<3,1>& v1, const Matrix<3,1>& v2 ) {
00838   return vector_product( v1, v2 );
00839 }
00840 
00841 inline
00842 Matrix<1,3> operator*( const Matrix<1,3>& v1, const Matrix<1,3>& v2 ) {
00843  return vector_product( v1, v2 );
00844 }
00845 
00847 inline
00848 void QR( const Matrix<3,3>& A, Matrix<3,3>& Q, Matrix<3,3>& R ) {
00849   Q = A;
00850 
00851   R(0,0) = sqrt(Q(0,0)*Q(0,0) + Q(1,0)*Q(1,0) + Q(2,0)*Q(2,0));
00852   double temp_dbl = 1.0/R(0,0);
00853   R(1,0) = 0.0;
00854   R(2,0) = 0.0;
00855   Q(0,0) *= temp_dbl;
00856   Q(1,0) *= temp_dbl;
00857   Q(2,0) *= temp_dbl;
00858 
00859 
00860   R(0,1)  = Q(0,0)*Q(0,1) + Q(1,0)*Q(1,1) + Q(2,0)*Q(2,1);
00861   Q(0,1) -= Q(0,0)*R(0,1);
00862   Q(1,1) -= Q(1,0)*R(0,1);
00863   Q(2,1) -= Q(2,0)*R(0,1);
00864 
00865   R(0,2)  = Q(0,0)*Q(0,2) + Q(1,0)*Q(1,2) + Q(2,0)*Q(2,2);
00866   Q(0,2) -= Q(0,0)*R(0,2);
00867   Q(1,2) -= Q(1,0)*R(0,2);
00868   Q(2,2) -= Q(2,0)*R(0,2);
00869 
00870   R(1,1) = sqrt(Q(0,1)*Q(0,1) + Q(1,1)*Q(1,1) + Q(2,1)*Q(2,1));
00871   temp_dbl = 1.0 / R(1,1);
00872   R(2,1) = 0.0;
00873   Q(0,1) *= temp_dbl;
00874   Q(1,1) *= temp_dbl;
00875   Q(2,1) *= temp_dbl;
00876 
00877 
00878   R(1,2)  = Q(0,1)*Q(0,2) + Q(1,1)*Q(1,2) + Q(2,1)*Q(2,2);
00879   Q(0,2) -= Q(0,1)*R(1,2);
00880   Q(1,2) -= Q(1,1)*R(1,2);
00881   Q(2,2) -= Q(2,1)*R(1,2);
00882 
00883   R(2,2) = sqrt(Q(0,2)*Q(0,2) + Q(1,2)*Q(1,2) + Q(2,2)*Q(2,2));
00884   temp_dbl = 1.0 / R(2,2);
00885   Q(0,2) *= temp_dbl;
00886   Q(1,2) *= temp_dbl;
00887   Q(2,2) *= temp_dbl;
00888 }
00889 
00890 
00892 inline
00893 void QR( const Matrix<2,2>& A, Matrix<2,2>& Q, Matrix<2,2>& R ) {
00894   R(0,0) = std::sqrt( A(0,0)*A(0,0) + A(1,0)*A(1,0) );
00895   const double r0inv = 1.0 / R(0,0);
00896   Q(0,0) = Q(1,1) = A(0,0) * r0inv;
00897   Q(1,0) = A(1,0) * r0inv;
00898   Q(0,1) = -Q(1,0);
00899   R(0,1) = r0inv * (A(0,0)*A(0,1) + A(1,0)*A(1,1));
00900   R(1,1) = r0inv * (A(0,0)*A(1,1) - A(0,1)*A(1,0));
00901   R(1,0) = 0.0;
00902 }
00903 
00904 } // namespace MeshKit
00905 
00906 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines