24 #ifndef MOAB_MATRIX3_HPP
25 #define MOAB_MATRIX3_HPP
39 #ifndef MOAB_HAVE_LAPACK
41 #ifndef MOAB_HAVE_EIGEN3
42 #error Need either Eigen3 or BLAS/LAPACK libraries
47 #pragma GCC diagnostic push
49 #pragma GCC diagnostic ignored "-Wshadow"
52 #define EIGEN_DEFAULT_TO_ROW_MAJOR
53 #define EIGEN_INITIALIZE_MATRICES_BY_ZERO
55 #include "Eigen/Dense"
59 #pragma GCC diagnostic pop
64 #if defined( MOAB_FC_FUNC_ )
65 #define MOAB_FC_WRAPPER MOAB_FC_FUNC_
66 #elif defined( MOAB_FC_FUNC )
67 #define MOAB_FC_WRAPPER MOAB_FC_FUNC
69 #define MOAB_FC_WRAPPER( name, NAME ) name##_
79 #define MOAB_dsyevd MOAB_FC_FUNC( dsyevd, DSYEVD )
80 #define MOAB_dsyevr MOAB_FC_FUNC( dsyevr, DSYEVR )
81 #define MOAB_dgeev MOAB_FC_FUNC( dgeev, DGEEV )
82 #define MOAB_dgetrf MOAB_FC_FUNC( dgetrf, DGETRF )
83 #define MOAB_dgetri MOAB_FC_FUNC( dgetri, DGETRI )
87 #define MOAB_dsyevd MOAB_FC_WRAPPER( dsyevd, DSYEVD )
88 #define MOAB_dsyevr MOAB_FC_WRAPPER( dsyevr, DSYEVR )
89 #define MOAB_dgeev MOAB_FC_WRAPPER( dgeev, DGEEV )
90 #define MOAB_dgetrf MOAB_FC_WRAPPER( dgetrf, DGETRF )
91 #define MOAB_dgetri MOAB_FC_WRAPPER( dgetri, DGETRI )
100 void MOAB_dsyevd(
char* jobz,
116 void MOAB_dsyevr(
char* jobz,
140 void MOAB_dgeev(
char* jobvl,
157 void MOAB_dgetrf(
int* M,
int* N,
double* A,
int* lda,
int* IPIV,
int* INFO );
161 void MOAB_dgetri(
int* N,
double* A,
int* lda,
int* IPIV,
double* WORK,
int* lwork,
int* INFO );
165 #define MOAB_DMEMZERO( a, b ) memset( a, 0, ( b ) * sizeof( double ) )
175 template <
typename Matrix >
176 inline Matrix
mmult3(
const Matrix& a,
const Matrix& b )
178 return Matrix( a( 0, 0 ) * b( 0, 0 ) + a( 0, 1 ) * b( 1, 0 ) + a( 0, 2 ) * b( 2, 0 ),
179 a( 0, 0 ) * b( 0, 1 ) + a( 0, 1 ) * b( 1, 1 ) + a( 0, 2 ) * b( 2, 1 ),
180 a( 0, 0 ) * b( 0, 2 ) + a( 0, 1 ) * b( 1, 2 ) + a( 0, 2 ) * b( 2, 2 ),
181 a( 1, 0 ) * b( 0, 0 ) + a( 1, 1 ) * b( 1, 0 ) + a( 1, 2 ) * b( 2, 0 ),
182 a( 1, 0 ) * b( 0, 1 ) + a( 1, 1 ) * b( 1, 1 ) + a( 1, 2 ) * b( 2, 1 ),
183 a( 1, 0 ) * b( 0, 2 ) + a( 1, 1 ) * b( 1, 2 ) + a( 1, 2 ) * b( 2, 2 ),
184 a( 2, 0 ) * b( 0, 0 ) + a( 2, 1 ) * b( 1, 0 ) + a( 2, 2 ) * b( 2, 0 ),
185 a( 2, 0 ) * b( 0, 1 ) + a( 2, 1 ) * b( 1, 1 ) + a( 2, 2 ) * b( 2, 1 ),
186 a( 2, 0 ) * b( 0, 2 ) + a( 2, 1 ) * b( 1, 2 ) + a( 2, 2 ) * b( 2, 2 ) );
189 template <
typename Matrix >
190 inline const Matrix
inverse(
const Matrix& d )
192 const double det = 1.0 / determinant3( d );
196 template <
typename Vector,
typename Matrix >
199 return Vector( v[0] * m( 0, 0 ) + v[1] * m( 1, 0 ) + v[2] * m( 2, 0 ),
200 v[0] * m( 0, 1 ) + v[1] * m( 1, 1 ) + v[2] * m( 2, 1 ),
201 v[0] * m( 0, 2 ) + v[1] * m( 1, 2 ) + v[2] * m( 2, 2 ) );
204 template <
typename Vector,
typename Matrix >
208 res[0] = v[0] * m( 0, 0 ) + v[1] * m( 0, 1 ) + v[2] * m( 0, 2 );
209 res[1] = v[0] * m( 1, 0 ) + v[1] * m( 1, 1 ) + v[2] * m( 1, 2 );
210 res[2] = v[0] * m( 2, 0 ) + v[1] * m( 2, 1 ) + v[2] * m( 2, 2 );
223 #ifndef MOAB_HAVE_LAPACK
233 #ifndef MOAB_HAVE_LAPACK
240 #ifndef MOAB_HAVE_LAPACK
248 #ifndef MOAB_HAVE_LAPACK
249 _mat << diagonal, 0.0, 0.0, 0.0, diagonal, 0.0, 0.0, 0.0, diagonal;
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];
262 _mat[0] = diagonal[0];
263 _mat[4] = diagonal[1];
264 _mat[8] = diagonal[2];
273 inline Matrix3(
const std::vector< double >& diagonal )
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];
279 _mat[0] = diagonal[0];
280 _mat[4] = diagonal[1];
281 _mat[8] = diagonal[2];
295 #ifndef MOAB_HAVE_LAPACK
296 _mat << v00, v01, v02, v10, v11, v12, v20, v21, v22;
314 #ifndef MOAB_HAVE_LAPACK
322 template <
typename Vector >
323 inline Matrix3(
const Vector& row0,
const Vector& row1,
const Vector& row2,
const bool isRow )
325 #ifndef MOAB_HAVE_LAPACK
328 _mat << row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], row2[0], row2[1], row2[2];
332 _mat << row0[0], row1[0], row2[0], row0[1], row1[1], row2[1], row0[2], row1[2], row2[2];
365 #define DEPRECATED __attribute__( ( deprecated ) )
367 #pragma message( "WARNING: You need to implement DEPRECATED for this compiler" )
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];
381 memcpy(
_mat, v,
size *
sizeof(
double ) );
387 #ifndef MOAB_HAVE_LAPACK
390 memcpy( v,
_mat,
size *
sizeof(
double ) );
396 #ifndef MOAB_HAVE_LAPACK
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];
409 memcpy(
_mat, v,
size *
sizeof(
double ) );
416 #ifndef MOAB_HAVE_LAPACK
417 return _mat.row( i ).data();
425 #ifndef MOAB_HAVE_LAPACK
426 return _mat.row( i ).data();
434 #ifndef MOAB_HAVE_LAPACK
437 return _mat[r * 3 + c];
443 #ifndef MOAB_HAVE_LAPACK
446 return _mat[r * 3 + c];
452 #ifndef MOAB_HAVE_LAPACK
461 #ifndef MOAB_HAVE_LAPACK
471 #ifndef MOAB_HAVE_LAPACK
480 #ifndef MOAB_HAVE_LAPACK
489 #ifndef MOAB_HAVE_LAPACK
500 #ifndef MOAB_HAVE_LAPACK
511 #ifndef MOAB_HAVE_LAPACK
522 #ifndef MOAB_HAVE_LAPACK
533 #ifndef MOAB_HAVE_LAPACK
538 std::vector< double > dmat;
555 const double EPS = 1e-13;
556 #ifndef MOAB_HAVE_LAPACK
571 #ifndef MOAB_HAVE_LAPACK
576 const double det =
_mat( 6 ) * subdet6 +
_mat( 7 ) * subdet7 +
_mat( 8 ) * subdet8;
577 return _mat( 0 ) > 0 && subdet8 > 0 && det > 0;
583 const double det =
_mat[6] * subdet6 +
_mat[7] * subdet7 +
_mat[8] * subdet8;
584 return _mat[0] > 0 && subdet8 > 0 && det > 0;
588 template <
typename Vector >
592 #ifndef MOAB_HAVE_LAPACK
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();
606 MB_CHK_SET_ERR( MB_FAILURE,
"Unsymmetric matrix implementation with Eigen3 is currently not provided." );
618 double devreal[3], drevecs[9];
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;
626 MOAB_dgeev( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, devimag, dlevecs, &NL, drevecs,
627 &NR, dwork, &LWORK, &info );
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";
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;
661 double query_work_size = 0;
662 int query_iwork_size = 0;
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";
673 MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, &dwork[0], &lwork, &iwork[0],
675 for(
int i = 0; i < 9; ++i )
676 drevecs[i] = devmat[i];
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];
698 std::cout <<
"Failure in LAPACK_" << ( bisSymmetric ?
"DSYEVD" :
"DGEEV" )
699 <<
" call for eigen decomposition.\n";
700 std::cout <<
"Failed with error = " << info <<
".\n";
708 #ifndef MOAB_HAVE_LAPACK
709 _mat.transposeInPlace();
723 #ifndef MOAB_HAVE_LAPACK
737 template <
typename Vector >
740 #ifndef MOAB_HAVE_LAPACK
741 _mat.col( index ).swap( vol );
764 inline void swapcol(
int srcindex,
int destindex )
768 #ifndef MOAB_HAVE_LAPACK
769 _mat.col( srcindex ).swap(
_mat.col( destindex ) );
771 CartVect svol = this->vcol< CartVect >( srcindex );
772 CartVect dvol = this->vcol< CartVect >( destindex );
812 template <
typename Vector >
813 inline Vector
vcol(
int index )
const
816 #ifndef MOAB_HAVE_LAPACK
817 return _mat.col( index );
828 return Vector( 0.0 );
835 #ifndef MOAB_HAVE_LAPACK
836 _mat.col( index ) *= scale;
862 #ifndef MOAB_HAVE_LAPACK
863 _mat.row( index ) *= scale;
889 #ifndef MOAB_HAVE_LAPACK
890 Eigen::Vector3d mvec =
_mat.col( index );
891 return CartVect( mvec[0], mvec[1], mvec[2] );
909 #ifndef MOAB_HAVE_LAPACK
910 Eigen::Vector3d mvec =
_mat.row( index );
911 return CartVect( mvec[0], mvec[1], mvec[2] );
932 #ifndef MOAB_HAVE_LAPACK
933 return _mat.determinant();
942 #ifndef MOAB_HAVE_LAPACK
947 const double d_determinant = 1.0 / this->
determinant();
963 bool invertible =
false;
964 double d_determinant;
965 #ifndef MOAB_HAVE_LAPACK
966 Eigen::Matrix3d invMat;
967 _mat.computeInverseAndDetWithCheck( invMat, d_determinant, invertible );
973 if( d_determinant > 1e-13 ) invertible =
true;
974 d_determinant = 1.0 / d_determinant;
975 std::vector< double > _m;
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] );
992 inline double subdet(
int r,
int c )
const
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;
1002 const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
1007 inline void print( std::ostream& s )
const
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 ) <<
" |";
1014 <<
" | " <<
_mat[6] <<
" " <<
_mat[7] <<
" " <<
_mat[8] <<
" |";
1020 template <
typename Vector >
1023 return Matrix3( u[0] * v[0], u[0] * v[1], u[0] * v[2], u[1] * v[0], u[1] * v[1], u[1] * v[2], u[2] * v[0],
1024 u[2] * v[1], u[2] * v[2] );
1029 #ifndef MOAB_HAVE_LAPACK
1034 s( i ) += b.
_mat[i];
1041 #ifndef MOAB_HAVE_LAPACK
1046 s( i ) -= b.
_mat[i];
1053 #ifndef MOAB_HAVE_LAPACK
1060 template <
typename T >
1066 template <
typename T >
1084 #ifdef MOAB_HAVE_LAPACK
1085 #undef MOAB_DMEMZERO
1088 #ifndef MOAB_MATRIX3_OPERATORLESS
1089 #define MOAB_MATRIX3_OPERATORLESS
1092 return s <<
"| " << m( 0, 0 ) <<
" " << m( 0, 1 ) <<
" " << m( 0, 2 ) <<
" | " << m( 1, 0 ) <<
" " << m( 1, 1 )
1093 <<
" " << m( 1, 2 ) <<
" | " << m( 2, 0 ) <<
" " << m( 2, 1 ) <<
" " << m( 2, 2 ) <<
" |";