MOAB: Mesh Oriented datABase  (version 5.3.1)
Matrix3.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines