ESyS-Particle  4.0.1
Matrix3.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 #ifndef __MATRIX3_HPP
00014 #define __MATRIX3_HPP
00015 
00016 #include "Foundation/Matrix3.h"
00017 #include "Foundation/vec3.h"
00018 
00019 
00020 MATRIX3_INLINE Matrix3::Matrix3()
00021 {
00022   for(int i=0;i<3;i++)
00023   {
00024     for(int j=0;j<3;j++)
00025     {
00026       m[i][j]=0;
00027     }
00028   }
00029 }
00030 
00031 MATRIX3_INLINE Matrix3::Matrix3(const Vec3& V1, const Vec3& V2,const Vec3& V3)
00032 {
00033   m[0][0]=V1.data[0];
00034   m[0][1]=V2.data[0];
00035   m[0][2]=V3.data[0];
00036   m[1][0]=V1.data[1];
00037   m[1][1]=V2.data[1];
00038   m[1][2]=V3.data[1];
00039   m[2][0]=V1.data[2];
00040   m[2][1]=V2.data[2];
00041   m[2][2]=V3.data[2];
00042   
00043 }
00044 
00045 MATRIX3_INLINE Matrix3::Matrix3(const double a[3][3])
00046 {
00047   m[0][0]=a[0][0];
00048   m[0][1]=a[0][1];
00049   m[0][2]=a[0][2];
00050   m[1][0]=a[1][0];
00051   m[1][1]=a[1][1];
00052   m[1][2]=a[1][2];
00053   m[2][0]=a[2][0];
00054   m[2][1]=a[2][1];
00055   m[2][2]=a[2][2];
00056 }
00057 
00058 MATRIX3_INLINE Matrix3::Matrix3(const Matrix3& rhs)
00059 {
00060   for(int i=0;i<3;i++)
00061   {
00062     for(int j=0;j<3;j++)
00063     {
00064       m[i][j]=rhs.m[i][j];
00065     }
00066   }
00067 }
00068 
00069 MATRIX3_INLINE Matrix3::~Matrix3()
00070 {
00071 }
00072 
00076 MATRIX3_INLINE double Matrix3::det()
00077 {
00078   return m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])+
00079     m[0][1]*(m[1][2]*m[2][0]-m[1][0]*m[2][2])+
00080     m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
00081 }
00082 
00083 
00084 MATRIX3_INLINE Matrix3 Matrix3::inv()
00085 {
00086   Matrix3 res=*this;
00087 
00088   res.invert();
00089   
00090   return res;
00091 }
00092 
00093 MATRIX3_INLINE void Matrix3::transpose()
00094 {
00095   double h;
00096 
00097   h=m[1][0];
00098   m[1][0]=m[0][1];
00099   m[0][1]=h;
00100   h=m[2][0];
00101   m[2][0]=m[0][2];
00102   m[0][2]=h;
00103   h=m[2][1];
00104   m[2][1]=m[1][2];
00105   m[1][2]=h;
00106 }
00107 
00108 MATRIX3_INLINE Matrix3 Matrix3::trans() const
00109 {
00110   Matrix3 res;
00111   
00112   res.m[0][0]=m[0][0];
00113   res.m[0][1]=m[1][0];
00114   res.m[0][2]=m[2][0];
00115   res.m[1][0]=m[0][1];
00116   res.m[1][1]=m[1][1];
00117   res.m[1][2]=m[2][1];
00118   res.m[2][0]=m[0][2];
00119   res.m[2][1]=m[1][2];
00120   res.m[2][2]=m[2][2];
00121 
00122   return res;
00123 }
00124 
00125 // 15 Flops ( 9 mult,6 add)
00126 MATRIX3_INLINE Vec3 Matrix3::operator *(const Vec3& V) const
00127 {
00128   double x=m[0][0]*V.data[0]+m[0][1]*V.data[1]+m[0][2]*V.data[2];
00129   double y=m[1][0]*V.data[0]+m[1][1]*V.data[1]+m[1][2]*V.data[2];
00130   double z=m[2][0]*V.data[0]+m[2][1]*V.data[1]+m[2][2]*V.data[2];
00131   
00132   return Vec3(x,y,z);
00133 }
00134 
00135 // 9 Flops
00136 MATRIX3_INLINE Matrix3 Matrix3::operator *(double d) const
00137 {
00138   Matrix3 res;
00139 
00140   res.m[0][0]=m[0][0]*d;
00141   res.m[0][1]=m[0][1]*d;
00142   res.m[0][2]=m[0][2]*d;
00143   res.m[1][0]=m[1][0]*d;
00144   res.m[1][1]=m[1][1]*d;
00145   res.m[1][2]=m[1][2]*d;
00146   res.m[2][0]=m[2][0]*d;
00147   res.m[2][1]=m[2][1]*d;
00148   res.m[2][2]=m[2][2]*d;
00149 
00150   return res;
00151 }
00152 
00153 // 9 Flops
00154 MATRIX3_INLINE Matrix3 Matrix3::operator /(double d) const
00155 {
00156   Matrix3 res;
00157 
00158   res.m[0][0]=m[0][0]/d;
00159   res.m[0][1]=m[0][1]/d;
00160   res.m[0][2]=m[0][2]/d;
00161   res.m[1][0]=m[1][0]/d;
00162   res.m[1][1]=m[1][1]/d;
00163   res.m[1][2]=m[1][2]/d;
00164   res.m[2][0]=m[2][0]/d;
00165   res.m[2][1]=m[2][1]/d;
00166   res.m[2][2]=m[2][2]/d;
00167 
00168   return res;
00169 }
00170 
00171 MATRIX3_INLINE Matrix3 &Matrix3::operator=(const Matrix3 &other)
00172 {
00173   for(int i = 0; i < 3; i++) {
00174     for(int j = 0; j < 3; j++) {
00175       m[i][j] = other.m[i][j];
00176     }
00177   }
00178   return *this;
00179 }
00180 
00181 MATRIX3_INLINE bool Matrix3::operator==(const Matrix3 &other) const
00182 {
00183   for(int i = 0; i < 3; i++) {
00184     for(int j = 0; j < 3; j++) {
00185       if (m[i][j] != other.m[i][j]) {
00186         return false;
00187       }
00188     }
00189   }
00190   return true;
00191 }
00192 
00193 MATRIX3_INLINE std::ostream &operator<<(std::ostream &oStream, const Matrix3 &m)
00194 {
00195   oStream << m(0,0);
00196   for(int i = 1; i < 9; i++) {
00197     oStream << " " << m(i/3, i%3);
00198   }
00199   return oStream;
00200 }
00201 
00202 // 45 Flops (27mult, 18add)
00203 MATRIX3_INLINE Matrix3 Matrix3::operator *(const Matrix3& rhs) const
00204 {
00205   Matrix3 res;
00206 
00207   for(int i=0;i<3;i++){
00208     for(int j=0;j<3;j++){
00209       res.m[i][j]=m[i][0]*rhs.m[0][j]+m[i][1]*rhs.m[1][j]+m[i][2]*rhs.m[2][j];
00210     }
00211   }
00212   return res;
00213 }
00214 
00215 // 9 Flops
00216 MATRIX3_INLINE Matrix3& Matrix3::operator+=(const Matrix3& rhs)
00217 {
00218   for(int i=0;i<3;i++){
00219     for(int j=0;j<3;j++){
00220       m[i][j]=m[i][j]+rhs.m[i][j];
00221     }
00222   }
00223   return *this;
00224 }
00225 
00229 MATRIX3_INLINE Matrix3 Matrix3::operator +(const Matrix3& rhs) const
00230 {
00231   Matrix3 res;
00232   
00233   for(int i=0;i<3;i++){
00234     for(int j=0;j<3;j++){
00235       res.m[i][j]=m[i][j]+rhs.m[i][j];
00236     }
00237   }
00238 
00239   return res;
00240 }
00241 
00245 MATRIX3_INLINE Matrix3 Matrix3::operator-(const Matrix3& rhs) const
00246 {
00247   Matrix3 res;
00248   
00249   for(int i=0;i<3;i++) {
00250     for(int j=0;j<3;j++) {
00251       res.m[i][j]=m[i][j]-rhs.m[i][j];
00252     }
00253   }
00254 
00255   return res;
00256 }
00257 
00261 MATRIX3_INLINE double Matrix3::trace() const
00262 {
00263   return m[0][0]+m[1][1]+m[2][2];
00264 }
00265 
00269 MATRIX3_INLINE double Matrix3::norm() const
00270 {
00271   double res=0.0;
00272 
00273   for(int i=0;i<3;i++){
00274     for(int j=0;j<3;j++){
00275       res+=m[i][j]*m[i][j];
00276     }
00277   }
00278   return res;
00279 }
00280 
00281 // matrix from vector
00282 MATRIX3_INLINE Matrix3 star(const Vec3& V)
00283 {
00284   Matrix3 res;
00285 
00286   res.m[0][1]=-V.Z();
00287   res.m[0][2]=V.Y();
00288   res.m[1][0]=V.Z();
00289   res.m[1][2]=-V.X();
00290   res.m[2][0]=-V.Y();
00291   res.m[2][1]=V.X();
00292 
00293   return res;
00294 }
00295 
00296 // unit matrix
00297 MATRIX3_INLINE Matrix3 Matrix3::Unit()
00298 {
00299   Matrix3 res;
00300 
00301   res.m[0][0]=1.0;
00302   res.m[1][1]=1.0;
00303   res.m[2][2]=1.0;
00304 
00305   return res;
00306 }
00307 
00311 MATRIX3_INLINE Matrix3 operator*(double d,const Matrix3& M)
00312 {
00313   Matrix3 res;
00314 
00315   for(int i=0;i<3;i++){
00316     for(int j=0;j<3;j++){
00317       res.m[i][j]=d*M.m[i][j];
00318     }
00319   }
00320 
00321   return res;
00322 }
00323 
00324 
00325 
00326 #endif // __MATRIX3_HPP