Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
Matrix3.hpp
Go to the documentation of this file.
1 /*
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  */
14 
15 /**\file Matrix3.hpp
16  *\author Jason Kraftcheck ([email protected])
17  *\date 2006-07-18
18  *\date 2012-08-2 Updated by rhl to be more generic. less code that does more!
19  * TODO: Remove all 'inline' keywords as it is only a suggestion to the compiler
20  * anyways, and it will ignore it or add it when it thinks its necessary.
21  *\date 2016-08-03 Updated to use Eigen3 support underneath to improve performance
22  */
23 
24 #ifndef MOAB_MATRIX3_HPP
25 #define MOAB_MATRIX3_HPP
26 
27 #include <iostream>
28 #include <iosfwd>
29 #include <limits>
30 #include <cmath>
31 #include <cassert>
32 
33 #include "moab/MOABConfig.h"
34 #include "moab/ErrorHandler.hpp"
35 #include "moab/Util.hpp"
36 #include "moab/Types.hpp"
37 #include "moab/CartVect.hpp"
38 
39 #ifndef MOAB_HAVE_LAPACK
40 
41 #ifndef MOAB_HAVE_EIGEN3
42 #error Need either Eigen3 or BLAS/LAPACK libraries
43 #endif
44 
45 #ifdef __GNUC__
46 // save diagnostic state
47 #pragma GCC diagnostic push
48 // turn off the specific warning. Can also use "-Wshadow"
49 #pragma GCC diagnostic ignored "-Wshadow"
50 #endif
51 
52 #define EIGEN_DEFAULT_TO_ROW_MAJOR
53 #define EIGEN_INITIALIZE_MATRICES_BY_ZERO
54 // #define EIGEN_NO_STATIC_ASSERT
55 #include "Eigen/Dense"
56 
57 #ifdef __GNUC__
58 // turn the warnings back on
59 #pragma GCC diagnostic pop
60 #endif
61 
62 #else
63 
64 #if defined( MOAB_FC_FUNC_ )
65 #define MOAB_FC_WRAPPER MOAB_FC_FUNC_
66 #elif defined( MOAB_FC_FUNC )
67 #define MOAB_FC_WRAPPER MOAB_FC_FUNC
68 #else
69 #define MOAB_FC_WRAPPER( name, NAME ) name##_
70 #endif
71 
72 // We will rely on LAPACK directly
73 
74 #ifdef WIN32
75 
76 // Should use second form below for windows but
77 // needed to do this to make it work.
78 // TODO: Need to clean this up
79 #define MOAB_dsyevd MOAB_FC_FUNC( dsyevd, DSYEVD )
80 #define MOAB_dsyevr MOAB_FC_FUNC( dsyevr, DSYEVR )
81 #define MOAB_dgeev MOAB_FC_FUNC( dgeev, DGEEV )
82 #define MOAB_dgetrf MOAB_FC_FUNC( dgetrf, DGETRF )
83 #define MOAB_dgetri MOAB_FC_FUNC( dgetri, DGETRI )
84 
85 #else
86 
87 #define MOAB_dsyevd MOAB_FC_WRAPPER( dsyevd, DSYEVD )
88 #define MOAB_dsyevr MOAB_FC_WRAPPER( dsyevr, DSYEVR )
89 #define MOAB_dgeev MOAB_FC_WRAPPER( dgeev, DGEEV )
90 #define MOAB_dgetrf MOAB_FC_WRAPPER( dgetrf, DGETRF )
91 #define MOAB_dgetri MOAB_FC_WRAPPER( dgetri, DGETRI )
92 
93 #endif
94 
95 extern "C" {
96 
97 // Computes all eigenvalues and, optionally, eigenvectors of a
98 // real symmetric matrix A. If eigenvectors are desired, it uses a
99 // divide and conquer algorithm.
100 void MOAB_dsyevd( char* jobz,
101  char* uplo,
102  int* n,
103  double a[],
104  int* lda,
105  double w[],
106  double work[],
107  int* lwork,
108  int iwork[],
109  int* liwork,
110  int* info );
111 
112 // Computes selected eigenvalues and, optionally, eigenvectors
113 // of a real symmetric matrix A. Eigenvalues and eigenvectors can be
114 // selected by specifying either a range of values or a range of
115 // indices for the desired eigenvalues.
116 void MOAB_dsyevr( char* jobz,
117  char* range,
118  char* uplo,
119  int* n,
120  double* a,
121  int* lda,
122  double* vl,
123  double* vu,
124  int* il,
125  int* iu,
126  double* abstol,
127  int* m,
128  double* w,
129  double* z,
130  int* ldz,
131  int* isuppz,
132  double* work,
133  int* lwork,
134  int* iwork,
135  int* liwork,
136  int* info );
137 
138 // Computes for an N-by-N real nonsymmetric matrix A, the
139 // eigenvalues and, optionally, the left and/or right eigenvectors.
140 void MOAB_dgeev( char* jobvl,
141  char* jobvr,
142  int* n,
143  double* a,
144  int* lda,
145  double* wr,
146  double* wi,
147  double* vl,
148  int* ldvl,
149  double* vr,
150  int* ldvr,
151  double* work,
152  int* lwork,
153  int* info );
154 
155 // Computes an LU factorization of a general M-by-N matrix A
156 // using partial pivoting with row interchanges.
157 void MOAB_dgetrf( int* M, int* N, double* A, int* lda, int* IPIV, int* INFO );
158 
159 // Computes the inverse of a matrix using the LU factorization
160 // computed by DGETRF.
161 void MOAB_dgetri( int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO );
162 }
163 
164 #include <cstring>
165 #define MOAB_DMEMZERO( a, b ) memset( a, 0, ( b ) * sizeof( double ) )
166 
167 #endif
168 
169 namespace moab
170 {
171 
172 namespace Matrix
173 {
174 
175  template < typename Matrix >
176  inline Matrix mmult3( const Matrix& a, const Matrix& b )
177  {
178  return Matrix( a( 0, 0 ) * b( 0, 0 ) + a( 0, 1 ) * b( 1, 0 ) + a( 0, 2 ) * b( 2, 0 ),
179  a( 0, 0 ) * b( 0, 1 ) + a( 0, 1 ) * b( 1, 1 ) + a( 0, 2 ) * b( 2, 1 ),
180  a( 0, 0 ) * b( 0, 2 ) + a( 0, 1 ) * b( 1, 2 ) + a( 0, 2 ) * b( 2, 2 ),
181  a( 1, 0 ) * b( 0, 0 ) + a( 1, 1 ) * b( 1, 0 ) + a( 1, 2 ) * b( 2, 0 ),
182  a( 1, 0 ) * b( 0, 1 ) + a( 1, 1 ) * b( 1, 1 ) + a( 1, 2 ) * b( 2, 1 ),
183  a( 1, 0 ) * b( 0, 2 ) + a( 1, 1 ) * b( 1, 2 ) + a( 1, 2 ) * b( 2, 2 ),
184  a( 2, 0 ) * b( 0, 0 ) + a( 2, 1 ) * b( 1, 0 ) + a( 2, 2 ) * b( 2, 0 ),
185  a( 2, 0 ) * b( 0, 1 ) + a( 2, 1 ) * b( 1, 1 ) + a( 2, 2 ) * b( 2, 1 ),
186  a( 2, 0 ) * b( 0, 2 ) + a( 2, 1 ) * b( 1, 2 ) + a( 2, 2 ) * b( 2, 2 ) );
187  }
188 
189  template < typename Matrix >
190  inline const Matrix inverse( const Matrix& d )
191  {
192  const double det = 1.0 / determinant3( d );
193  return inverse( d, det );
194  }
195 
196  template < typename Vector, typename Matrix >
197  inline Vector vector_matrix( const Vector& v, const Matrix& m )
198  {
199  return Vector( v[0] * m( 0, 0 ) + v[1] * m( 1, 0 ) + v[2] * m( 2, 0 ),
200  v[0] * m( 0, 1 ) + v[1] * m( 1, 1 ) + v[2] * m( 2, 1 ),
201  v[0] * m( 0, 2 ) + v[1] * m( 1, 2 ) + v[2] * m( 2, 2 ) );
202  }
203 
204  template < typename Vector, typename Matrix >
205  inline Vector matrix_vector( const Matrix& m, const Vector& v )
206  {
207  Vector res = v;
208  res[0] = v[0] * m( 0, 0 ) + v[1] * m( 0, 1 ) + v[2] * m( 0, 2 );
209  res[1] = v[0] * m( 1, 0 ) + v[1] * m( 1, 1 ) + v[2] * m( 1, 2 );
210  res[2] = v[0] * m( 2, 0 ) + v[1] * m( 2, 1 ) + v[2] * m( 2, 2 );
211  return res;
212  }
213 
214 } // namespace Matrix
215 
216 class Matrix3
217 {
218 
219  public:
220  const static int size = 9;
221 
222  private:
223 #ifndef MOAB_HAVE_LAPACK
224  Eigen::Matrix3d _mat;
225 #else
226  double _mat[size];
227 #endif
228 
229  public:
230  // Default Constructor
231  inline Matrix3()
232  {
233 #ifndef MOAB_HAVE_LAPACK
234  _mat.fill( 0.0 );
235 #else
236  MOAB_DMEMZERO( _mat, Matrix3::size );
237 #endif
238  }
239 
240 #ifndef MOAB_HAVE_LAPACK
241  inline Matrix3( Eigen::Matrix3d mat ) : _mat( mat ) {}
242 #endif
243 
244  // TODO: Deprecate this.
245  // Then we can go from three Constructors to one.
246  inline Matrix3( double diagonal )
247  {
248 #ifndef MOAB_HAVE_LAPACK
249  _mat << diagonal, 0.0, 0.0, 0.0, diagonal, 0.0, 0.0, 0.0, diagonal;
250 #else
251  MOAB_DMEMZERO( _mat, Matrix3::size );
252  _mat[0] = _mat[4] = _mat[8] = diagonal;
253 #endif
254  }
255 
256  inline Matrix3( const CartVect& diagonal )
257  {
258 #ifndef MOAB_HAVE_LAPACK
259  _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
260 #else
261  MOAB_DMEMZERO( _mat, Matrix3::size );
262  _mat[0] = diagonal[0];
263  _mat[4] = diagonal[1];
264  _mat[8] = diagonal[2];
265 #endif
266  }
267 
268  // TODO: not strictly correct as the Matrix3 object
269  // is a double d[ 9] so the only valid model of T is
270  // double, or any refinement (int, float)
271  //*but* it doesn't really matter anything else
272  // will fail to compile.
273  inline Matrix3( const std::vector< double >& diagonal )
274  {
275 #ifndef MOAB_HAVE_LAPACK
276  _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
277 #else
278  MOAB_DMEMZERO( _mat, Matrix3::size );
279  _mat[0] = diagonal[0];
280  _mat[4] = diagonal[1];
281  _mat[8] = diagonal[2];
282 #endif
283  }
284 
285  inline Matrix3( double v00,
286  double v01,
287  double v02,
288  double v10,
289  double v11,
290  double v12,
291  double v20,
292  double v21,
293  double v22 )
294  {
295 #ifndef MOAB_HAVE_LAPACK
296  _mat << v00, v01, v02, v10, v11, v12, v20, v21, v22;
297 #else
298  MOAB_DMEMZERO( _mat, Matrix3::size );
299  _mat[0] = v00;
300  _mat[1] = v01;
301  _mat[2] = v02;
302  _mat[3] = v10;
303  _mat[4] = v11;
304  _mat[5] = v12;
305  _mat[6] = v20;
306  _mat[7] = v21;
307  _mat[8] = v22;
308 #endif
309  }
310 
311  // Copy constructor
312  Matrix3( const Matrix3& f )
313  {
314 #ifndef MOAB_HAVE_LAPACK
315  _mat = f._mat;
316 #else
317  memcpy( _mat, f._mat, size * sizeof( double ) );
318 #endif
319  }
320 
321  // Weird constructors
322  template < typename Vector >
323  inline Matrix3( const Vector& row0, const Vector& row1, const Vector& row2, const bool isRow )
324  {
325 #ifndef MOAB_HAVE_LAPACK
326  if( isRow )
327  {
328  _mat << row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], row2[0], row2[1], row2[2];
329  }
330  else
331  {
332  _mat << row0[0], row1[0], row2[0], row0[1], row1[1], row2[1], row0[2], row1[2], row2[2];
333  }
334 #else
335  MOAB_DMEMZERO( _mat, Matrix3::size );
336  if( isRow )
337  {
338  _mat[0] = row0[0];
339  _mat[1] = row0[1];
340  _mat[2] = row0[2];
341  _mat[3] = row1[0];
342  _mat[4] = row1[1];
343  _mat[5] = row1[2];
344  _mat[6] = row2[0];
345  _mat[7] = row2[1];
346  _mat[8] = row2[2];
347  }
348  else
349  {
350  _mat[0] = row0[0];
351  _mat[1] = row1[0];
352  _mat[2] = row2[0];
353  _mat[3] = row0[1];
354  _mat[4] = row1[1];
355  _mat[5] = row2[1];
356  _mat[6] = row0[2];
357  _mat[7] = row1[2];
358  _mat[8] = row2[2];
359  }
360 #endif
361  }
362 
363 #ifndef DEPRECATED
364 #ifdef __GNUC__
365 #define DEPRECATED __attribute__( ( deprecated ) )
366 #else
367 #pragma message( "WARNING: You need to implement DEPRECATED for this compiler" )
368 #define DEPRECATED
369 #endif
370 #endif
371 
372  /*
373  * \deprecated { Use instead the constructor with explicit fourth argument, bool isRow, above }
374  *
375  */
376  inline Matrix3( const double v[size] )
377  {
378 #ifndef MOAB_HAVE_LAPACK
379  _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
380 #else
381  memcpy( _mat, v, size * sizeof( double ) );
382 #endif
383  }
384 
385  inline void copyto( double v[Matrix3::size] )
386  {
387 #ifndef MOAB_HAVE_LAPACK
388  std::copy( _mat.data(), _mat.data() + size, v );
389 #else
390  memcpy( v, _mat, size * sizeof( double ) );
391 #endif
392  }
393 
394  inline Matrix3& operator=( const Matrix3& m )
395  {
396 #ifndef MOAB_HAVE_LAPACK
397  _mat = m._mat;
398 #else
399  memcpy( _mat, m._mat, size * sizeof( double ) );
400 #endif
401  return *this;
402  }
403 
404  inline Matrix3& operator=( const double v[size] )
405  {
406 #ifndef MOAB_HAVE_LAPACK
407  _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
408 #else
409  memcpy( _mat, v, size * sizeof( double ) );
410 #endif
411  return *this;
412  }
413 
414  inline double* operator[]( unsigned i )
415  {
416 #ifndef MOAB_HAVE_LAPACK
417  return _mat.row( i ).data();
418 #else
419  return &_mat[i * 3]; // Row Major
420 #endif
421  }
422 
423  inline const double* operator[]( unsigned i ) const
424  {
425 #ifndef MOAB_HAVE_LAPACK
426  return _mat.row( i ).data();
427 #else
428  return &_mat[i * 3];
429 #endif
430  }
431 
432  inline double& operator()( unsigned r, unsigned c )
433  {
434 #ifndef MOAB_HAVE_LAPACK
435  return _mat( r, c );
436 #else
437  return _mat[r * 3 + c];
438 #endif
439  }
440 
441  inline double operator()( unsigned r, unsigned c ) const
442  {
443 #ifndef MOAB_HAVE_LAPACK
444  return _mat( r, c );
445 #else
446  return _mat[r * 3 + c];
447 #endif
448  }
449 
450  inline double& operator()( unsigned i )
451  {
452 #ifndef MOAB_HAVE_LAPACK
453  return _mat( i );
454 #else
455  return _mat[i];
456 #endif
457  }
458 
459  inline double operator()( unsigned i ) const
460  {
461 #ifndef MOAB_HAVE_LAPACK
462  return _mat( i );
463 #else
464  return _mat[i];
465 #endif
466  }
467 
468  // get pointer to array of nine doubles
469  inline double* array()
470  {
471 #ifndef MOAB_HAVE_LAPACK
472  return _mat.data();
473 #else
474  return _mat;
475 #endif
476  }
477 
478  inline const double* array() const
479  {
480 #ifndef MOAB_HAVE_LAPACK
481  return _mat.data();
482 #else
483  return _mat;
484 #endif
485  }
486 
487  inline Matrix3& operator+=( const Matrix3& m )
488  {
489 #ifndef MOAB_HAVE_LAPACK
490  _mat += m._mat;
491 #else
492  for( int i = 0; i < Matrix3::size; ++i )
493  _mat[i] += m._mat[i];
494 #endif
495  return *this;
496  }
497 
498  inline Matrix3& operator-=( const Matrix3& m )
499  {
500 #ifndef MOAB_HAVE_LAPACK
501  _mat -= m._mat;
502 #else
503  for( int i = 0; i < Matrix3::size; ++i )
504  _mat[i] -= m._mat[i];
505 #endif
506  return *this;
507  }
508 
509  inline Matrix3& operator*=( double s )
510  {
511 #ifndef MOAB_HAVE_LAPACK
512  _mat *= s;
513 #else
514  for( int i = 0; i < Matrix3::size; ++i )
515  _mat[i] *= s;
516 #endif
517  return *this;
518  }
519 
520  inline Matrix3& operator/=( double s )
521  {
522 #ifndef MOAB_HAVE_LAPACK
523  _mat /= s;
524 #else
525  for( int i = 0; i < Matrix3::size; ++i )
526  _mat[i] /= s;
527 #endif
528  return *this;
529  }
530 
531  inline Matrix3& operator*=( const Matrix3& m )
532  {
533 #ifndef MOAB_HAVE_LAPACK
534  _mat *= m._mat;
535 #else
536  // Uncomment below if you want point-wise multiplication instead (.*)
537  // for (int i=0; i < Matrix3::size; ++i) _mat[i] *= m._mat[i];
538  std::vector< double > dmat;
539  dmat.assign( _mat, _mat + size );
540  _mat[0] = dmat[0] * m._mat[0] + dmat[1] * m._mat[3] + dmat[2] * m._mat[6];
541  _mat[1] = dmat[0] * m._mat[1] + dmat[1] * m._mat[4] + dmat[2] * m._mat[7];
542  _mat[2] = dmat[0] * m._mat[2] + dmat[1] * m._mat[5] + dmat[2] * m._mat[8];
543  _mat[3] = dmat[3] * m._mat[0] + dmat[4] * m._mat[3] + dmat[5] * m._mat[6];
544  _mat[4] = dmat[3] * m._mat[1] + dmat[4] * m._mat[4] + dmat[5] * m._mat[7];
545  _mat[5] = dmat[3] * m._mat[2] + dmat[4] * m._mat[5] + dmat[5] * m._mat[8];
546  _mat[6] = dmat[6] * m._mat[0] + dmat[7] * m._mat[3] + dmat[8] * m._mat[6];
547  _mat[7] = dmat[6] * m._mat[1] + dmat[7] * m._mat[4] + dmat[8] * m._mat[7];
548  _mat[8] = dmat[6] * m._mat[2] + dmat[7] * m._mat[5] + dmat[8] * m._mat[8];
549 #endif
550  return *this;
551  }
552 
553  inline bool is_symmetric()
554  {
555  const double EPS = 1e-13;
556 #ifndef MOAB_HAVE_LAPACK
557  if( ( fabs( _mat( 1 ) - _mat( 3 ) ) < EPS ) && ( fabs( _mat( 2 ) - _mat( 6 ) ) < EPS ) &&
558  ( fabs( _mat( 5 ) - _mat( 7 ) ) < EPS ) )
559  return true;
560 #else
561  if( ( fabs( _mat[1] - _mat[3] ) < EPS ) && ( fabs( _mat[2] - _mat[6] ) < EPS ) &&
562  ( fabs( _mat[5] - _mat[7] ) < EPS ) )
563  return true;
564 #endif
565  else
566  return false;
567  }
568 
569  inline bool is_positive_definite()
570  {
571 #ifndef MOAB_HAVE_LAPACK
572  double subdet6 = _mat( 1 ) * _mat( 5 ) - _mat( 2 ) * _mat( 4 );
573  double subdet7 = _mat( 2 ) * _mat( 3 ) - _mat( 0 ) * _mat( 5 );
574  double subdet8 = _mat( 0 ) * _mat( 4 ) - _mat( 1 ) * _mat( 3 );
575  // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
576  const double det = _mat( 6 ) * subdet6 + _mat( 7 ) * subdet7 + _mat( 8 ) * subdet8;
577  return _mat( 0 ) > 0 && subdet8 > 0 && det > 0;
578 #else
579  double subdet6 = _mat[1] * _mat[5] - _mat[2] * _mat[4];
580  double subdet7 = _mat[2] * _mat[3] - _mat[0] * _mat[5];
581  double subdet8 = _mat[0] * _mat[4] - _mat[1] * _mat[3];
582  // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
583  const double det = _mat[6] * subdet6 + _mat[7] * subdet7 + _mat[8] * subdet8;
584  return _mat[0] > 0 && subdet8 > 0 && det > 0;
585 #endif
586  }
587 
588  template < typename Vector >
589  inline ErrorCode eigen_decomposition( Vector& evals, Matrix3& evecs )
590  {
591  const bool bisSymmetric = this->is_symmetric();
592 #ifndef MOAB_HAVE_LAPACK
593  if( bisSymmetric )
594  {
595  Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d > eigensolver( this->_mat );
596  if( eigensolver.info() != Eigen::Success ) return MB_FAILURE;
597  const Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d >::RealVectorType& e3evals = eigensolver.eigenvalues();
598  evals[0] = e3evals( 0 );
599  evals[1] = e3evals( 1 );
600  evals[2] = e3evals( 2 );
601  evecs._mat = eigensolver.eigenvectors(); //.col(1)
602  return MB_SUCCESS;
603  }
604  else
605  {
606  MB_CHK_SET_ERR( MB_FAILURE, "Unsymmetric matrix implementation with Eigen3 is currently not provided." );
607  // Eigen::EigenSolver<Eigen::Matrix3d> eigensolver(this->_mat, true);
608  // if (eigensolver.info() != Eigen::Success)
609  // return MB_FAILURE;
610  // const Eigen::EigenSolver<Eigen::Matrix3d>::EigenvalueType& e3evals =
611  // eigensolver.eigenvalues().real(); evals[0] = e3evals(0); evals[1] = e3evals(1);
612  // evals[2] = e3evals(2); evecs._mat = eigensolver.eigenvectors().real(); //.col(1)
613  // return MB_SUCCESS;
614  }
615 #else
616  int info;
617  /* Solve eigenproblem */
618  double devreal[3], drevecs[9];
619  if( !bisSymmetric )
620  {
621  double devimag[3], dlevecs[9], dwork[102];
622  char dgeev_opts[2] = { 'N', 'V' };
623  int N = 3, LWORK = 102, NL = 1, NR = N;
624  std::vector< double > devmat;
625  devmat.assign( _mat, _mat + size );
626  MOAB_dgeev( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, devimag, dlevecs, &NL, drevecs,
627  &NR, dwork, &LWORK, &info );
628  // The result eigenvalues are ordered as high-->low
629  evals[0] = devreal[2];
630  evals[1] = devreal[1];
631  evals[2] = devreal[0];
632  evecs._mat[0] = drevecs[6];
633  evecs._mat[1] = drevecs[3];
634  evecs._mat[2] = drevecs[0];
635  evecs._mat[3] = drevecs[7];
636  evecs._mat[4] = drevecs[4];
637  evecs._mat[5] = drevecs[1];
638  evecs._mat[6] = drevecs[8];
639  evecs._mat[7] = drevecs[5];
640  evecs._mat[8] = drevecs[2];
641  std::cout << "DGEEV: Optimal work vector: dsize = " << dwork[0] << ".\n";
642  }
643  else
644  {
645  char dgeev_opts[2] = { 'V', 'L' };
646  const bool find_optimal = false;
647  std::vector< int > iwork( 18 );
648  std::vector< double > devmat( 9, 0.0 );
649  std::vector< double > dwork( 38 );
650  int N = 3, lwork = 38, liwork = 18;
651  devmat[0] = _mat[0];
652  devmat[1] = _mat[1];
653  devmat[2] = _mat[2];
654  devmat[4] = _mat[4];
655  devmat[5] = _mat[5];
656  devmat[8] = _mat[8];
657  if( find_optimal )
658  {
659  int _lwork = -1;
660  int _liwork = -1;
661  double query_work_size = 0;
662  int query_iwork_size = 0;
663  // Make an empty call to find the optimal work vector size
664  MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, NULL, &N, NULL, &query_work_size, &_lwork,
665  &query_iwork_size, &_liwork, &info );
666  lwork = (int)query_work_size;
667  dwork.resize( lwork );
668  liwork = query_iwork_size;
669  iwork.resize( liwork );
670  std::cout << "DSYEVD: Optimal work vector: dsize = " << lwork << ", and isize = " << liwork << ".\n";
671  }
672 
673  MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, &dwork[0], &lwork, &iwork[0],
674  &liwork, &info );
675  for( int i = 0; i < 9; ++i )
676  drevecs[i] = devmat[i];
677  // The result eigenvalues are ordered as low-->high, but vectors are in rows of A.
678  evals[0] = devreal[0];
679  evals[1] = devreal[1];
680  evals[2] = devreal[2];
681  evecs._mat[0] = drevecs[0];
682  evecs._mat[3] = drevecs[1];
683  evecs._mat[6] = drevecs[2];
684  evecs._mat[1] = drevecs[3];
685  evecs._mat[4] = drevecs[4];
686  evecs._mat[7] = drevecs[5];
687  evecs._mat[2] = drevecs[6];
688  evecs._mat[5] = drevecs[7];
689  evecs._mat[8] = drevecs[8];
690  }
691 
692  if( !info )
693  {
694  return MB_SUCCESS;
695  }
696  else
697  {
698  std::cout << "Failure in LAPACK_" << ( bisSymmetric ? "DSYEVD" : "DGEEV" )
699  << " call for eigen decomposition.\n";
700  std::cout << "Failed with error = " << info << ".\n";
701  return MB_FAILURE;
702  }
703 #endif
704  }
705 
706  inline void transpose_inplace()
707  {
708 #ifndef MOAB_HAVE_LAPACK
709  _mat.transposeInPlace();
710 #else
711  Matrix3 mtmp( *this );
712  _mat[1] = mtmp._mat[3];
713  _mat[3] = mtmp._mat[1];
714  _mat[2] = mtmp._mat[6];
715  _mat[6] = mtmp._mat[2];
716  _mat[5] = mtmp._mat[7];
717  _mat[7] = mtmp._mat[5];
718 #endif
719  }
720 
721  inline Matrix3 transpose() const
722  {
723 #ifndef MOAB_HAVE_LAPACK
724  return Matrix3( _mat.transpose() );
725 #else
726  Matrix3 mtmp( *this );
727  mtmp._mat[1] = _mat[3];
728  mtmp._mat[3] = _mat[1];
729  mtmp._mat[2] = _mat[6];
730  mtmp._mat[6] = _mat[2];
731  mtmp._mat[5] = _mat[7];
732  mtmp._mat[7] = _mat[5];
733  return mtmp;
734 #endif
735  }
736 
737  template < typename Vector >
738  inline void copycol( int index, Vector& vol )
739  {
740 #ifndef MOAB_HAVE_LAPACK
741  _mat.col( index ).swap( vol );
742 #else
743  switch( index )
744  {
745  case 0:
746  _mat[0] = vol[0];
747  _mat[3] = vol[1];
748  _mat[6] = vol[2];
749  break;
750  case 1:
751  _mat[1] = vol[0];
752  _mat[4] = vol[1];
753  _mat[7] = vol[2];
754  break;
755  case 2:
756  _mat[2] = vol[0];
757  _mat[5] = vol[1];
758  _mat[8] = vol[2];
759  break;
760  }
761 #endif
762  }
763 
764  inline void swapcol( int srcindex, int destindex )
765  {
766  assert( srcindex < Matrix3::size );
767  assert( destindex < Matrix3::size );
768 #ifndef MOAB_HAVE_LAPACK
769  _mat.col( srcindex ).swap( _mat.col( destindex ) );
770 #else
771  CartVect svol = this->vcol< CartVect >( srcindex );
772  CartVect dvol = this->vcol< CartVect >( destindex );
773  switch( srcindex )
774  {
775  case 0:
776  _mat[0] = dvol[0];
777  _mat[3] = dvol[1];
778  _mat[6] = dvol[2];
779  break;
780  case 1:
781  _mat[1] = dvol[0];
782  _mat[4] = dvol[1];
783  _mat[7] = dvol[2];
784  break;
785  case 2:
786  _mat[2] = dvol[0];
787  _mat[5] = dvol[1];
788  _mat[8] = dvol[2];
789  break;
790  }
791  switch( destindex )
792  {
793  case 0:
794  _mat[0] = svol[0];
795  _mat[3] = svol[1];
796  _mat[6] = svol[2];
797  break;
798  case 1:
799  _mat[1] = svol[0];
800  _mat[4] = svol[1];
801  _mat[7] = svol[2];
802  break;
803  case 2:
804  _mat[2] = svol[0];
805  _mat[5] = svol[1];
806  _mat[8] = svol[2];
807  break;
808  }
809 #endif
810  }
811 
812  template < typename Vector >
813  inline Vector vcol( int index ) const
814  {
815  assert( index < Matrix3::size );
816 #ifndef MOAB_HAVE_LAPACK
817  return _mat.col( index );
818 #else
819  switch( index )
820  {
821  case 0:
822  return Vector( _mat[0], _mat[3], _mat[6] );
823  case 1:
824  return Vector( _mat[1], _mat[4], _mat[7] );
825  case 2:
826  return Vector( _mat[2], _mat[5], _mat[8] );
827  }
828  return Vector( 0.0 );
829 #endif
830  }
831 
832  inline void colscale( int index, double scale )
833  {
834  assert( index < Matrix3::size );
835 #ifndef MOAB_HAVE_LAPACK
836  _mat.col( index ) *= scale;
837 #else
838  switch( index )
839  {
840  case 0:
841  _mat[0] *= scale;
842  _mat[3] *= scale;
843  _mat[6] *= scale;
844  break;
845  case 1:
846  _mat[1] *= scale;
847  _mat[4] *= scale;
848  _mat[7] *= scale;
849  break;
850  case 2:
851  _mat[2] *= scale;
852  _mat[5] *= scale;
853  _mat[8] *= scale;
854  break;
855  }
856 #endif
857  }
858 
859  inline void rowscale( int index, double scale )
860  {
861  assert( index < Matrix3::size );
862 #ifndef MOAB_HAVE_LAPACK
863  _mat.row( index ) *= scale;
864 #else
865  switch( index )
866  {
867  case 0:
868  _mat[0] *= scale;
869  _mat[1] *= scale;
870  _mat[2] *= scale;
871  break;
872  case 1:
873  _mat[3] *= scale;
874  _mat[4] *= scale;
875  _mat[5] *= scale;
876  break;
877  case 2:
878  _mat[6] *= scale;
879  _mat[7] *= scale;
880  _mat[8] *= scale;
881  break;
882  }
883 #endif
884  }
885 
886  inline CartVect col( int index ) const
887  {
888  assert( index < Matrix3::size );
889 #ifndef MOAB_HAVE_LAPACK
890  Eigen::Vector3d mvec = _mat.col( index );
891  return CartVect( mvec[0], mvec[1], mvec[2] );
892 #else
893  switch( index )
894  {
895  case 0:
896  return CartVect( _mat[0], _mat[3], _mat[6] );
897  case 1:
898  return CartVect( _mat[1], _mat[4], _mat[7] );
899  case 2:
900  return CartVect( _mat[2], _mat[5], _mat[8] );
901  }
902  return CartVect( 0.0 );
903 #endif
904  }
905 
906  inline CartVect row( int index ) const
907  {
908  assert( index < Matrix3::size );
909 #ifndef MOAB_HAVE_LAPACK
910  Eigen::Vector3d mvec = _mat.row( index );
911  return CartVect( mvec[0], mvec[1], mvec[2] );
912 #else
913  switch( index )
914  {
915  case 0:
916  return CartVect( _mat[0], _mat[1], _mat[2] );
917  case 1:
918  return CartVect( _mat[3], _mat[4], _mat[5] );
919  case 2:
920  return CartVect( _mat[6], _mat[7], _mat[8] );
921  }
922  return CartVect( 0.0 );
923 #endif
924  }
925 
926  friend Matrix3 operator+( const Matrix3& a, const Matrix3& b );
927  friend Matrix3 operator-( const Matrix3& a, const Matrix3& b );
928  friend Matrix3 operator*( const Matrix3& a, const Matrix3& b );
929 
930  inline double determinant() const
931  {
932 #ifndef MOAB_HAVE_LAPACK
933  return _mat.determinant();
934 #else
935  return ( _mat[0] * _mat[4] * _mat[8] + _mat[1] * _mat[5] * _mat[6] + _mat[2] * _mat[3] * _mat[7] -
936  _mat[0] * _mat[5] * _mat[7] - _mat[1] * _mat[3] * _mat[8] - _mat[2] * _mat[4] * _mat[6] );
937 #endif
938  }
939 
940  inline Matrix3 inverse() const
941  {
942 #ifndef MOAB_HAVE_LAPACK
943  return Matrix3( _mat.inverse() );
944 #else
945  // return Matrix::compute_inverse( *this, this->determinant() );
946  Matrix3 m( 0.0 );
947  const double d_determinant = 1.0 / this->determinant();
948  m._mat[0] = d_determinant * ( _mat[4] * _mat[8] - _mat[5] * _mat[7] );
949  m._mat[1] = d_determinant * ( _mat[2] * _mat[7] - _mat[8] * _mat[1] );
950  m._mat[2] = d_determinant * ( _mat[1] * _mat[5] - _mat[4] * _mat[2] );
951  m._mat[3] = d_determinant * ( _mat[5] * _mat[6] - _mat[8] * _mat[3] );
952  m._mat[4] = d_determinant * ( _mat[0] * _mat[8] - _mat[6] * _mat[2] );
953  m._mat[5] = d_determinant * ( _mat[2] * _mat[3] - _mat[5] * _mat[0] );
954  m._mat[6] = d_determinant * ( _mat[3] * _mat[7] - _mat[6] * _mat[4] );
955  m._mat[7] = d_determinant * ( _mat[1] * _mat[6] - _mat[7] * _mat[0] );
956  m._mat[8] = d_determinant * ( _mat[0] * _mat[4] - _mat[3] * _mat[1] );
957  return m;
958 #endif
959  }
960 
961  inline bool invert()
962  {
963  bool invertible = false;
964  double d_determinant;
965 #ifndef MOAB_HAVE_LAPACK
966  Eigen::Matrix3d invMat;
967  _mat.computeInverseAndDetWithCheck( invMat, d_determinant, invertible );
968  if( !Util::is_finite( d_determinant ) ) return false;
969  _mat = invMat;
970  return invertible;
971 #else
972  d_determinant = this->determinant();
973  if( d_determinant > 1e-13 ) invertible = true;
974  d_determinant = 1.0 / d_determinant; // invert the determinant
975  std::vector< double > _m;
976  _m.assign( _mat, _mat + size );
977  _mat[0] = d_determinant * ( _m[4] * _m[8] - _m[5] * _m[7] );
978  _mat[1] = d_determinant * ( _m[2] * _m[7] - _m[8] * _m[1] );
979  _mat[2] = d_determinant * ( _m[1] * _m[5] - _m[4] * _m[2] );
980  _mat[3] = d_determinant * ( _m[5] * _m[6] - _m[8] * _m[3] );
981  _mat[4] = d_determinant * ( _m[0] * _m[8] - _m[6] * _m[2] );
982  _mat[5] = d_determinant * ( _m[2] * _m[3] - _m[5] * _m[0] );
983  _mat[6] = d_determinant * ( _m[3] * _m[7] - _m[6] * _m[4] );
984  _mat[7] = d_determinant * ( _m[1] * _m[6] - _m[7] * _m[0] );
985  _mat[8] = d_determinant * ( _m[0] * _m[4] - _m[3] * _m[1] );
986 #endif
987  return invertible;
988  }
989 
990  // Calculate determinant of 2x2 submatrix composed of the
991  // elements not in the passed row or column.
992  inline double subdet( int r, int c ) const
993  {
994  assert( r >= 0 && c >= 0 );
995  if( r < 0 || c < 0 ) return DBL_MAX;
996 #ifndef MOAB_HAVE_LAPACK
997  const int r1 = ( r + 1 ) % 3, r2 = ( r + 2 ) % 3;
998  const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
999  return _mat( r1, c1 ) * _mat( r2, c2 ) - _mat( r1, c2 ) * _mat( r2, c1 );
1000 #else
1001  const int r1 = Matrix3::size * ( ( r + 1 ) % 3 ), r2 = Matrix3::size * ( ( r + 2 ) % 3 );
1002  const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
1003  return _mat[r1 + c1] * _mat[r2 + c2] - _mat[r1 + c2] * _mat[r2 + c1];
1004 #endif
1005  }
1006 
1007  inline void print( std::ostream& s ) const
1008  {
1009 #ifndef MOAB_HAVE_LAPACK
1010  s << "| " << _mat( 0 ) << " " << _mat( 1 ) << " " << _mat( 2 ) << " | " << _mat( 3 ) << " " << _mat( 4 ) << " "
1011  << _mat( 5 ) << " | " << _mat( 6 ) << " " << _mat( 7 ) << " " << _mat( 8 ) << " |";
1012 #else
1013  s << "| " << _mat[0] << " " << _mat[1] << " " << _mat[2] << " | " << _mat[3] << " " << _mat[4] << " " << _mat[5]
1014  << " | " << _mat[6] << " " << _mat[7] << " " << _mat[8] << " |";
1015 #endif
1016  }
1017 
1018 }; // class Matrix3
1019 
1020 template < typename Vector >
1021 inline Matrix3 outer_product( const Vector& u, const Vector& v )
1022 {
1023  return Matrix3( u[0] * v[0], u[0] * v[1], u[0] * v[2], u[1] * v[0], u[1] * v[1], u[1] * v[2], u[2] * v[0],
1024  u[2] * v[1], u[2] * v[2] );
1025 }
1026 
1027 inline Matrix3 operator+( const Matrix3& a, const Matrix3& b )
1028 {
1029 #ifndef MOAB_HAVE_LAPACK
1030  return Matrix3( a._mat + b._mat );
1031 #else
1032  Matrix3 s( a );
1033  for( int i = 0; i < Matrix3::size; ++i )
1034  s( i ) += b._mat[i];
1035  return s;
1036 #endif
1037 }
1038 
1039 inline Matrix3 operator-( const Matrix3& a, const Matrix3& b )
1040 {
1041 #ifndef MOAB_HAVE_LAPACK
1042  return Matrix3( a._mat - b._mat );
1043 #else
1044  Matrix3 s( a );
1045  for( int i = 0; i < Matrix3::size; ++i )
1046  s( i ) -= b._mat[i];
1047  return s;
1048 #endif
1049 }
1050 
1051 inline Matrix3 operator*( const Matrix3& a, const Matrix3& b )
1052 {
1053 #ifndef MOAB_HAVE_LAPACK
1054  return Matrix3( a._mat * b._mat );
1055 #else
1056  return Matrix::mmult3( a, b );
1057 #endif
1058 }
1059 
1060 template < typename T >
1061 inline std::vector< T > operator*( const Matrix3& m, const std::vector< T >& v )
1062 {
1063  return moab::Matrix::matrix_vector( m, v );
1064 }
1065 
1066 template < typename T >
1067 inline std::vector< T > operator*( const std::vector< T >& v, const Matrix3& m )
1068 {
1069  return moab::Matrix::vector_matrix( v, m );
1070 }
1071 
1072 inline CartVect operator*( const Matrix3& m, const CartVect& v )
1073 {
1074  return moab::Matrix::matrix_vector( m, v );
1075 }
1076 
1077 inline CartVect operator*( const CartVect& v, const Matrix3& m )
1078 {
1079  return moab::Matrix::vector_matrix( v, m );
1080 }
1081 
1082 } // namespace moab
1083 
1084 #ifdef MOAB_HAVE_LAPACK
1085 #undef MOAB_DMEMZERO
1086 #endif
1087 
1088 #ifndef MOAB_MATRIX3_OPERATORLESS
1089 #define MOAB_MATRIX3_OPERATORLESS
1090 inline std::ostream& operator<<( std::ostream& s, const moab::Matrix3& m )
1091 {
1092  return s << "| " << m( 0, 0 ) << " " << m( 0, 1 ) << " " << m( 0, 2 ) << " | " << m( 1, 0 ) << " " << m( 1, 1 )
1093  << " " << m( 1, 2 ) << " | " << m( 2, 0 ) << " " << m( 2, 1 ) << " " << m( 2, 2 ) << " |";
1094 }
1095 #endif // MOAB_MATRIX3_OPERATORLESS
1096 #endif // MOAB_MATRIX3_HPP