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