Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
Public Member Functions | Static Public Attributes | Private Attributes | Friends | List of all members
moab::Matrix3 Class Reference

#include <Matrix3.hpp>

Public Member Functions

 Matrix3 ()
 
 Matrix3 (Eigen::Matrix3d mat)
 
 Matrix3 (double diagonal)
 
 Matrix3 (const CartVect &diagonal)
 
 Matrix3 (const std::vector< double > &diagonal)
 
 Matrix3 (double v00, double v01, double v02, double v10, double v11, double v12, double v20, double v21, double v22)
 
 Matrix3 (const Matrix3 &f)
 
template<typename Vector >
 Matrix3 (const Vector &row0, const Vector &row1, const Vector &row2, const bool isRow)
 
 Matrix3 (const double v[size])
 
void copyto (double v[Matrix3::size])
 
Matrix3 & operator= (const Matrix3 &m)
 
Matrix3 & operator= (const double v[size])
 
double * operator[] (unsigned i)
 
const double * operator[] (unsigned i) const
 
double & operator() (unsigned r, unsigned c)
 
double operator() (unsigned r, unsigned c) const
 
double & operator() (unsigned i)
 
double operator() (unsigned i) const
 
double * array ()
 
const double * array () const
 
Matrix3 & operator+= (const Matrix3 &m)
 
Matrix3 & operator-= (const Matrix3 &m)
 
Matrix3 & operator*= (double s)
 
Matrix3 & operator/= (double s)
 
Matrix3 & operator*= (const Matrix3 &m)
 
bool is_symmetric ()
 
bool is_positive_definite ()
 
template<typename Vector >
ErrorCode eigen_decomposition (Vector &evals, Matrix3 &evecs)
 
void transpose_inplace ()
 
Matrix3 transpose () const
 
template<typename Vector >
void copycol (int index, Vector &vol)
 
void swapcol (int srcindex, int destindex)
 
template<typename Vector >
Vector vcol (int index) const
 
void colscale (int index, double scale)
 
void rowscale (int index, double scale)
 
CartVect col (int index) const
 
CartVect row (int index) const
 
double determinant () const
 
Matrix3 inverse () const
 
bool invert ()
 
double subdet (int r, int c) const
 
void print (std::ostream &s) const
 

Static Public Attributes

static const int size = 9
 

Private Attributes

Eigen::Matrix3d _mat
 

Friends

Matrix3 operator+ (const Matrix3 &a, const Matrix3 &b)
 
Matrix3 operator- (const Matrix3 &a, const Matrix3 &b)
 
Matrix3 operator* (const Matrix3 &a, const Matrix3 &b)
 

Detailed Description

Definition at line 216 of file Matrix3.hpp.

Constructor & Destructor Documentation

◆ Matrix3() [1/9]

moab::Matrix3::Matrix3 ( )
inline

Definition at line 231 of file Matrix3.hpp.

232  {
233 #ifndef MOAB_HAVE_LAPACK
234  _mat.fill( 0.0 );
235 #else
236  MOAB_DMEMZERO( _mat, Matrix3::size );
237 #endif
238  }

References _mat, and size.

Referenced by inverse(), and transpose().

◆ Matrix3() [2/9]

moab::Matrix3::Matrix3 ( Eigen::Matrix3d  mat)
inline

Definition at line 241 of file Matrix3.hpp.

241 : _mat( mat ) {}

◆ Matrix3() [3/9]

moab::Matrix3::Matrix3 ( double  diagonal)
inline

Definition at line 246 of file Matrix3.hpp.

247  {
248 #ifndef MOAB_HAVE_LAPACK
249  _mat << diagonal, 0.0, 0.0, 0.0, diagonal, 0.0, 0.0, 0.0, diagonal;
250 #else
251  MOAB_DMEMZERO( _mat, Matrix3::size );
252  _mat[0] = _mat[4] = _mat[8] = diagonal;
253 #endif
254  }

References _mat, and size.

◆ Matrix3() [4/9]

moab::Matrix3::Matrix3 ( const CartVect &  diagonal)
inline

Definition at line 256 of file Matrix3.hpp.

257  {
258 #ifndef MOAB_HAVE_LAPACK
259  _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
260 #else
261  MOAB_DMEMZERO( _mat, Matrix3::size );
262  _mat[0] = diagonal[0];
263  _mat[4] = diagonal[1];
264  _mat[8] = diagonal[2];
265 #endif
266  }

References _mat, and size.

◆ Matrix3() [5/9]

moab::Matrix3::Matrix3 ( const std::vector< double > &  diagonal)
inline

Definition at line 273 of file Matrix3.hpp.

274  {
275 #ifndef MOAB_HAVE_LAPACK
276  _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
277 #else
278  MOAB_DMEMZERO( _mat, Matrix3::size );
279  _mat[0] = diagonal[0];
280  _mat[4] = diagonal[1];
281  _mat[8] = diagonal[2];
282 #endif
283  }

References _mat, and size.

◆ Matrix3() [6/9]

moab::Matrix3::Matrix3 ( double  v00,
double  v01,
double  v02,
double  v10,
double  v11,
double  v12,
double  v20,
double  v21,
double  v22 
)
inline

Definition at line 285 of file Matrix3.hpp.

294  {
295 #ifndef MOAB_HAVE_LAPACK
296  _mat << v00, v01, v02, v10, v11, v12, v20, v21, v22;
297 #else
298  MOAB_DMEMZERO( _mat, Matrix3::size );
299  _mat[0] = v00;
300  _mat[1] = v01;
301  _mat[2] = v02;
302  _mat[3] = v10;
303  _mat[4] = v11;
304  _mat[5] = v12;
305  _mat[6] = v20;
306  _mat[7] = v21;
307  _mat[8] = v22;
308 #endif
309  }

References _mat, and size.

◆ Matrix3() [7/9]

moab::Matrix3::Matrix3 ( const Matrix3 &  f)
inline

Definition at line 312 of file Matrix3.hpp.

313  {
314 #ifndef MOAB_HAVE_LAPACK
315  _mat = f._mat;
316 #else
317  memcpy( _mat, f._mat, size * sizeof( double ) );
318 #endif
319  }

References _mat, and size.

◆ Matrix3() [8/9]

template<typename Vector >
moab::Matrix3::Matrix3 ( const Vector &  row0,
const Vector &  row1,
const Vector &  row2,
const bool  isRow 
)
inline

Definition at line 323 of file Matrix3.hpp.

324  {
325 #ifndef MOAB_HAVE_LAPACK
326  if( isRow )
327  {
328  _mat << row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], row2[0], row2[1], row2[2];
329  }
330  else
331  {
332  _mat << row0[0], row1[0], row2[0], row0[1], row1[1], row2[1], row0[2], row1[2], row2[2];
333  }
334 #else
335  MOAB_DMEMZERO( _mat, Matrix3::size );
336  if( isRow )
337  {
338  _mat[0] = row0[0];
339  _mat[1] = row0[1];
340  _mat[2] = row0[2];
341  _mat[3] = row1[0];
342  _mat[4] = row1[1];
343  _mat[5] = row1[2];
344  _mat[6] = row2[0];
345  _mat[7] = row2[1];
346  _mat[8] = row2[2];
347  }
348  else
349  {
350  _mat[0] = row0[0];
351  _mat[1] = row1[0];
352  _mat[2] = row2[0];
353  _mat[3] = row0[1];
354  _mat[4] = row1[1];
355  _mat[5] = row2[1];
356  _mat[6] = row0[2];
357  _mat[7] = row1[2];
358  _mat[8] = row2[2];
359  }
360 #endif
361  }

References _mat, and size.

◆ Matrix3() [9/9]

moab::Matrix3::Matrix3 ( const double  v[size])
inline

Definition at line 376 of file Matrix3.hpp.

377  {
378 #ifndef MOAB_HAVE_LAPACK
379  _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
380 #else
381  memcpy( _mat, v, size * sizeof( double ) );
382 #endif
383  }

References _mat, and size.

Member Function Documentation

◆ array() [1/2]

double* moab::Matrix3::array ( )
inline

Definition at line 469 of file Matrix3.hpp.

470  {
471 #ifndef MOAB_HAVE_LAPACK
472  return _mat.data();
473 #else
474  return _mat;
475 #endif
476  }

References _mat.

Referenced by moab::EvalSet::evaluate_reverse(), and moab::LinearTet::evaluate_reverse().

◆ array() [2/2]

const double* moab::Matrix3::array ( ) const
inline

Definition at line 478 of file Matrix3.hpp.

479  {
480 #ifndef MOAB_HAVE_LAPACK
481  return _mat.data();
482 #else
483  return _mat;
484 #endif
485  }

References _mat.

◆ col()

CartVect moab::Matrix3::col ( int  index) const
inline

Definition at line 886 of file Matrix3.hpp.

887  {
888  assert( index < Matrix3::size );
889 #ifndef MOAB_HAVE_LAPACK
890  Eigen::Vector3d mvec = _mat.col( index );
891  return CartVect( mvec[0], mvec[1], mvec[2] );
892 #else
893  switch( index )
894  {
895  case 0:
896  return CartVect( _mat[0], _mat[3], _mat[6] );
897  case 1:
898  return CartVect( _mat[1], _mat[4], _mat[7] );
899  case 2:
900  return CartVect( _mat[2], _mat[5], _mat[8] );
901  }
902  return CartVect( 0.0 );
903 #endif
904  }

References _mat, and size.

Referenced by moab::OrientedBox::area(), moab::OrientedBox::axis(), moab::box_from_axes(), moab::OrientedBox::closest_location_in_box(), moab::OrientedBox::contained(), moab::OrientedBox::dimensions(), moab::OrientedBox::inner_radius(), moab::OrientedBox::inner_radius_squared(), moab::OrientedBox::make_hex(), moab::operator<<(), moab::OrientedBox::OrientedBox(), moab::OrientedBox::outer_radius(), moab::OrientedBox::outer_radius_squared(), moab::OrientedBox::scaled_axis(), and moab::OrientedBox::volume().

◆ colscale()

void moab::Matrix3::colscale ( int  index,
double  scale 
)
inline

Definition at line 832 of file Matrix3.hpp.

833  {
834  assert( index < Matrix3::size );
835 #ifndef MOAB_HAVE_LAPACK
836  _mat.col( index ) *= scale;
837 #else
838  switch( index )
839  {
840  case 0:
841  _mat[0] *= scale;
842  _mat[3] *= scale;
843  _mat[6] *= scale;
844  break;
845  case 1:
846  _mat[1] *= scale;
847  _mat[4] *= scale;
848  _mat[7] *= scale;
849  break;
850  case 2:
851  _mat[2] *= scale;
852  _mat[5] *= scale;
853  _mat[8] *= scale;
854  break;
855  }
856 #endif
857  }

References _mat, and size.

Referenced by moab::box_from_axes(), and moab::OrientedBox::order_axes_by_length().

◆ copycol()

template<typename Vector >
void moab::Matrix3::copycol ( int  index,
Vector &  vol 
)
inline

Definition at line 738 of file Matrix3.hpp.

739  {
740 #ifndef MOAB_HAVE_LAPACK
741  _mat.col( index ).swap( vol );
742 #else
743  switch( index )
744  {
745  case 0:
746  _mat[0] = vol[0];
747  _mat[3] = vol[1];
748  _mat[6] = vol[2];
749  break;
750  case 1:
751  _mat[1] = vol[0];
752  _mat[4] = vol[1];
753  _mat[7] = vol[2];
754  break;
755  case 2:
756  _mat[2] = vol[0];
757  _mat[5] = vol[1];
758  _mat[8] = vol[2];
759  break;
760  }
761 #endif
762  }

References _mat.

◆ copyto()

void moab::Matrix3::copyto ( double  v[Matrix3::size])
inline

Definition at line 385 of file Matrix3.hpp.

386  {
387 #ifndef MOAB_HAVE_LAPACK
388  std::copy( _mat.data(), _mat.data() + size, v );
389 #else
390  memcpy( v, _mat, size * sizeof( double ) );
391 #endif
392  }

References _mat, and size.

Referenced by moab::LinearTet::initFcn(), and moab::LinearTri::initFcn().

◆ determinant()

double moab::Matrix3::determinant ( ) const
inline

Definition at line 930 of file Matrix3.hpp.

931  {
932 #ifndef MOAB_HAVE_LAPACK
933  return _mat.determinant();
934 #else
935  return ( _mat[0] * _mat[4] * _mat[8] + _mat[1] * _mat[5] * _mat[6] + _mat[2] * _mat[3] * _mat[7] -
936  _mat[0] * _mat[5] * _mat[7] - _mat[1] * _mat[3] * _mat[8] - _mat[2] * _mat[4] * _mat[6] );
937 #endif
938  }

References _mat.

Referenced by moab::Element::Map::det_ijacobian(), moab::Element::Map::det_jacobian(), moab::EvalSet::evaluate_reverse(), moab::LinearTet::evaluate_reverse(), moab::LinearTri::evaluate_reverse(), moab::Element::Map::ievaluate(), moab::LinearTet::initFcn(), moab::LinearTri::initFcn(), moab::Element::SpectralHex::integrate_scalar_field(), moab::element_utility::Spectral_hex_map< _Matrix >::integrate_scalar_field(), moab::ElemUtil::integrate_trilinear_hex(), moab::LinearHex::integrateFcn(), moab::LinearQuad::integrateFcn(), inverse(), invert(), moab::AffineXform::reflection(), moab::AffineXform::scale(), moab::Element::LinearTet::set_vertices(), moab::Element::LinearTri::set_vertices(), moab::GeomUtil::VolMap::solve_inverse(), moab::ElemUtil::VolMap::solve_inverse(), and test_spectral_hex().

◆ eigen_decomposition()

template<typename Vector >
ErrorCode moab::Matrix3::eigen_decomposition ( Vector &  evals,
Matrix3 &  evecs 
)
inline

Definition at line 589 of file Matrix3.hpp.

590  {
591  const bool bisSymmetric = this->is_symmetric();
592 #ifndef MOAB_HAVE_LAPACK
593  if( bisSymmetric )
594  {
595  Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d > eigensolver( this->_mat );
596  if( eigensolver.info() != Eigen::Success ) return MB_FAILURE;
597  const Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d >::RealVectorType& e3evals = eigensolver.eigenvalues();
598  evals[0] = e3evals( 0 );
599  evals[1] = e3evals( 1 );
600  evals[2] = e3evals( 2 );
601  evecs._mat = eigensolver.eigenvectors(); //.col(1)
602  return MB_SUCCESS;
603  }
604  else
605  {
606  MB_CHK_SET_ERR( MB_FAILURE, "Unsymmetric matrix implementation with Eigen3 is currently not provided." );
607  // Eigen::EigenSolver<Eigen::Matrix3d> eigensolver(this->_mat, true);
608  // if (eigensolver.info() != Eigen::Success)
609  // return MB_FAILURE;
610  // const Eigen::EigenSolver<Eigen::Matrix3d>::EigenvalueType& e3evals =
611  // eigensolver.eigenvalues().real(); evals[0] = e3evals(0); evals[1] = e3evals(1);
612  // evals[2] = e3evals(2); evecs._mat = eigensolver.eigenvectors().real(); //.col(1)
613  // return MB_SUCCESS;
614  }
615 #else
616  int info;
617  /* Solve eigenproblem */
618  double devreal[3], drevecs[9];
619  if( !bisSymmetric )
620  {
621  double devimag[3], dlevecs[9], dwork[102];
622  char dgeev_opts[2] = { 'N', 'V' };
623  int N = 3, LWORK = 102, NL = 1, NR = N;
624  std::vector< double > devmat;
625  devmat.assign( _mat, _mat + size );
626  MOAB_dgeev( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, devimag, dlevecs, &NL, drevecs,
627  &NR, dwork, &LWORK, &info );
628  // The result eigenvalues are ordered as high-->low
629  evals[0] = devreal[2];
630  evals[1] = devreal[1];
631  evals[2] = devreal[0];
632  evecs._mat[0] = drevecs[6];
633  evecs._mat[1] = drevecs[3];
634  evecs._mat[2] = drevecs[0];
635  evecs._mat[3] = drevecs[7];
636  evecs._mat[4] = drevecs[4];
637  evecs._mat[5] = drevecs[1];
638  evecs._mat[6] = drevecs[8];
639  evecs._mat[7] = drevecs[5];
640  evecs._mat[8] = drevecs[2];
641  std::cout << "DGEEV: Optimal work vector: dsize = " << dwork[0] << ".\n";
642  }
643  else
644  {
645  char dgeev_opts[2] = { 'V', 'L' };
646  const bool find_optimal = false;
647  std::vector< int > iwork( 18 );
648  std::vector< double > devmat( 9, 0.0 );
649  std::vector< double > dwork( 38 );
650  int N = 3, lwork = 38, liwork = 18;
651  devmat[0] = _mat[0];
652  devmat[1] = _mat[1];
653  devmat[2] = _mat[2];
654  devmat[4] = _mat[4];
655  devmat[5] = _mat[5];
656  devmat[8] = _mat[8];
657  if( find_optimal )
658  {
659  int _lwork = -1;
660  int _liwork = -1;
661  double query_work_size = 0;
662  int query_iwork_size = 0;
663  // Make an empty call to find the optimal work vector size
664  MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, NULL, &N, NULL, &query_work_size, &_lwork,
665  &query_iwork_size, &_liwork, &info );
666  lwork = (int)query_work_size;
667  dwork.resize( lwork );
668  liwork = query_iwork_size;
669  iwork.resize( liwork );
670  std::cout << "DSYEVD: Optimal work vector: dsize = " << lwork << ", and isize = " << liwork << ".\n";
671  }
672 
673  MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, &dwork[0], &lwork, &iwork[0],
674  &liwork, &info );
675  for( int i = 0; i < 9; ++i )
676  drevecs[i] = devmat[i];
677  // The result eigenvalues are ordered as low-->high, but vectors are in rows of A.
678  evals[0] = devreal[0];
679  evals[1] = devreal[1];
680  evals[2] = devreal[2];
681  evecs._mat[0] = drevecs[0];
682  evecs._mat[3] = drevecs[1];
683  evecs._mat[6] = drevecs[2];
684  evecs._mat[1] = drevecs[3];
685  evecs._mat[4] = drevecs[4];
686  evecs._mat[7] = drevecs[5];
687  evecs._mat[2] = drevecs[6];
688  evecs._mat[5] = drevecs[7];
689  evecs._mat[8] = drevecs[8];
690  }
691 
692  if( !info )
693  {
694  return MB_SUCCESS;
695  }
696  else
697  {
698  std::cout << "Failure in LAPACK_" << ( bisSymmetric ? "DSYEVD" : "DGEEV" )
699  << " call for eigen decomposition.\n";
700  std::cout << "Failed with error = " << info << ".\n";
701  return MB_FAILURE;
702  }
703 #endif
704  }

References _mat, is_symmetric(), MB_CHK_SET_ERR, MB_SUCCESS, and size.

Referenced by moab::OrientedBox::compute_from_covariance_data(), and moab::OrientedBox::compute_from_vertices().

◆ inverse()

Matrix3 moab::Matrix3::inverse ( ) const
inline

Definition at line 940 of file Matrix3.hpp.

941  {
942 #ifndef MOAB_HAVE_LAPACK
943  return Matrix3( _mat.inverse() );
944 #else
945  // return Matrix::compute_inverse( *this, this->determinant() );
946  Matrix3 m( 0.0 );
947  const double d_determinant = 1.0 / this->determinant();
948  m._mat[0] = d_determinant * ( _mat[4] * _mat[8] - _mat[5] * _mat[7] );
949  m._mat[1] = d_determinant * ( _mat[2] * _mat[7] - _mat[8] * _mat[1] );
950  m._mat[2] = d_determinant * ( _mat[1] * _mat[5] - _mat[4] * _mat[2] );
951  m._mat[3] = d_determinant * ( _mat[5] * _mat[6] - _mat[8] * _mat[3] );
952  m._mat[4] = d_determinant * ( _mat[0] * _mat[8] - _mat[6] * _mat[2] );
953  m._mat[5] = d_determinant * ( _mat[2] * _mat[3] - _mat[5] * _mat[0] );
954  m._mat[6] = d_determinant * ( _mat[3] * _mat[7] - _mat[6] * _mat[4] );
955  m._mat[7] = d_determinant * ( _mat[1] * _mat[6] - _mat[7] * _mat[0] );
956  m._mat[8] = d_determinant * ( _mat[0] * _mat[4] - _mat[3] * _mat[1] );
957  return m;
958 #endif
959  }

References _mat, determinant(), and Matrix3().

Referenced by moab::Element::Map::det_ijacobian(), moab::EvalSet::evaluate_reverse(), moab::LinearTet::evaluate_reverse(), moab::LinearTri::evaluate_reverse(), moab::Element::Map::ievaluate(), moab::Element::Map::ijacobian(), moab::LinearTet::initFcn(), moab::LinearTri::initFcn(), moab::AffineXform::inverse(), moab::Element::LinearTet::set_vertices(), moab::Element::LinearTri::set_vertices(), moab::element_utility::Linear_hex_map::solve_inverse(), moab::GeomUtil::VolMap::solve_inverse(), and moab::ElemUtil::VolMap::solve_inverse().

◆ invert()

bool moab::Matrix3::invert ( )
inline

Definition at line 961 of file Matrix3.hpp.

962  {
963  bool invertible = false;
964  double d_determinant;
965 #ifndef MOAB_HAVE_LAPACK
966  Eigen::Matrix3d invMat;
967  _mat.computeInverseAndDetWithCheck( invMat, d_determinant, invertible );
968  if( !Util::is_finite( d_determinant ) ) return false;
969  _mat = invMat;
970  return invertible;
971 #else
972  d_determinant = this->determinant();
973  if( d_determinant > 1e-13 ) invertible = true;
974  d_determinant = 1.0 / d_determinant; // invert the determinant
975  std::vector< double > _m;
976  _m.assign( _mat, _mat + size );
977  _mat[0] = d_determinant * ( _m[4] * _m[8] - _m[5] * _m[7] );
978  _mat[1] = d_determinant * ( _m[2] * _m[7] - _m[8] * _m[1] );
979  _mat[2] = d_determinant * ( _m[1] * _m[5] - _m[4] * _m[2] );
980  _mat[3] = d_determinant * ( _m[5] * _m[6] - _m[8] * _m[3] );
981  _mat[4] = d_determinant * ( _m[0] * _m[8] - _m[6] * _m[2] );
982  _mat[5] = d_determinant * ( _m[2] * _m[3] - _m[5] * _m[0] );
983  _mat[6] = d_determinant * ( _m[3] * _m[7] - _m[6] * _m[4] );
984  _mat[7] = d_determinant * ( _m[1] * _m[6] - _m[7] * _m[0] );
985  _mat[8] = d_determinant * ( _m[0] * _m[4] - _m[3] * _m[1] );
986 #endif
987  return invertible;
988  }

References _mat, determinant(), moab::Util::is_finite(), and size.

◆ is_positive_definite()

bool moab::Matrix3::is_positive_definite ( )
inline

Definition at line 569 of file Matrix3.hpp.

570  {
571 #ifndef MOAB_HAVE_LAPACK
572  double subdet6 = _mat( 1 ) * _mat( 5 ) - _mat( 2 ) * _mat( 4 );
573  double subdet7 = _mat( 2 ) * _mat( 3 ) - _mat( 0 ) * _mat( 5 );
574  double subdet8 = _mat( 0 ) * _mat( 4 ) - _mat( 1 ) * _mat( 3 );
575  // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
576  const double det = _mat( 6 ) * subdet6 + _mat( 7 ) * subdet7 + _mat( 8 ) * subdet8;
577  return _mat( 0 ) > 0 && subdet8 > 0 && det > 0;
578 #else
579  double subdet6 = _mat[1] * _mat[5] - _mat[2] * _mat[4];
580  double subdet7 = _mat[2] * _mat[3] - _mat[0] * _mat[5];
581  double subdet8 = _mat[0] * _mat[4] - _mat[1] * _mat[3];
582  // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
583  const double det = _mat[6] * subdet6 + _mat[7] * subdet7 + _mat[8] * subdet8;
584  return _mat[0] > 0 && subdet8 > 0 && det > 0;
585 #endif
586  }

References _mat.

◆ is_symmetric()

bool moab::Matrix3::is_symmetric ( )
inline

Definition at line 553 of file Matrix3.hpp.

554  {
555  const double EPS = 1e-13;
556 #ifndef MOAB_HAVE_LAPACK
557  if( ( fabs( _mat( 1 ) - _mat( 3 ) ) < EPS ) && ( fabs( _mat( 2 ) - _mat( 6 ) ) < EPS ) &&
558  ( fabs( _mat( 5 ) - _mat( 7 ) ) < EPS ) )
559  return true;
560 #else
561  if( ( fabs( _mat[1] - _mat[3] ) < EPS ) && ( fabs( _mat[2] - _mat[6] ) < EPS ) &&
562  ( fabs( _mat[5] - _mat[7] ) < EPS ) )
563  return true;
564 #endif
565  else
566  return false;
567  }

References _mat, and EPS.

Referenced by eigen_decomposition().

◆ operator()() [1/4]

double& moab::Matrix3::operator() ( unsigned  i)
inline

Definition at line 450 of file Matrix3.hpp.

451  {
452 #ifndef MOAB_HAVE_LAPACK
453  return _mat( i );
454 #else
455  return _mat[i];
456 #endif
457  }

References _mat.

◆ operator()() [2/4]

double moab::Matrix3::operator() ( unsigned  i) const
inline

Definition at line 459 of file Matrix3.hpp.

460  {
461 #ifndef MOAB_HAVE_LAPACK
462  return _mat( i );
463 #else
464  return _mat[i];
465 #endif
466  }

References _mat.

◆ operator()() [3/4]

double& moab::Matrix3::operator() ( unsigned  r,
unsigned  c 
)
inline

Definition at line 432 of file Matrix3.hpp.

433  {
434 #ifndef MOAB_HAVE_LAPACK
435  return _mat( r, c );
436 #else
437  return _mat[r * 3 + c];
438 #endif
439  }

References _mat.

◆ operator()() [4/4]

double moab::Matrix3::operator() ( unsigned  r,
unsigned  c 
) const
inline

Definition at line 441 of file Matrix3.hpp.

442  {
443 #ifndef MOAB_HAVE_LAPACK
444  return _mat( r, c );
445 #else
446  return _mat[r * 3 + c];
447 #endif
448  }

References _mat.

◆ operator*=() [1/2]

Matrix3& moab::Matrix3::operator*= ( const Matrix3 &  m)
inline

Definition at line 531 of file Matrix3.hpp.

532  {
533 #ifndef MOAB_HAVE_LAPACK
534  _mat *= m._mat;
535 #else
536  // Uncomment below if you want point-wise multiplication instead (.*)
537  // for (int i=0; i < Matrix3::size; ++i) _mat[i] *= m._mat[i];
538  std::vector< double > dmat;
539  dmat.assign( _mat, _mat + size );
540  _mat[0] = dmat[0] * m._mat[0] + dmat[1] * m._mat[3] + dmat[2] * m._mat[6];
541  _mat[1] = dmat[0] * m._mat[1] + dmat[1] * m._mat[4] + dmat[2] * m._mat[7];
542  _mat[2] = dmat[0] * m._mat[2] + dmat[1] * m._mat[5] + dmat[2] * m._mat[8];
543  _mat[3] = dmat[3] * m._mat[0] + dmat[4] * m._mat[3] + dmat[5] * m._mat[6];
544  _mat[4] = dmat[3] * m._mat[1] + dmat[4] * m._mat[4] + dmat[5] * m._mat[7];
545  _mat[5] = dmat[3] * m._mat[2] + dmat[4] * m._mat[5] + dmat[5] * m._mat[8];
546  _mat[6] = dmat[6] * m._mat[0] + dmat[7] * m._mat[3] + dmat[8] * m._mat[6];
547  _mat[7] = dmat[6] * m._mat[1] + dmat[7] * m._mat[4] + dmat[8] * m._mat[7];
548  _mat[8] = dmat[6] * m._mat[2] + dmat[7] * m._mat[5] + dmat[8] * m._mat[8];
549 #endif
550  return *this;
551  }

References _mat, and size.

◆ operator*=() [2/2]

Matrix3& moab::Matrix3::operator*= ( double  s)
inline

Definition at line 509 of file Matrix3.hpp.

510  {
511 #ifndef MOAB_HAVE_LAPACK
512  _mat *= s;
513 #else
514  for( int i = 0; i < Matrix3::size; ++i )
515  _mat[i] *= s;
516 #endif
517  return *this;
518  }

References _mat, and size.

◆ operator+=()

Matrix3& moab::Matrix3::operator+= ( const Matrix3 &  m)
inline

Definition at line 487 of file Matrix3.hpp.

488  {
489 #ifndef MOAB_HAVE_LAPACK
490  _mat += m._mat;
491 #else
492  for( int i = 0; i < Matrix3::size; ++i )
493  _mat[i] += m._mat[i];
494 #endif
495  return *this;
496  }

References _mat, and size.

◆ operator-=()

Matrix3& moab::Matrix3::operator-= ( const Matrix3 &  m)
inline

Definition at line 498 of file Matrix3.hpp.

499  {
500 #ifndef MOAB_HAVE_LAPACK
501  _mat -= m._mat;
502 #else
503  for( int i = 0; i < Matrix3::size; ++i )
504  _mat[i] -= m._mat[i];
505 #endif
506  return *this;
507  }

References _mat, and size.

◆ operator/=()

Matrix3& moab::Matrix3::operator/= ( double  s)
inline

Definition at line 520 of file Matrix3.hpp.

521  {
522 #ifndef MOAB_HAVE_LAPACK
523  _mat /= s;
524 #else
525  for( int i = 0; i < Matrix3::size; ++i )
526  _mat[i] /= s;
527 #endif
528  return *this;
529  }

References _mat, and size.

◆ operator=() [1/2]

Matrix3& moab::Matrix3::operator= ( const double  v[size])
inline

Definition at line 404 of file Matrix3.hpp.

405  {
406 #ifndef MOAB_HAVE_LAPACK
407  _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
408 #else
409  memcpy( _mat, v, size * sizeof( double ) );
410 #endif
411  return *this;
412  }

References _mat, and size.

◆ operator=() [2/2]

Matrix3& moab::Matrix3::operator= ( const Matrix3 &  m)
inline

Definition at line 394 of file Matrix3.hpp.

395  {
396 #ifndef MOAB_HAVE_LAPACK
397  _mat = m._mat;
398 #else
399  memcpy( _mat, m._mat, size * sizeof( double ) );
400 #endif
401  return *this;
402  }

References _mat, and size.

◆ operator[]() [1/2]

double* moab::Matrix3::operator[] ( unsigned  i)
inline

Definition at line 414 of file Matrix3.hpp.

415  {
416 #ifndef MOAB_HAVE_LAPACK
417  return _mat.row( i ).data();
418 #else
419  return &_mat[i * 3]; // Row Major
420 #endif
421  }

References _mat.

◆ operator[]() [2/2]

const double* moab::Matrix3::operator[] ( unsigned  i) const
inline

Definition at line 423 of file Matrix3.hpp.

424  {
425 #ifndef MOAB_HAVE_LAPACK
426  return _mat.row( i ).data();
427 #else
428  return &_mat[i * 3];
429 #endif
430  }

References _mat.

◆ print()

void moab::Matrix3::print ( std::ostream &  s) const
inline

Definition at line 1007 of file Matrix3.hpp.

1008  {
1009 #ifndef MOAB_HAVE_LAPACK
1010  s << "| " << _mat( 0 ) << " " << _mat( 1 ) << " " << _mat( 2 ) << " | " << _mat( 3 ) << " " << _mat( 4 ) << " "
1011  << _mat( 5 ) << " | " << _mat( 6 ) << " " << _mat( 7 ) << " " << _mat( 8 ) << " |";
1012 #else
1013  s << "| " << _mat[0] << " " << _mat[1] << " " << _mat[2] << " | " << _mat[3] << " " << _mat[4] << " " << _mat[5]
1014  << " | " << _mat[6] << " " << _mat[7] << " " << _mat[8] << " |";
1015 #endif
1016  }

References _mat.

◆ row()

CartVect moab::Matrix3::row ( int  index) const
inline

Definition at line 906 of file Matrix3.hpp.

907  {
908  assert( index < Matrix3::size );
909 #ifndef MOAB_HAVE_LAPACK
910  Eigen::Vector3d mvec = _mat.row( index );
911  return CartVect( mvec[0], mvec[1], mvec[2] );
912 #else
913  switch( index )
914  {
915  case 0:
916  return CartVect( _mat[0], _mat[1], _mat[2] );
917  case 1:
918  return CartVect( _mat[3], _mat[4], _mat[5] );
919  case 2:
920  return CartVect( _mat[6], _mat[7], _mat[8] );
921  }
922  return CartVect( 0.0 );
923 #endif
924  }

References _mat, and size.

◆ rowscale()

void moab::Matrix3::rowscale ( int  index,
double  scale 
)
inline

Definition at line 859 of file Matrix3.hpp.

860  {
861  assert( index < Matrix3::size );
862 #ifndef MOAB_HAVE_LAPACK
863  _mat.row( index ) *= scale;
864 #else
865  switch( index )
866  {
867  case 0:
868  _mat[0] *= scale;
869  _mat[1] *= scale;
870  _mat[2] *= scale;
871  break;
872  case 1:
873  _mat[3] *= scale;
874  _mat[4] *= scale;
875  _mat[5] *= scale;
876  break;
877  case 2:
878  _mat[6] *= scale;
879  _mat[7] *= scale;
880  _mat[8] *= scale;
881  break;
882  }
883 #endif
884  }

References _mat, and size.

◆ subdet()

double moab::Matrix3::subdet ( int  r,
int  c 
) const
inline

Definition at line 992 of file Matrix3.hpp.

993  {
994  assert( r >= 0 && c >= 0 );
995  if( r < 0 || c < 0 ) return DBL_MAX;
996 #ifndef MOAB_HAVE_LAPACK
997  const int r1 = ( r + 1 ) % 3, r2 = ( r + 2 ) % 3;
998  const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
999  return _mat( r1, c1 ) * _mat( r2, c2 ) - _mat( r1, c2 ) * _mat( r2, c1 );
1000 #else
1001  const int r1 = Matrix3::size * ( ( r + 1 ) % 3 ), r2 = Matrix3::size * ( ( r + 2 ) % 3 );
1002  const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
1003  return _mat[r1 + c1] * _mat[r2 + c2] - _mat[r1 + c2] * _mat[r2 + c1];
1004 #endif
1005  }

References _mat, and size.

◆ swapcol()

void moab::Matrix3::swapcol ( int  srcindex,
int  destindex 
)
inline

Definition at line 764 of file Matrix3.hpp.

765  {
766  assert( srcindex < Matrix3::size );
767  assert( destindex < Matrix3::size );
768 #ifndef MOAB_HAVE_LAPACK
769  _mat.col( srcindex ).swap( _mat.col( destindex ) );
770 #else
771  CartVect svol = this->vcol< CartVect >( srcindex );
772  CartVect dvol = this->vcol< CartVect >( destindex );
773  switch( srcindex )
774  {
775  case 0:
776  _mat[0] = dvol[0];
777  _mat[3] = dvol[1];
778  _mat[6] = dvol[2];
779  break;
780  case 1:
781  _mat[1] = dvol[0];
782  _mat[4] = dvol[1];
783  _mat[7] = dvol[2];
784  break;
785  case 2:
786  _mat[2] = dvol[0];
787  _mat[5] = dvol[1];
788  _mat[8] = dvol[2];
789  break;
790  }
791  switch( destindex )
792  {
793  case 0:
794  _mat[0] = svol[0];
795  _mat[3] = svol[1];
796  _mat[6] = svol[2];
797  break;
798  case 1:
799  _mat[1] = svol[0];
800  _mat[4] = svol[1];
801  _mat[7] = svol[2];
802  break;
803  case 2:
804  _mat[2] = svol[0];
805  _mat[5] = svol[1];
806  _mat[8] = svol[2];
807  break;
808  }
809 #endif
810  }

References _mat, and size.

Referenced by moab::box_from_axes(), and moab::OrientedBox::order_axes_by_length().

◆ transpose()

Matrix3 moab::Matrix3::transpose ( ) const
inline

Definition at line 721 of file Matrix3.hpp.

722  {
723 #ifndef MOAB_HAVE_LAPACK
724  return Matrix3( _mat.transpose() );
725 #else
726  Matrix3 mtmp( *this );
727  mtmp._mat[1] = _mat[3];
728  mtmp._mat[3] = _mat[1];
729  mtmp._mat[2] = _mat[6];
730  mtmp._mat[6] = _mat[2];
731  mtmp._mat[5] = _mat[7];
732  mtmp._mat[7] = _mat[5];
733  return mtmp;
734 #endif
735  }

References _mat, and Matrix3().

Referenced by moab::OrientedBox::intersect_ray().

◆ transpose_inplace()

void moab::Matrix3::transpose_inplace ( )
inline

Definition at line 706 of file Matrix3.hpp.

707  {
708 #ifndef MOAB_HAVE_LAPACK
709  _mat.transposeInPlace();
710 #else
711  Matrix3 mtmp( *this );
712  _mat[1] = mtmp._mat[3];
713  _mat[3] = mtmp._mat[1];
714  _mat[2] = mtmp._mat[6];
715  _mat[6] = mtmp._mat[2];
716  _mat[5] = mtmp._mat[7];
717  _mat[7] = mtmp._mat[5];
718 #endif
719  }

References _mat.

◆ vcol()

template<typename Vector >
Vector moab::Matrix3::vcol ( int  index) const
inline

Definition at line 813 of file Matrix3.hpp.

814  {
815  assert( index < Matrix3::size );
816 #ifndef MOAB_HAVE_LAPACK
817  return _mat.col( index );
818 #else
819  switch( index )
820  {
821  case 0:
822  return Vector( _mat[0], _mat[3], _mat[6] );
823  case 1:
824  return Vector( _mat[1], _mat[4], _mat[7] );
825  case 2:
826  return Vector( _mat[2], _mat[5], _mat[8] );
827  }
828  return Vector( 0.0 );
829 #endif
830  }

References _mat, and size.

Friends And Related Function Documentation

◆ operator*

Matrix3 operator* ( const Matrix3 &  a,
const Matrix3 &  b 
)
friend

Definition at line 1051 of file Matrix3.hpp.

1052 {
1053 #ifndef MOAB_HAVE_LAPACK
1054  return Matrix3( a._mat * b._mat );
1055 #else
1056  return Matrix::mmult3( a, b );
1057 #endif
1058 }

◆ operator+

Matrix3 operator+ ( const Matrix3 &  a,
const Matrix3 &  b 
)
friend

Definition at line 1027 of file Matrix3.hpp.

1028 {
1029 #ifndef MOAB_HAVE_LAPACK
1030  return Matrix3( a._mat + b._mat );
1031 #else
1032  Matrix3 s( a );
1033  for( int i = 0; i < Matrix3::size; ++i )
1034  s( i ) += b._mat[i];
1035  return s;
1036 #endif
1037 }

◆ operator-

Matrix3 operator- ( const Matrix3 &  a,
const Matrix3 &  b 
)
friend

Definition at line 1039 of file Matrix3.hpp.

1040 {
1041 #ifndef MOAB_HAVE_LAPACK
1042  return Matrix3( a._mat - b._mat );
1043 #else
1044  Matrix3 s( a );
1045  for( int i = 0; i < Matrix3::size; ++i )
1046  s( i ) -= b._mat[i];
1047  return s;
1048 #endif
1049 }

Member Data Documentation

◆ _mat

Eigen::Matrix3d moab::Matrix3::_mat
private

Definition at line 224 of file Matrix3.hpp.

Referenced by array(), col(), colscale(), copycol(), copyto(), determinant(), eigen_decomposition(), inverse(), invert(), is_positive_definite(), is_symmetric(), Matrix3(), operator()(), operator*=(), operator+=(), operator-=(), operator/=(), operator=(), operator[](), print(), row(), rowscale(), subdet(), swapcol(), transpose(), transpose_inplace(), and vcol().

◆ size

const int moab::Matrix3::size = 9
static

Definition at line 220 of file Matrix3.hpp.

Referenced by col(), colscale(), copyto(), eigen_decomposition(), moab::LinearTet::initFcn(), moab::LinearTri::initFcn(), invert(), Matrix3(), operator*=(), operator+=(), operator-=(), operator/=(), operator=(), row(), rowscale(), subdet(), swapcol(), and vcol().

The documentation for this class was generated from the following file: