MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* 00002 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00003 * storing and accessing finite element mesh data. 00004 * 00005 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00006 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00007 * retains certain rights in 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 00015 /**\file Matrix3.hpp 00016 *\author Jason Kraftcheck (kraftche@cae.wisc.edu) 00017 *\date 2006-07-18 00018 *\date 2012-08-2 Updated by rhl to be more generic. less code that does more! 00019 * TODO: Remove all 'inline' keywords as it is only a suggestion to the compiler 00020 * anyways, and it will ignore it or add it when it thinks its necessary. 00021 *\date 2016-08-03 Updated to use Eigen3 support underneath to improve performance 00022 */ 00023 00024 #ifndef MOAB_MATRIX3_HPP 00025 #define MOAB_MATRIX3_HPP 00026 00027 #include <iostream> 00028 #include <iosfwd> 00029 #include <limits> 00030 #include <cmath> 00031 #include <cassert> 00032 00033 #include "moab/MOABConfig.h" 00034 #include "moab/ErrorHandler.hpp" 00035 #include "moab/Util.hpp" 00036 #include "moab/Types.hpp" 00037 #include "moab/CartVect.hpp" 00038 00039 #ifndef MOAB_HAVE_LAPACK 00040 00041 #ifndef MOAB_HAVE_EIGEN3 00042 #error Need either Eigen3 or BLAS/LAPACK libraries 00043 #endif 00044 00045 #ifdef __GNUC__ 00046 // save diagnostic state 00047 #pragma GCC diagnostic push 00048 // turn off the specific warning. Can also use "-Wshadow" 00049 #pragma GCC diagnostic ignored "-Wshadow" 00050 #endif 00051 00052 #define EIGEN_DEFAULT_TO_ROW_MAJOR 00053 #define EIGEN_INITIALIZE_MATRICES_BY_ZERO 00054 // #define EIGEN_NO_STATIC_ASSERT 00055 #include "Eigen/Dense" 00056 00057 #ifdef __GNUC__ 00058 // turn the warnings back on 00059 #pragma GCC diagnostic pop 00060 #endif 00061 00062 #else 00063 00064 #if defined( MOAB_FC_FUNC_ ) 00065 #define MOAB_FC_WRAPPER MOAB_FC_FUNC_ 00066 #elif defined( MOAB_FC_FUNC ) 00067 #define MOAB_FC_WRAPPER MOAB_FC_FUNC 00068 #else 00069 #define MOAB_FC_WRAPPER( name, NAME ) name##_ 00070 #endif 00071 00072 // We will rely on LAPACK directly 00073 00074 #ifdef WIN32 00075 00076 // Should use second form below for windows but 00077 // needed to do this to make it work. 00078 // TODO: Need to clean this up 00079 #define MOAB_dsyevd MOAB_FC_FUNC( dsyevd, DSYEVD ) 00080 #define MOAB_dsyevr MOAB_FC_FUNC( dsyevr, DSYEVR ) 00081 #define MOAB_dgeev MOAB_FC_FUNC( dgeev, DGEEV ) 00082 #define MOAB_dgetrf MOAB_FC_FUNC( dgetrf, DGETRF ) 00083 #define MOAB_dgetri MOAB_FC_FUNC( dgetri, DGETRI ) 00084 00085 #else 00086 00087 #define MOAB_dsyevd MOAB_FC_WRAPPER( dsyevd, DSYEVD ) 00088 #define MOAB_dsyevr MOAB_FC_WRAPPER( dsyevr, DSYEVR ) 00089 #define MOAB_dgeev MOAB_FC_WRAPPER( dgeev, DGEEV ) 00090 #define MOAB_dgetrf MOAB_FC_WRAPPER( dgetrf, DGETRF ) 00091 #define MOAB_dgetri MOAB_FC_WRAPPER( dgetri, DGETRI ) 00092 00093 #endif 00094 00095 extern "C" { 00096 00097 // Computes all eigenvalues and, optionally, eigenvectors of a 00098 // real symmetric matrix A. If eigenvectors are desired, it uses a 00099 // divide and conquer algorithm. 00100 void MOAB_dsyevd( char* jobz, char* uplo, int* n, double a[], int* lda, double w[], double work[], int* lwork, 00101 int iwork[], int* liwork, int* info ); 00102 00103 // Computes selected eigenvalues and, optionally, eigenvectors 00104 // of a real symmetric matrix A. Eigenvalues and eigenvectors can be 00105 // selected by specifying either a range of values or a range of 00106 // indices for the desired eigenvalues. 00107 void MOAB_dsyevr( char* jobz, char* range, char* uplo, int* n, double* a, int* lda, double* vl, double* vu, int* il, 00108 int* iu, double* abstol, int* m, double* w, double* z, int* ldz, int* isuppz, double* work, 00109 int* lwork, int* iwork, int* liwork, int* info ); 00110 00111 // Computes for an N-by-N real nonsymmetric matrix A, the 00112 // eigenvalues and, optionally, the left and/or right eigenvectors. 00113 void MOAB_dgeev( char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, 00114 double* vr, int* ldvr, double* work, int* lwork, int* info ); 00115 00116 // Computes an LU factorization of a general M-by-N matrix A 00117 // using partial pivoting with row interchanges. 00118 void MOAB_dgetrf( int* M, int* N, double* A, int* lda, int* IPIV, int* INFO ); 00119 00120 // Computes the inverse of a matrix using the LU factorization 00121 // computed by DGETRF. 00122 void MOAB_dgetri( int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO ); 00123 } 00124 00125 #include <cstring> 00126 #define MOAB_DMEMZERO( a, b ) memset( a, 0, b * sizeof( double ) ) 00127 00128 #endif 00129 00130 namespace moab 00131 { 00132 00133 namespace Matrix 00134 { 00135 00136 template < typename Matrix > 00137 inline Matrix mmult3( const Matrix& a, const Matrix& b ) 00138 { 00139 return Matrix( a( 0, 0 ) * b( 0, 0 ) + a( 0, 1 ) * b( 1, 0 ) + a( 0, 2 ) * b( 2, 0 ), 00140 a( 0, 0 ) * b( 0, 1 ) + a( 0, 1 ) * b( 1, 1 ) + a( 0, 2 ) * b( 2, 1 ), 00141 a( 0, 0 ) * b( 0, 2 ) + a( 0, 1 ) * b( 1, 2 ) + a( 0, 2 ) * b( 2, 2 ), 00142 a( 1, 0 ) * b( 0, 0 ) + a( 1, 1 ) * b( 1, 0 ) + a( 1, 2 ) * b( 2, 0 ), 00143 a( 1, 0 ) * b( 0, 1 ) + a( 1, 1 ) * b( 1, 1 ) + a( 1, 2 ) * b( 2, 1 ), 00144 a( 1, 0 ) * b( 0, 2 ) + a( 1, 1 ) * b( 1, 2 ) + a( 1, 2 ) * b( 2, 2 ), 00145 a( 2, 0 ) * b( 0, 0 ) + a( 2, 1 ) * b( 1, 0 ) + a( 2, 2 ) * b( 2, 0 ), 00146 a( 2, 0 ) * b( 0, 1 ) + a( 2, 1 ) * b( 1, 1 ) + a( 2, 2 ) * b( 2, 1 ), 00147 a( 2, 0 ) * b( 0, 2 ) + a( 2, 1 ) * b( 1, 2 ) + a( 2, 2 ) * b( 2, 2 ) ); 00148 } 00149 00150 template < typename Matrix > 00151 inline const Matrix inverse( const Matrix& d ) 00152 { 00153 const double det = 1.0 / determinant3( d ); 00154 return inverse( d, det ); 00155 } 00156 00157 template < typename Vector, typename Matrix > 00158 inline Vector vector_matrix( const Vector& v, const Matrix& m ) 00159 { 00160 return Vector( v[0] * m( 0, 0 ) + v[1] * m( 1, 0 ) + v[2] * m( 2, 0 ), 00161 v[0] * m( 0, 1 ) + v[1] * m( 1, 1 ) + v[2] * m( 2, 1 ), 00162 v[0] * m( 0, 2 ) + v[1] * m( 1, 2 ) + v[2] * m( 2, 2 ) ); 00163 } 00164 00165 template < typename Vector, typename Matrix > 00166 inline Vector matrix_vector( const Matrix& m, const Vector& v ) 00167 { 00168 Vector res = v; 00169 res[0] = v[0] * m( 0, 0 ) + v[1] * m( 0, 1 ) + v[2] * m( 0, 2 ); 00170 res[1] = v[0] * m( 1, 0 ) + v[1] * m( 1, 1 ) + v[2] * m( 1, 2 ); 00171 res[2] = v[0] * m( 2, 0 ) + v[1] * m( 2, 1 ) + v[2] * m( 2, 2 ); 00172 return res; 00173 } 00174 00175 } // namespace Matrix 00176 00177 class Matrix3 00178 { 00179 00180 public: 00181 const static int size = 9; 00182 00183 private: 00184 #ifndef MOAB_HAVE_LAPACK 00185 Eigen::Matrix3d _mat; 00186 #else 00187 double _mat[size]; 00188 #endif 00189 00190 public: 00191 // Default Constructor 00192 inline Matrix3() 00193 { 00194 #ifndef MOAB_HAVE_LAPACK 00195 _mat.fill( 0.0 ); 00196 #else 00197 MOAB_DMEMZERO( _mat, Matrix3::size ); 00198 #endif 00199 } 00200 00201 #ifndef MOAB_HAVE_LAPACK 00202 inline Matrix3( Eigen::Matrix3d mat ) : _mat( mat ) {} 00203 #endif 00204 00205 // TODO: Deprecate this. 00206 // Then we can go from three Constructors to one. 00207 inline Matrix3( double diagonal ) 00208 { 00209 #ifndef MOAB_HAVE_LAPACK 00210 _mat << diagonal, 0.0, 0.0, 0.0, diagonal, 0.0, 0.0, 0.0, diagonal; 00211 #else 00212 MOAB_DMEMZERO( _mat, Matrix3::size ); 00213 _mat[0] = _mat[4] = _mat[8] = diagonal; 00214 #endif 00215 } 00216 00217 inline Matrix3( const CartVect& diagonal ) 00218 { 00219 #ifndef MOAB_HAVE_LAPACK 00220 _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2]; 00221 #else 00222 MOAB_DMEMZERO( _mat, Matrix3::size ); 00223 _mat[0] = diagonal[0]; 00224 _mat[4] = diagonal[1]; 00225 _mat[8] = diagonal[2]; 00226 #endif 00227 } 00228 00229 // TODO: not strictly correct as the Matrix3 object 00230 // is a double d[ 9] so the only valid model of T is 00231 // double, or any refinement (int, float) 00232 //*but* it doesn't really matter anything else 00233 // will fail to compile. 00234 inline Matrix3( const std::vector< double >& diagonal ) 00235 { 00236 #ifndef MOAB_HAVE_LAPACK 00237 _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2]; 00238 #else 00239 MOAB_DMEMZERO( _mat, Matrix3::size ); 00240 _mat[0] = diagonal[0]; 00241 _mat[4] = diagonal[1]; 00242 _mat[8] = diagonal[2]; 00243 #endif 00244 } 00245 00246 inline Matrix3( double v00, double v01, double v02, double v10, double v11, double v12, double v20, double v21, 00247 double v22 ) 00248 { 00249 #ifndef MOAB_HAVE_LAPACK 00250 _mat << v00, v01, v02, v10, v11, v12, v20, v21, v22; 00251 #else 00252 MOAB_DMEMZERO( _mat, Matrix3::size ); 00253 _mat[0] = v00; 00254 _mat[1] = v01; 00255 _mat[2] = v02; 00256 _mat[3] = v10; 00257 _mat[4] = v11; 00258 _mat[5] = v12; 00259 _mat[6] = v20; 00260 _mat[7] = v21; 00261 _mat[8] = v22; 00262 #endif 00263 } 00264 00265 // Copy constructor 00266 Matrix3( const Matrix3& f ) 00267 { 00268 #ifndef MOAB_HAVE_LAPACK 00269 _mat = f._mat; 00270 #else 00271 memcpy( _mat, f._mat, size * sizeof( double ) ); 00272 #endif 00273 } 00274 00275 // Weird constructors 00276 template < typename Vector > 00277 inline Matrix3( const Vector& row0, const Vector& row1, const Vector& row2, const bool isRow ) 00278 { 00279 #ifndef MOAB_HAVE_LAPACK 00280 if( isRow ) { _mat << row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], row2[0], row2[1], row2[2]; } 00281 else 00282 { 00283 _mat << row0[0], row1[0], row2[0], row0[1], row1[1], row2[1], row0[2], row1[2], row2[2]; 00284 } 00285 #else 00286 MOAB_DMEMZERO( _mat, Matrix3::size ); 00287 if( isRow ) 00288 { 00289 _mat[0] = row0[0]; 00290 _mat[1] = row0[1]; 00291 _mat[2] = row0[2]; 00292 _mat[3] = row1[0]; 00293 _mat[4] = row1[1]; 00294 _mat[5] = row1[2]; 00295 _mat[6] = row2[0]; 00296 _mat[7] = row2[1]; 00297 _mat[8] = row2[2]; 00298 } 00299 else 00300 { 00301 _mat[0] = row0[0]; 00302 _mat[1] = row1[0]; 00303 _mat[2] = row2[0]; 00304 _mat[3] = row0[1]; 00305 _mat[4] = row1[1]; 00306 _mat[5] = row2[1]; 00307 _mat[6] = row0[2]; 00308 _mat[7] = row1[2]; 00309 _mat[8] = row2[2]; 00310 } 00311 #endif 00312 } 00313 00314 #ifndef DEPRECATED 00315 #ifdef __GNUC__ 00316 #define DEPRECATED __attribute__( ( deprecated ) ) 00317 #else 00318 #pragma message( "WARNING: You need to implement DEPRECATED for this compiler" ) 00319 #define DEPRECATED 00320 #endif 00321 #endif 00322 00323 /* 00324 * \deprecated { Use instead the constructor with explicit fourth argument, bool isRow, above } 00325 * 00326 */ 00327 inline Matrix3( const double v[size] ) 00328 { 00329 #ifndef MOAB_HAVE_LAPACK 00330 _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]; 00331 #else 00332 memcpy( _mat, v, size * sizeof( double ) ); 00333 #endif 00334 } 00335 00336 inline void copyto( double v[Matrix3::size] ) 00337 { 00338 #ifndef MOAB_HAVE_LAPACK 00339 std::copy( _mat.data(), _mat.data() + size, v ); 00340 #else 00341 memcpy( v, _mat, size * sizeof( double ) ); 00342 #endif 00343 } 00344 00345 inline Matrix3& operator=( const Matrix3& m ) 00346 { 00347 #ifndef MOAB_HAVE_LAPACK 00348 _mat = m._mat; 00349 #else 00350 memcpy( _mat, m._mat, size * sizeof( double ) ); 00351 #endif 00352 return *this; 00353 } 00354 00355 inline Matrix3& operator=( const double v[size] ) 00356 { 00357 #ifndef MOAB_HAVE_LAPACK 00358 _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8]; 00359 #else 00360 memcpy( _mat, v, size * sizeof( double ) ); 00361 #endif 00362 return *this; 00363 } 00364 00365 inline double* operator[]( unsigned i ) 00366 { 00367 #ifndef MOAB_HAVE_LAPACK 00368 return _mat.row( i ).data(); 00369 #else 00370 return &_mat[i * 3]; // Row Major 00371 #endif 00372 } 00373 00374 inline const double* operator[]( unsigned i ) const 00375 { 00376 #ifndef MOAB_HAVE_LAPACK 00377 return _mat.row( i ).data(); 00378 #else 00379 return &_mat[i * 3]; 00380 #endif 00381 } 00382 00383 inline double& operator()( unsigned r, unsigned c ) 00384 { 00385 #ifndef MOAB_HAVE_LAPACK 00386 return _mat( r, c ); 00387 #else 00388 return _mat[r * 3 + c]; 00389 #endif 00390 } 00391 00392 inline double operator()( unsigned r, unsigned c ) const 00393 { 00394 #ifndef MOAB_HAVE_LAPACK 00395 return _mat( r, c ); 00396 #else 00397 return _mat[r * 3 + c]; 00398 #endif 00399 } 00400 00401 inline double& operator()( unsigned i ) 00402 { 00403 #ifndef MOAB_HAVE_LAPACK 00404 return _mat( i ); 00405 #else 00406 return _mat[i]; 00407 #endif 00408 } 00409 00410 inline double operator()( unsigned i ) const 00411 { 00412 #ifndef MOAB_HAVE_LAPACK 00413 return _mat( i ); 00414 #else 00415 return _mat[i]; 00416 #endif 00417 } 00418 00419 // get pointer to array of nine doubles 00420 inline double* array() 00421 { 00422 #ifndef MOAB_HAVE_LAPACK 00423 return _mat.data(); 00424 #else 00425 return _mat; 00426 #endif 00427 } 00428 00429 inline const double* array() const 00430 { 00431 #ifndef MOAB_HAVE_LAPACK 00432 return _mat.data(); 00433 #else 00434 return _mat; 00435 #endif 00436 } 00437 00438 inline Matrix3& operator+=( const Matrix3& m ) 00439 { 00440 #ifndef MOAB_HAVE_LAPACK 00441 _mat += m._mat; 00442 #else 00443 for( int i = 0; i < Matrix3::size; ++i ) 00444 _mat[i] += m._mat[i]; 00445 #endif 00446 return *this; 00447 } 00448 00449 inline Matrix3& operator-=( const Matrix3& m ) 00450 { 00451 #ifndef MOAB_HAVE_LAPACK 00452 _mat -= m._mat; 00453 #else 00454 for( int i = 0; i < Matrix3::size; ++i ) 00455 _mat[i] -= m._mat[i]; 00456 #endif 00457 return *this; 00458 } 00459 00460 inline Matrix3& operator*=( double s ) 00461 { 00462 #ifndef MOAB_HAVE_LAPACK 00463 _mat *= s; 00464 #else 00465 for( int i = 0; i < Matrix3::size; ++i ) 00466 _mat[i] *= s; 00467 #endif 00468 return *this; 00469 } 00470 00471 inline Matrix3& operator/=( double s ) 00472 { 00473 #ifndef MOAB_HAVE_LAPACK 00474 _mat /= s; 00475 #else 00476 for( int i = 0; i < Matrix3::size; ++i ) 00477 _mat[i] /= s; 00478 #endif 00479 return *this; 00480 } 00481 00482 inline Matrix3& operator*=( const Matrix3& m ) 00483 { 00484 #ifndef MOAB_HAVE_LAPACK 00485 _mat *= m._mat; 00486 #else 00487 // Uncomment below if you want point-wise multiplication instead (.*) 00488 // for (int i=0; i < Matrix3::size; ++i) _mat[i] *= m._mat[i]; 00489 std::vector< double > dmat; 00490 dmat.assign( _mat, _mat + size ); 00491 _mat[0] = dmat[0] * m._mat[0] + dmat[1] * m._mat[3] + dmat[2] * m._mat[6]; 00492 _mat[1] = dmat[0] * m._mat[1] + dmat[1] * m._mat[4] + dmat[2] * m._mat[7]; 00493 _mat[2] = dmat[0] * m._mat[2] + dmat[1] * m._mat[5] + dmat[2] * m._mat[8]; 00494 _mat[3] = dmat[3] * m._mat[0] + dmat[4] * m._mat[3] + dmat[5] * m._mat[6]; 00495 _mat[4] = dmat[3] * m._mat[1] + dmat[4] * m._mat[4] + dmat[5] * m._mat[7]; 00496 _mat[5] = dmat[3] * m._mat[2] + dmat[4] * m._mat[5] + dmat[5] * m._mat[8]; 00497 _mat[6] = dmat[6] * m._mat[0] + dmat[7] * m._mat[3] + dmat[8] * m._mat[6]; 00498 _mat[7] = dmat[6] * m._mat[1] + dmat[7] * m._mat[4] + dmat[8] * m._mat[7]; 00499 _mat[8] = dmat[6] * m._mat[2] + dmat[7] * m._mat[5] + dmat[8] * m._mat[8]; 00500 #endif 00501 return *this; 00502 } 00503 00504 inline bool is_symmetric() 00505 { 00506 const double EPS = 1e-13; 00507 #ifndef MOAB_HAVE_LAPACK 00508 if( ( fabs( _mat( 1 ) - _mat( 3 ) ) < EPS ) && ( fabs( _mat( 2 ) - _mat( 6 ) ) < EPS ) && 00509 ( fabs( _mat( 5 ) - _mat( 7 ) ) < EPS ) ) 00510 return true; 00511 #else 00512 if( ( fabs( _mat[1] - _mat[3] ) < EPS ) && ( fabs( _mat[2] - _mat[6] ) < EPS ) && 00513 ( fabs( _mat[5] - _mat[7] ) < EPS ) ) 00514 return true; 00515 #endif 00516 else 00517 return false; 00518 } 00519 00520 inline bool is_positive_definite() 00521 { 00522 #ifndef MOAB_HAVE_LAPACK 00523 double subdet6 = _mat( 1 ) * _mat( 5 ) - _mat( 2 ) * _mat( 4 ); 00524 double subdet7 = _mat( 2 ) * _mat( 3 ) - _mat( 0 ) * _mat( 5 ); 00525 double subdet8 = _mat( 0 ) * _mat( 4 ) - _mat( 1 ) * _mat( 3 ); 00526 // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8; 00527 const double det = _mat( 6 ) * subdet6 + _mat( 7 ) * subdet7 + _mat( 8 ) * subdet8; 00528 return _mat( 0 ) > 0 && subdet8 > 0 && det > 0; 00529 #else 00530 double subdet6 = _mat[1] * _mat[5] - _mat[2] * _mat[4]; 00531 double subdet7 = _mat[2] * _mat[3] - _mat[0] * _mat[5]; 00532 double subdet8 = _mat[0] * _mat[4] - _mat[1] * _mat[3]; 00533 // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8; 00534 const double det = _mat[6] * subdet6 + _mat[7] * subdet7 + _mat[8] * subdet8; 00535 return _mat[0] > 0 && subdet8 > 0 && det > 0; 00536 #endif 00537 } 00538 00539 template < typename Vector > 00540 inline ErrorCode eigen_decomposition( Vector& evals, Matrix3& evecs ) 00541 { 00542 const bool bisSymmetric = this->is_symmetric(); 00543 #ifndef MOAB_HAVE_LAPACK 00544 if( bisSymmetric ) 00545 { 00546 Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d > eigensolver( this->_mat ); 00547 if( eigensolver.info() != Eigen::Success ) return MB_FAILURE; 00548 const Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d >::RealVectorType& e3evals = eigensolver.eigenvalues(); 00549 evals[0] = e3evals( 0 ); 00550 evals[1] = e3evals( 1 ); 00551 evals[2] = e3evals( 2 ); 00552 evecs._mat = eigensolver.eigenvectors(); //.col(1) 00553 return MB_SUCCESS; 00554 } 00555 else 00556 { 00557 MB_CHK_SET_ERR( MB_FAILURE, "Unsymmetric matrix implementation with Eigen3 is currently not provided." ); 00558 // Eigen::EigenSolver<Eigen::Matrix3d> eigensolver(this->_mat, true); 00559 // if (eigensolver.info() != Eigen::Success) 00560 // return MB_FAILURE; 00561 // const Eigen::EigenSolver<Eigen::Matrix3d>::EigenvalueType& e3evals = 00562 // eigensolver.eigenvalues().real(); evals[0] = e3evals(0); evals[1] = e3evals(1); 00563 // evals[2] = e3evals(2); evecs._mat = eigensolver.eigenvectors().real(); //.col(1) 00564 // return MB_SUCCESS; 00565 } 00566 #else 00567 int info; 00568 /* Solve eigenproblem */ 00569 double devreal[3], drevecs[9]; 00570 if( !bisSymmetric ) 00571 { 00572 double devimag[3], dlevecs[9], dwork[102]; 00573 char dgeev_opts[2] = { 'N', 'V' }; 00574 int N = 3, LWORK = 102, NL = 1, NR = N; 00575 std::vector< double > devmat; 00576 devmat.assign( _mat, _mat + size ); 00577 MOAB_dgeev( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, devimag, dlevecs, &NL, drevecs, 00578 &NR, dwork, &LWORK, &info ); 00579 // The result eigenvalues are ordered as high-->low 00580 evals[0] = devreal[2]; 00581 evals[1] = devreal[1]; 00582 evals[2] = devreal[0]; 00583 evecs._mat[0] = drevecs[6]; 00584 evecs._mat[1] = drevecs[3]; 00585 evecs._mat[2] = drevecs[0]; 00586 evecs._mat[3] = drevecs[7]; 00587 evecs._mat[4] = drevecs[4]; 00588 evecs._mat[5] = drevecs[1]; 00589 evecs._mat[6] = drevecs[8]; 00590 evecs._mat[7] = drevecs[5]; 00591 evecs._mat[8] = drevecs[2]; 00592 std::cout << "DGEEV: Optimal work vector: dsize = " << dwork[0] << ".\n"; 00593 } 00594 else 00595 { 00596 char dgeev_opts[2] = { 'V', 'L' }; 00597 const bool find_optimal = false; 00598 std::vector< int > iwork( 18 ); 00599 std::vector< double > devmat( 9, 0.0 ); 00600 std::vector< double > dwork( 38 ); 00601 int N = 3, lwork = 38, liwork = 18; 00602 devmat[0] = _mat[0]; 00603 devmat[1] = _mat[1]; 00604 devmat[2] = _mat[2]; 00605 devmat[4] = _mat[4]; 00606 devmat[5] = _mat[5]; 00607 devmat[8] = _mat[8]; 00608 if( find_optimal ) 00609 { 00610 int _lwork = -1; 00611 int _liwork = -1; 00612 double query_work_size = 0; 00613 int query_iwork_size = 0; 00614 // Make an empty call to find the optimal work vector size 00615 MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, NULL, &N, NULL, &query_work_size, &_lwork, 00616 &query_iwork_size, &_liwork, &info ); 00617 lwork = (int)query_work_size; 00618 dwork.resize( lwork ); 00619 liwork = query_iwork_size; 00620 iwork.resize( liwork ); 00621 std::cout << "DSYEVD: Optimal work vector: dsize = " << lwork << ", and isize = " << liwork << ".\n"; 00622 } 00623 00624 MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, &dwork[0], &lwork, &iwork[0], 00625 &liwork, &info ); 00626 for( int i = 0; i < 9; ++i ) 00627 drevecs[i] = devmat[i]; 00628 // The result eigenvalues are ordered as low-->high, but vectors are in rows of A. 00629 evals[0] = devreal[0]; 00630 evals[1] = devreal[1]; 00631 evals[2] = devreal[2]; 00632 evecs._mat[0] = drevecs[0]; 00633 evecs._mat[3] = drevecs[1]; 00634 evecs._mat[6] = drevecs[2]; 00635 evecs._mat[1] = drevecs[3]; 00636 evecs._mat[4] = drevecs[4]; 00637 evecs._mat[7] = drevecs[5]; 00638 evecs._mat[2] = drevecs[6]; 00639 evecs._mat[5] = drevecs[7]; 00640 evecs._mat[8] = drevecs[8]; 00641 } 00642 00643 if( !info ) { return MB_SUCCESS; } 00644 else 00645 { 00646 std::cout << "Failure in LAPACK_" << ( bisSymmetric ? "DSYEVD" : "DGEEV" ) 00647 << " call for eigen decomposition.\n"; 00648 std::cout << "Failed with error = " << info << ".\n"; 00649 return MB_FAILURE; 00650 } 00651 #endif 00652 } 00653 00654 inline void transpose_inplace() 00655 { 00656 #ifndef MOAB_HAVE_LAPACK 00657 _mat.transposeInPlace(); 00658 #else 00659 Matrix3 mtmp( *this ); 00660 _mat[1] = mtmp._mat[3]; 00661 _mat[3] = mtmp._mat[1]; 00662 _mat[2] = mtmp._mat[6]; 00663 _mat[6] = mtmp._mat[2]; 00664 _mat[5] = mtmp._mat[7]; 00665 _mat[7] = mtmp._mat[5]; 00666 #endif 00667 } 00668 00669 inline Matrix3 transpose() const 00670 { 00671 #ifndef MOAB_HAVE_LAPACK 00672 return Matrix3( _mat.transpose() ); 00673 #else 00674 Matrix3 mtmp( *this ); 00675 mtmp._mat[1] = _mat[3]; 00676 mtmp._mat[3] = _mat[1]; 00677 mtmp._mat[2] = _mat[6]; 00678 mtmp._mat[6] = _mat[2]; 00679 mtmp._mat[5] = _mat[7]; 00680 mtmp._mat[7] = _mat[5]; 00681 return mtmp; 00682 #endif 00683 } 00684 00685 template < typename Vector > 00686 inline void copycol( int index, Vector& vol ) 00687 { 00688 #ifndef MOAB_HAVE_LAPACK 00689 _mat.col( index ).swap( vol ); 00690 #else 00691 switch( index ) 00692 { 00693 case 0: 00694 _mat[0] = vol[0]; 00695 _mat[3] = vol[1]; 00696 _mat[6] = vol[2]; 00697 break; 00698 case 1: 00699 _mat[1] = vol[0]; 00700 _mat[4] = vol[1]; 00701 _mat[7] = vol[2]; 00702 break; 00703 case 2: 00704 _mat[2] = vol[0]; 00705 _mat[5] = vol[1]; 00706 _mat[8] = vol[2]; 00707 break; 00708 } 00709 #endif 00710 } 00711 00712 inline void swapcol( int srcindex, int destindex ) 00713 { 00714 assert( srcindex < Matrix3::size ); 00715 assert( destindex < Matrix3::size ); 00716 #ifndef MOAB_HAVE_LAPACK 00717 _mat.col( srcindex ).swap( _mat.col( destindex ) ); 00718 #else 00719 CartVect svol = this->vcol< CartVect >( srcindex ); 00720 CartVect dvol = this->vcol< CartVect >( destindex ); 00721 switch( srcindex ) 00722 { 00723 case 0: 00724 _mat[0] = dvol[0]; 00725 _mat[3] = dvol[1]; 00726 _mat[6] = dvol[2]; 00727 break; 00728 case 1: 00729 _mat[1] = dvol[0]; 00730 _mat[4] = dvol[1]; 00731 _mat[7] = dvol[2]; 00732 break; 00733 case 2: 00734 _mat[2] = dvol[0]; 00735 _mat[5] = dvol[1]; 00736 _mat[8] = dvol[2]; 00737 break; 00738 } 00739 switch( destindex ) 00740 { 00741 case 0: 00742 _mat[0] = svol[0]; 00743 _mat[3] = svol[1]; 00744 _mat[6] = svol[2]; 00745 break; 00746 case 1: 00747 _mat[1] = svol[0]; 00748 _mat[4] = svol[1]; 00749 _mat[7] = svol[2]; 00750 break; 00751 case 2: 00752 _mat[2] = svol[0]; 00753 _mat[5] = svol[1]; 00754 _mat[8] = svol[2]; 00755 break; 00756 } 00757 #endif 00758 } 00759 00760 template < typename Vector > 00761 inline Vector vcol( int index ) const 00762 { 00763 assert( index < Matrix3::size ); 00764 #ifndef MOAB_HAVE_LAPACK 00765 return _mat.col( index ); 00766 #else 00767 switch( index ) 00768 { 00769 case 0: 00770 return Vector( _mat[0], _mat[3], _mat[6] ); 00771 case 1: 00772 return Vector( _mat[1], _mat[4], _mat[7] ); 00773 case 2: 00774 return Vector( _mat[2], _mat[5], _mat[8] ); 00775 } 00776 return Vector( 0.0 ); 00777 #endif 00778 } 00779 00780 inline void colscale( int index, double scale ) 00781 { 00782 assert( index < Matrix3::size ); 00783 #ifndef MOAB_HAVE_LAPACK 00784 _mat.col( index ) *= scale; 00785 #else 00786 switch( index ) 00787 { 00788 case 0: 00789 _mat[0] *= scale; 00790 _mat[3] *= scale; 00791 _mat[6] *= scale; 00792 break; 00793 case 1: 00794 _mat[1] *= scale; 00795 _mat[4] *= scale; 00796 _mat[7] *= scale; 00797 break; 00798 case 2: 00799 _mat[2] *= scale; 00800 _mat[5] *= scale; 00801 _mat[8] *= scale; 00802 break; 00803 } 00804 #endif 00805 } 00806 00807 inline void rowscale( int index, double scale ) 00808 { 00809 assert( index < Matrix3::size ); 00810 #ifndef MOAB_HAVE_LAPACK 00811 _mat.row( index ) *= scale; 00812 #else 00813 switch( index ) 00814 { 00815 case 0: 00816 _mat[0] *= scale; 00817 _mat[1] *= scale; 00818 _mat[2] *= scale; 00819 break; 00820 case 1: 00821 _mat[3] *= scale; 00822 _mat[4] *= scale; 00823 _mat[5] *= scale; 00824 break; 00825 case 2: 00826 _mat[6] *= scale; 00827 _mat[7] *= scale; 00828 _mat[8] *= scale; 00829 break; 00830 } 00831 #endif 00832 } 00833 00834 inline CartVect col( int index ) const 00835 { 00836 assert( index < Matrix3::size ); 00837 #ifndef MOAB_HAVE_LAPACK 00838 Eigen::Vector3d mvec = _mat.col( index ); 00839 return CartVect( mvec[0], mvec[1], mvec[2] ); 00840 #else 00841 switch( index ) 00842 { 00843 case 0: 00844 return CartVect( _mat[0], _mat[3], _mat[6] ); 00845 case 1: 00846 return CartVect( _mat[1], _mat[4], _mat[7] ); 00847 case 2: 00848 return CartVect( _mat[2], _mat[5], _mat[8] ); 00849 } 00850 return CartVect( 0.0 ); 00851 #endif 00852 } 00853 00854 inline CartVect row( int index ) const 00855 { 00856 assert( index < Matrix3::size ); 00857 #ifndef MOAB_HAVE_LAPACK 00858 Eigen::Vector3d mvec = _mat.row( index ); 00859 return CartVect( mvec[0], mvec[1], mvec[2] ); 00860 #else 00861 switch( index ) 00862 { 00863 case 0: 00864 return CartVect( _mat[0], _mat[1], _mat[2] ); 00865 case 1: 00866 return CartVect( _mat[3], _mat[4], _mat[5] ); 00867 case 2: 00868 return CartVect( _mat[6], _mat[7], _mat[8] ); 00869 } 00870 return CartVect( 0.0 ); 00871 #endif 00872 } 00873 00874 friend Matrix3 operator+( const Matrix3& a, const Matrix3& b ); 00875 friend Matrix3 operator-( const Matrix3& a, const Matrix3& b ); 00876 friend Matrix3 operator*( const Matrix3& a, const Matrix3& b ); 00877 00878 inline double determinant() const 00879 { 00880 #ifndef MOAB_HAVE_LAPACK 00881 return _mat.determinant(); 00882 #else 00883 return ( _mat[0] * _mat[4] * _mat[8] + _mat[1] * _mat[5] * _mat[6] + _mat[2] * _mat[3] * _mat[7] - 00884 _mat[0] * _mat[5] * _mat[7] - _mat[1] * _mat[3] * _mat[8] - _mat[2] * _mat[4] * _mat[6] ); 00885 #endif 00886 } 00887 00888 inline Matrix3 inverse() const 00889 { 00890 #ifndef MOAB_HAVE_LAPACK 00891 return Matrix3( _mat.inverse() ); 00892 #else 00893 // return Matrix::compute_inverse( *this, this->determinant() ); 00894 Matrix3 m( 0.0 ); 00895 const double d_determinant = 1.0 / this->determinant(); 00896 m._mat[0] = d_determinant * ( _mat[4] * _mat[8] - _mat[5] * _mat[7] ); 00897 m._mat[1] = d_determinant * ( _mat[2] * _mat[7] - _mat[8] * _mat[1] ); 00898 m._mat[2] = d_determinant * ( _mat[1] * _mat[5] - _mat[4] * _mat[2] ); 00899 m._mat[3] = d_determinant * ( _mat[5] * _mat[6] - _mat[8] * _mat[3] ); 00900 m._mat[4] = d_determinant * ( _mat[0] * _mat[8] - _mat[6] * _mat[2] ); 00901 m._mat[5] = d_determinant * ( _mat[2] * _mat[3] - _mat[5] * _mat[0] ); 00902 m._mat[6] = d_determinant * ( _mat[3] * _mat[7] - _mat[6] * _mat[4] ); 00903 m._mat[7] = d_determinant * ( _mat[1] * _mat[6] - _mat[7] * _mat[0] ); 00904 m._mat[8] = d_determinant * ( _mat[0] * _mat[4] - _mat[3] * _mat[1] ); 00905 return m; 00906 #endif 00907 } 00908 00909 inline bool invert() 00910 { 00911 bool invertible = false; 00912 double d_determinant; 00913 #ifndef MOAB_HAVE_LAPACK 00914 Eigen::Matrix3d invMat; 00915 _mat.computeInverseAndDetWithCheck( invMat, d_determinant, invertible ); 00916 if( !Util::is_finite( d_determinant ) ) return false; 00917 _mat = invMat; 00918 return invertible; 00919 #else 00920 d_determinant = this->determinant(); 00921 if( d_determinant > 1e-13 ) invertible = true; 00922 d_determinant = 1.0 / d_determinant; // invert the determinant 00923 std::vector< double > _m; 00924 _m.assign( _mat, _mat + size ); 00925 _mat[0] = d_determinant * ( _m[4] * _m[8] - _m[5] * _m[7] ); 00926 _mat[1] = d_determinant * ( _m[2] * _m[7] - _m[8] * _m[1] ); 00927 _mat[2] = d_determinant * ( _m[1] * _m[5] - _m[4] * _m[2] ); 00928 _mat[3] = d_determinant * ( _m[5] * _m[6] - _m[8] * _m[3] ); 00929 _mat[4] = d_determinant * ( _m[0] * _m[8] - _m[6] * _m[2] ); 00930 _mat[5] = d_determinant * ( _m[2] * _m[3] - _m[5] * _m[0] ); 00931 _mat[6] = d_determinant * ( _m[3] * _m[7] - _m[6] * _m[4] ); 00932 _mat[7] = d_determinant * ( _m[1] * _m[6] - _m[7] * _m[0] ); 00933 _mat[8] = d_determinant * ( _m[0] * _m[4] - _m[3] * _m[1] ); 00934 #endif 00935 return invertible; 00936 } 00937 00938 // Calculate determinant of 2x2 submatrix composed of the 00939 // elements not in the passed row or column. 00940 inline double subdet( int r, int c ) const 00941 { 00942 assert( r >= 0 && c >= 0 ); 00943 if( r < 0 || c < 0 ) return DBL_MAX; 00944 #ifndef MOAB_HAVE_LAPACK 00945 const int r1 = ( r + 1 ) % 3, r2 = ( r + 2 ) % 3; 00946 const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3; 00947 return _mat( r1, c1 ) * _mat( r2, c2 ) - _mat( r1, c2 ) * _mat( r2, c1 ); 00948 #else 00949 const int r1 = Matrix3::size * ( ( r + 1 ) % 3 ), r2 = Matrix3::size * ( ( r + 2 ) % 3 ); 00950 const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3; 00951 return _mat[r1 + c1] * _mat[r2 + c2] - _mat[r1 + c2] * _mat[r2 + c1]; 00952 #endif 00953 } 00954 00955 inline void print( std::ostream& s ) const 00956 { 00957 #ifndef MOAB_HAVE_LAPACK 00958 s << "| " << _mat( 0 ) << " " << _mat( 1 ) << " " << _mat( 2 ) << " | " << _mat( 3 ) << " " << _mat( 4 ) << " " 00959 << _mat( 5 ) << " | " << _mat( 6 ) << " " << _mat( 7 ) << " " << _mat( 8 ) << " |"; 00960 #else 00961 s << "| " << _mat[0] << " " << _mat[1] << " " << _mat[2] << " | " << _mat[3] << " " << _mat[4] << " " << _mat[5] 00962 << " | " << _mat[6] << " " << _mat[7] << " " << _mat[8] << " |"; 00963 #endif 00964 } 00965 00966 }; // class Matrix3 00967 00968 template < typename Vector > 00969 inline Matrix3 outer_product( const Vector& u, const Vector& v ) 00970 { 00971 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], 00972 u[2] * v[1], u[2] * v[2] ); 00973 } 00974 00975 inline Matrix3 operator+( const Matrix3& a, const Matrix3& b ) 00976 { 00977 #ifndef MOAB_HAVE_LAPACK 00978 return Matrix3( a._mat + b._mat ); 00979 #else 00980 Matrix3 s( a ); 00981 for( int i = 0; i < Matrix3::size; ++i ) 00982 s( i ) += b._mat[i]; 00983 return s; 00984 #endif 00985 } 00986 00987 inline Matrix3 operator-( const Matrix3& a, const Matrix3& b ) 00988 { 00989 #ifndef MOAB_HAVE_LAPACK 00990 return Matrix3( a._mat - b._mat ); 00991 #else 00992 Matrix3 s( a ); 00993 for( int i = 0; i < Matrix3::size; ++i ) 00994 s( i ) -= b._mat[i]; 00995 return s; 00996 #endif 00997 } 00998 00999 inline Matrix3 operator*( const Matrix3& a, const Matrix3& b ) 01000 { 01001 #ifndef MOAB_HAVE_LAPACK 01002 return Matrix3( a._mat * b._mat ); 01003 #else 01004 return Matrix::mmult3( a, b ); 01005 #endif 01006 } 01007 01008 template < typename T > 01009 inline std::vector< T > operator*( const Matrix3& m, const std::vector< T >& v ) 01010 { 01011 return moab::Matrix::matrix_vector( m, v ); 01012 } 01013 01014 template < typename T > 01015 inline std::vector< T > operator*( const std::vector< T >& v, const Matrix3& m ) 01016 { 01017 return moab::Matrix::vector_matrix( v, m ); 01018 } 01019 01020 inline CartVect operator*( const Matrix3& m, const CartVect& v ) 01021 { 01022 return moab::Matrix::matrix_vector( m, v ); 01023 } 01024 01025 inline CartVect operator*( const CartVect& v, const Matrix3& m ) 01026 { 01027 return moab::Matrix::vector_matrix( v, m ); 01028 } 01029 01030 } // namespace moab 01031 01032 #ifdef MOAB_HAVE_LAPACK 01033 #undef MOAB_DMEMZERO 01034 #endif 01035 01036 #ifndef MOAB_MATRIX3_OPERATORLESS 01037 #define MOAB_MATRIX3_OPERATORLESS 01038 inline std::ostream& operator<<( std::ostream& s, const moab::Matrix3& m ) 01039 { 01040 return s << "| " << m( 0, 0 ) << " " << m( 0, 1 ) << " " << m( 0, 2 ) << " | " << m( 1, 0 ) << " " << m( 1, 1 ) 01041 << " " << m( 1, 2 ) << " | " << m( 2, 0 ) << " " << m( 2, 1 ) << " " << m( 2, 2 ) << " |"; 01042 } 01043 #endif // MOAB_MATRIX3_OPERATORLESS 01044 #endif // MOAB_MATRIX3_HPP