Branch data Line data Source code
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_EIGEN
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, char* uplo, int* n, double a[], int* lda, double w[], double work[], int* lwork,
101 : : int iwork[], int* liwork, int* info );
102 : :
103 : : // Computes selected eigenvalues and, optionally, eigenvectors
104 : : // of a real symmetric matrix A. Eigenvalues and eigenvectors can be
105 : : // selected by specifying either a range of values or a range of
106 : : // indices for the desired eigenvalues.
107 : : void MOAB_dsyevr( char* jobz, char* range, char* uplo, int* n, double* a, int* lda, double* vl, double* vu, int* il,
108 : : int* iu, double* abstol, int* m, double* w, double* z, int* ldz, int* isuppz, double* work,
109 : : int* lwork, int* iwork, int* liwork, int* info );
110 : :
111 : : // Computes for an N-by-N real nonsymmetric matrix A, the
112 : : // eigenvalues and, optionally, the left and/or right eigenvectors.
113 : : void MOAB_dgeev( char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl,
114 : : double* vr, int* ldvr, double* work, int* lwork, int* info );
115 : :
116 : : // Computes an LU factorization of a general M-by-N matrix A
117 : : // using partial pivoting with row interchanges.
118 : : void MOAB_dgetrf( int* M, int* N, double* A, int* lda, int* IPIV, int* INFO );
119 : :
120 : : // Computes the inverse of a matrix using the LU factorization
121 : : // computed by DGETRF.
122 : : void MOAB_dgetri( int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO );
123 : : }
124 : :
125 : : #include <cstring>
126 : : #define MOAB_DMEMZERO( a, b ) memset( a, 0, b * sizeof( double ) )
127 : :
128 : : #endif
129 : :
130 : : namespace moab
131 : : {
132 : :
133 : : namespace Matrix
134 : : {
135 : :
136 : : template < typename Matrix >
137 : 1424945 : inline Matrix mmult3( const Matrix& a, const Matrix& b )
138 : : {
139 : 1424945 : return Matrix( a( 0, 0 ) * b( 0, 0 ) + a( 0, 1 ) * b( 1, 0 ) + a( 0, 2 ) * b( 2, 0 ),
140 : 1424945 : a( 0, 0 ) * b( 0, 1 ) + a( 0, 1 ) * b( 1, 1 ) + a( 0, 2 ) * b( 2, 1 ),
141 : 1424945 : a( 0, 0 ) * b( 0, 2 ) + a( 0, 1 ) * b( 1, 2 ) + a( 0, 2 ) * b( 2, 2 ),
142 : 1424945 : a( 1, 0 ) * b( 0, 0 ) + a( 1, 1 ) * b( 1, 0 ) + a( 1, 2 ) * b( 2, 0 ),
143 : 1424945 : a( 1, 0 ) * b( 0, 1 ) + a( 1, 1 ) * b( 1, 1 ) + a( 1, 2 ) * b( 2, 1 ),
144 : 1424945 : a( 1, 0 ) * b( 0, 2 ) + a( 1, 1 ) * b( 1, 2 ) + a( 1, 2 ) * b( 2, 2 ),
145 : 1424945 : a( 2, 0 ) * b( 0, 0 ) + a( 2, 1 ) * b( 1, 0 ) + a( 2, 2 ) * b( 2, 0 ),
146 : 1424945 : a( 2, 0 ) * b( 0, 1 ) + a( 2, 1 ) * b( 1, 1 ) + a( 2, 2 ) * b( 2, 1 ),
147 : 2849890 : a( 2, 0 ) * b( 0, 2 ) + a( 2, 1 ) * b( 1, 2 ) + a( 2, 2 ) * b( 2, 2 ) );
148 : : }
149 : :
150 : : template < typename Matrix >
151 : : inline const Matrix inverse( const Matrix& d )
152 : : {
153 : : const double det = 1.0 / determinant3( d );
154 : : return inverse( d, det );
155 : : }
156 : :
157 : : template < typename Vector, typename Matrix >
158 : : inline Vector vector_matrix( const Vector& v, const Matrix& m )
159 : : {
160 : : return Vector( v[0] * m( 0, 0 ) + v[1] * m( 1, 0 ) + v[2] * m( 2, 0 ),
161 : : v[0] * m( 0, 1 ) + v[1] * m( 1, 1 ) + v[2] * m( 2, 1 ),
162 : : v[0] * m( 0, 2 ) + v[1] * m( 1, 2 ) + v[2] * m( 2, 2 ) );
163 : : }
164 : :
165 : : template < typename Vector, typename Matrix >
166 : 165857 : inline Vector matrix_vector( const Matrix& m, const Vector& v )
167 : : {
168 : 165857 : Vector res = v;
169 : 165857 : res[0] = v[0] * m( 0, 0 ) + v[1] * m( 0, 1 ) + v[2] * m( 0, 2 );
170 : 165857 : res[1] = v[0] * m( 1, 0 ) + v[1] * m( 1, 1 ) + v[2] * m( 1, 2 );
171 : 165857 : res[2] = v[0] * m( 2, 0 ) + v[1] * m( 2, 1 ) + v[2] * m( 2, 2 );
172 : 165857 : return res;
173 : : }
174 : :
175 : : } // namespace Matrix
176 : :
177 : : class Matrix3
178 : : {
179 : :
180 : : public:
181 : : const static int size = 9;
182 : :
183 : : private:
184 : : #ifndef MOAB_HAVE_LAPACK
185 : : Eigen::Matrix3d _mat;
186 : : #else
187 : : double _mat[size];
188 : : #endif
189 : :
190 : : public:
191 : : // Default Constructor
192 : 208151 : inline Matrix3()
193 : : {
194 : : #ifndef MOAB_HAVE_LAPACK
195 : : _mat.fill( 0.0 );
196 : : #else
197 : 208151 : MOAB_DMEMZERO( _mat, Matrix3::size );
198 : : #endif
199 : 208151 : }
200 : :
201 : : #ifndef MOAB_HAVE_LAPACK
202 : : inline Matrix3( Eigen::Matrix3d mat ) : _mat( mat ) {}
203 : : #endif
204 : :
205 : : // TODO: Deprecate this.
206 : : // Then we can go from three Constructors to one.
207 : 1478023 : inline Matrix3( double diagonal )
208 : : {
209 : : #ifndef MOAB_HAVE_LAPACK
210 : : _mat << diagonal, 0.0, 0.0, 0.0, diagonal, 0.0, 0.0, 0.0, diagonal;
211 : : #else
212 : 1478023 : MOAB_DMEMZERO( _mat, Matrix3::size );
213 : 1478023 : _mat[0] = _mat[4] = _mat[8] = diagonal;
214 : : #endif
215 : 1478023 : }
216 : :
217 : 9 : inline Matrix3( const CartVect& diagonal )
218 : : {
219 : : #ifndef MOAB_HAVE_LAPACK
220 : : _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
221 : : #else
222 : 9 : MOAB_DMEMZERO( _mat, Matrix3::size );
223 : 9 : _mat[0] = diagonal[0];
224 : 9 : _mat[4] = diagonal[1];
225 : 9 : _mat[8] = diagonal[2];
226 : : #endif
227 : 9 : }
228 : :
229 : : // TODO: not strictly correct as the Matrix3 object
230 : : // is a double d[ 9] so the only valid model of T is
231 : : // double, or any refinement (int, float)
232 : : //*but* it doesn't really matter anything else
233 : : // will fail to compile.
234 : : inline Matrix3( const std::vector< double >& diagonal )
235 : : {
236 : : #ifndef MOAB_HAVE_LAPACK
237 : : _mat << diagonal[0], 0.0, 0.0, 0.0, diagonal[1], 0.0, 0.0, 0.0, diagonal[2];
238 : : #else
239 : : MOAB_DMEMZERO( _mat, Matrix3::size );
240 : : _mat[0] = diagonal[0];
241 : : _mat[4] = diagonal[1];
242 : : _mat[8] = diagonal[2];
243 : : #endif
244 : : }
245 : :
246 : 4297547 : inline Matrix3( double v00, double v01, double v02, double v10, double v11, double v12, double v20, double v21,
247 : : double v22 )
248 : : {
249 : : #ifndef MOAB_HAVE_LAPACK
250 : : _mat << v00, v01, v02, v10, v11, v12, v20, v21, v22;
251 : : #else
252 : 4297547 : MOAB_DMEMZERO( _mat, Matrix3::size );
253 : 4297547 : _mat[0] = v00;
254 : 4297547 : _mat[1] = v01;
255 : 4297547 : _mat[2] = v02;
256 : 4297547 : _mat[3] = v10;
257 : 4297547 : _mat[4] = v11;
258 : 4297547 : _mat[5] = v12;
259 : 4297547 : _mat[6] = v20;
260 : 4297547 : _mat[7] = v21;
261 : 4297547 : _mat[8] = v22;
262 : : #endif
263 : 4297547 : }
264 : :
265 : : // Copy constructor
266 : 2162941 : Matrix3( const Matrix3& f )
267 : : {
268 : : #ifndef MOAB_HAVE_LAPACK
269 : : _mat = f._mat;
270 : : #else
271 : 2162941 : memcpy( _mat, f._mat, size * sizeof( double ) );
272 : : #endif
273 : 2162941 : }
274 : :
275 : : // Weird constructors
276 : : template < typename Vector >
277 : 5 : inline Matrix3( const Vector& row0, const Vector& row1, const Vector& row2, const bool isRow )
278 : : {
279 : : #ifndef MOAB_HAVE_LAPACK
280 : : if( isRow ) { _mat << row0[0], row0[1], row0[2], row1[0], row1[1], row1[2], row2[0], row2[1], row2[2]; }
281 : : else
282 : : {
283 : : _mat << row0[0], row1[0], row2[0], row0[1], row1[1], row2[1], row0[2], row1[2], row2[2];
284 : : }
285 : : #else
286 : 5 : MOAB_DMEMZERO( _mat, Matrix3::size );
287 [ + + ]: 5 : if( isRow )
288 : : {
289 : 2 : _mat[0] = row0[0];
290 : 2 : _mat[1] = row0[1];
291 : 2 : _mat[2] = row0[2];
292 : 2 : _mat[3] = row1[0];
293 : 2 : _mat[4] = row1[1];
294 : 2 : _mat[5] = row1[2];
295 : 2 : _mat[6] = row2[0];
296 : 2 : _mat[7] = row2[1];
297 : 2 : _mat[8] = row2[2];
298 : : }
299 : : else
300 : : {
301 : 3 : _mat[0] = row0[0];
302 : 3 : _mat[1] = row1[0];
303 : 3 : _mat[2] = row2[0];
304 : 3 : _mat[3] = row0[1];
305 : 3 : _mat[4] = row1[1];
306 : 3 : _mat[5] = row2[1];
307 : 3 : _mat[6] = row0[2];
308 : 3 : _mat[7] = row1[2];
309 : 3 : _mat[8] = row2[2];
310 : : }
311 : : #endif
312 : 5 : }
313 : :
314 : : #ifndef DEPRECATED
315 : : #ifdef __GNUC__
316 : : #define DEPRECATED __attribute__( ( deprecated ) )
317 : : #else
318 : : #pragma message( "WARNING: You need to implement DEPRECATED for this compiler" )
319 : : #define DEPRECATED
320 : : #endif
321 : : #endif
322 : :
323 : : /*
324 : : * \deprecated { Use instead the constructor with explicit fourth argument, bool isRow, above }
325 : : *
326 : : */
327 : 0 : inline Matrix3( const double v[size] )
328 : : {
329 : : #ifndef MOAB_HAVE_LAPACK
330 : : _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
331 : : #else
332 : 0 : memcpy( _mat, v, size * sizeof( double ) );
333 : : #endif
334 : 0 : }
335 : :
336 : 40 : inline void copyto( double v[Matrix3::size] )
337 : : {
338 : : #ifndef MOAB_HAVE_LAPACK
339 : : std::copy( _mat.data(), _mat.data() + size, v );
340 : : #else
341 : 40 : memcpy( v, _mat, size * sizeof( double ) );
342 : : #endif
343 : 40 : }
344 : :
345 : 34661 : inline Matrix3& operator=( const Matrix3& m )
346 : : {
347 : : #ifndef MOAB_HAVE_LAPACK
348 : : _mat = m._mat;
349 : : #else
350 : 34661 : memcpy( _mat, m._mat, size * sizeof( double ) );
351 : : #endif
352 : 34661 : return *this;
353 : : }
354 : :
355 : : inline Matrix3& operator=( const double v[size] )
356 : : {
357 : : #ifndef MOAB_HAVE_LAPACK
358 : : _mat << v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], v[8];
359 : : #else
360 : : memcpy( _mat, v, size * sizeof( double ) );
361 : : #endif
362 : : return *this;
363 : : }
364 : :
365 : 1170 : inline double* operator[]( unsigned i )
366 : : {
367 : : #ifndef MOAB_HAVE_LAPACK
368 : : return _mat.row( i ).data();
369 : : #else
370 : 1170 : return &_mat[i * 3]; // Row Major
371 : : #endif
372 : : }
373 : :
374 : 24201 : inline const double* operator[]( unsigned i ) const
375 : : {
376 : : #ifndef MOAB_HAVE_LAPACK
377 : : return _mat.row( i ).data();
378 : : #else
379 : 24201 : return &_mat[i * 3];
380 : : #endif
381 : : }
382 : :
383 : 1381484 : inline double& operator()( unsigned r, unsigned c )
384 : : {
385 : : #ifndef MOAB_HAVE_LAPACK
386 : : return _mat( r, c );
387 : : #else
388 : 1381484 : return _mat[r * 3 + c];
389 : : #endif
390 : : }
391 : :
392 : 78439798 : inline double operator()( unsigned r, unsigned c ) const
393 : : {
394 : : #ifndef MOAB_HAVE_LAPACK
395 : : return _mat( r, c );
396 : : #else
397 : 78439798 : return _mat[r * 3 + c];
398 : : #endif
399 : : }
400 : :
401 : 19236557 : inline double& operator()( unsigned i )
402 : : {
403 : : #ifndef MOAB_HAVE_LAPACK
404 : : return _mat( i );
405 : : #else
406 : 19236557 : return _mat[i];
407 : : #endif
408 : : }
409 : :
410 : : inline double operator()( unsigned i ) const
411 : : {
412 : : #ifndef MOAB_HAVE_LAPACK
413 : : return _mat( i );
414 : : #else
415 : : return _mat[i];
416 : : #endif
417 : : }
418 : :
419 : : // get pointer to array of nine doubles
420 : 27166 : inline double* array()
421 : : {
422 : : #ifndef MOAB_HAVE_LAPACK
423 : : return _mat.data();
424 : : #else
425 : 27166 : return _mat;
426 : : #endif
427 : : }
428 : :
429 : : inline const double* array() const
430 : : {
431 : : #ifndef MOAB_HAVE_LAPACK
432 : : return _mat.data();
433 : : #else
434 : : return _mat;
435 : : #endif
436 : : }
437 : :
438 : 714245 : inline Matrix3& operator+=( const Matrix3& m )
439 : : {
440 : : #ifndef MOAB_HAVE_LAPACK
441 : : _mat += m._mat;
442 : : #else
443 [ + + ]: 7142450 : for( int i = 0; i < Matrix3::size; ++i )
444 : 6428205 : _mat[i] += m._mat[i];
445 : : #endif
446 : 714245 : return *this;
447 : : }
448 : :
449 : 22701 : inline Matrix3& operator-=( const Matrix3& m )
450 : : {
451 : : #ifndef MOAB_HAVE_LAPACK
452 : : _mat -= m._mat;
453 : : #else
454 [ + + ]: 227010 : for( int i = 0; i < Matrix3::size; ++i )
455 : 204309 : _mat[i] -= m._mat[i];
456 : : #endif
457 : 22701 : return *this;
458 : : }
459 : :
460 : 12245 : inline Matrix3& operator*=( double s )
461 : : {
462 : : #ifndef MOAB_HAVE_LAPACK
463 : : _mat *= s;
464 : : #else
465 [ + + ]: 122450 : for( int i = 0; i < Matrix3::size; ++i )
466 : 110205 : _mat[i] *= s;
467 : : #endif
468 : 12245 : return *this;
469 : : }
470 : :
471 : 22702 : inline Matrix3& operator/=( double s )
472 : : {
473 : : #ifndef MOAB_HAVE_LAPACK
474 : : _mat /= s;
475 : : #else
476 [ + + ]: 227020 : for( int i = 0; i < Matrix3::size; ++i )
477 : 204318 : _mat[i] /= s;
478 : : #endif
479 : 22702 : return *this;
480 : : }
481 : :
482 : : inline Matrix3& operator*=( const Matrix3& m )
483 : : {
484 : : #ifndef MOAB_HAVE_LAPACK
485 : : _mat *= m._mat;
486 : : #else
487 : : // Uncomment below if you want point-wise multiplication instead (.*)
488 : : // for (int i=0; i < Matrix3::size; ++i) _mat[i] *= m._mat[i];
489 : : std::vector< double > dmat;
490 : : dmat.assign( _mat, _mat + size );
491 : : _mat[0] = dmat[0] * m._mat[0] + dmat[1] * m._mat[3] + dmat[2] * m._mat[6];
492 : : _mat[1] = dmat[0] * m._mat[1] + dmat[1] * m._mat[4] + dmat[2] * m._mat[7];
493 : : _mat[2] = dmat[0] * m._mat[2] + dmat[1] * m._mat[5] + dmat[2] * m._mat[8];
494 : : _mat[3] = dmat[3] * m._mat[0] + dmat[4] * m._mat[3] + dmat[5] * m._mat[6];
495 : : _mat[4] = dmat[3] * m._mat[1] + dmat[4] * m._mat[4] + dmat[5] * m._mat[7];
496 : : _mat[5] = dmat[3] * m._mat[2] + dmat[4] * m._mat[5] + dmat[5] * m._mat[8];
497 : : _mat[6] = dmat[6] * m._mat[0] + dmat[7] * m._mat[3] + dmat[8] * m._mat[6];
498 : : _mat[7] = dmat[6] * m._mat[1] + dmat[7] * m._mat[4] + dmat[8] * m._mat[7];
499 : : _mat[8] = dmat[6] * m._mat[2] + dmat[7] * m._mat[5] + dmat[8] * m._mat[8];
500 : : #endif
501 : : return *this;
502 : : }
503 : :
504 : 22703 : inline bool is_symmetric()
505 : : {
506 : 22703 : const double EPS = 1e-13;
507 : : #ifndef MOAB_HAVE_LAPACK
508 : : if( ( fabs( _mat( 1 ) - _mat( 3 ) ) < EPS ) && ( fabs( _mat( 2 ) - _mat( 6 ) ) < EPS ) &&
509 : : ( fabs( _mat( 5 ) - _mat( 7 ) ) < EPS ) )
510 : : return true;
511 : : #else
512 [ + - ][ + - ]: 22703 : if( ( fabs( _mat[1] - _mat[3] ) < EPS ) && ( fabs( _mat[2] - _mat[6] ) < EPS ) &&
[ + - ]
513 : 22703 : ( fabs( _mat[5] - _mat[7] ) < EPS ) )
514 : 22703 : return true;
515 : : #endif
516 : : else
517 : 0 : return false;
518 : : }
519 : :
520 : : inline bool is_positive_definite()
521 : : {
522 : : #ifndef MOAB_HAVE_LAPACK
523 : : double subdet6 = _mat( 1 ) * _mat( 5 ) - _mat( 2 ) * _mat( 4 );
524 : : double subdet7 = _mat( 2 ) * _mat( 3 ) - _mat( 0 ) * _mat( 5 );
525 : : double subdet8 = _mat( 0 ) * _mat( 4 ) - _mat( 1 ) * _mat( 3 );
526 : : // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
527 : : const double det = _mat( 6 ) * subdet6 + _mat( 7 ) * subdet7 + _mat( 8 ) * subdet8;
528 : : return _mat( 0 ) > 0 && subdet8 > 0 && det > 0;
529 : : #else
530 : : double subdet6 = _mat[1] * _mat[5] - _mat[2] * _mat[4];
531 : : double subdet7 = _mat[2] * _mat[3] - _mat[0] * _mat[5];
532 : : double subdet8 = _mat[0] * _mat[4] - _mat[1] * _mat[3];
533 : : // Determinant:= d(6)*subdet6 + d(7)*subdet7 + d(8)*subdet8;
534 : : const double det = _mat[6] * subdet6 + _mat[7] * subdet7 + _mat[8] * subdet8;
535 : : return _mat[0] > 0 && subdet8 > 0 && det > 0;
536 : : #endif
537 : : }
538 : :
539 : : template < typename Vector >
540 : 22703 : inline ErrorCode eigen_decomposition( Vector& evals, Matrix3& evecs )
541 : : {
542 [ + - ]: 22703 : const bool bisSymmetric = this->is_symmetric();
543 : : #ifndef MOAB_HAVE_LAPACK
544 : : if( bisSymmetric )
545 : : {
546 : : Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d > eigensolver( this->_mat );
547 : : if( eigensolver.info() != Eigen::Success ) return MB_FAILURE;
548 : : const Eigen::SelfAdjointEigenSolver< Eigen::Matrix3d >::RealVectorType& e3evals = eigensolver.eigenvalues();
549 : : evals[0] = e3evals( 0 );
550 : : evals[1] = e3evals( 1 );
551 : : evals[2] = e3evals( 2 );
552 : : evecs._mat = eigensolver.eigenvectors(); //.col(1)
553 : : return MB_SUCCESS;
554 : : }
555 : : else
556 : : {
557 : : MB_CHK_SET_ERR( MB_FAILURE, "Unsymmetric matrix implementation with Eigen3 is currently not provided." );
558 : : // Eigen::EigenSolver<Eigen::Matrix3d> eigensolver(this->_mat, true);
559 : : // if (eigensolver.info() != Eigen::Success)
560 : : // return MB_FAILURE;
561 : : // const Eigen::EigenSolver<Eigen::Matrix3d>::EigenvalueType& e3evals =
562 : : // eigensolver.eigenvalues().real(); evals[0] = e3evals(0); evals[1] = e3evals(1);
563 : : // evals[2] = e3evals(2); evecs._mat = eigensolver.eigenvectors().real(); //.col(1)
564 : : // return MB_SUCCESS;
565 : : }
566 : : #else
567 : : int info;
568 : : /* Solve eigenproblem */
569 : : double devreal[3], drevecs[9];
570 [ - + ]: 22703 : if( !bisSymmetric )
571 : : {
572 : : double devimag[3], dlevecs[9], dwork[102];
573 : 0 : char dgeev_opts[2] = { 'N', 'V' };
574 : 0 : int N = 3, LWORK = 102, NL = 1, NR = N;
575 [ # # ]: 0 : std::vector< double > devmat;
576 [ # # ]: 0 : devmat.assign( _mat, _mat + size );
577 [ # # ]: 0 : MOAB_dgeev( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, devimag, dlevecs, &NL, drevecs,
578 [ # # ]: 0 : &NR, dwork, &LWORK, &info );
579 : : // The result eigenvalues are ordered as high-->low
580 [ # # ]: 0 : evals[0] = devreal[2];
581 [ # # ]: 0 : evals[1] = devreal[1];
582 [ # # ]: 0 : evals[2] = devreal[0];
583 : 0 : evecs._mat[0] = drevecs[6];
584 : 0 : evecs._mat[1] = drevecs[3];
585 : 0 : evecs._mat[2] = drevecs[0];
586 : 0 : evecs._mat[3] = drevecs[7];
587 : 0 : evecs._mat[4] = drevecs[4];
588 : 0 : evecs._mat[5] = drevecs[1];
589 : 0 : evecs._mat[6] = drevecs[8];
590 : 0 : evecs._mat[7] = drevecs[5];
591 : 0 : evecs._mat[8] = drevecs[2];
592 [ # # ][ # # ]: 0 : std::cout << "DGEEV: Optimal work vector: dsize = " << dwork[0] << ".\n";
[ # # ]
593 : : }
594 : : else
595 : : {
596 : 22703 : char dgeev_opts[2] = { 'V', 'L' };
597 : 22703 : const bool find_optimal = false;
598 [ + - ]: 22703 : std::vector< int > iwork( 18 );
599 [ + - ]: 45406 : std::vector< double > devmat( 9, 0.0 );
600 [ + - ]: 45406 : std::vector< double > dwork( 38 );
601 : 22703 : int N = 3, lwork = 38, liwork = 18;
602 [ + - ]: 22703 : devmat[0] = _mat[0];
603 [ + - ]: 22703 : devmat[1] = _mat[1];
604 [ + - ]: 22703 : devmat[2] = _mat[2];
605 [ + - ]: 22703 : devmat[4] = _mat[4];
606 [ + - ]: 22703 : devmat[5] = _mat[5];
607 [ + - ]: 22703 : devmat[8] = _mat[8];
608 : : if( find_optimal )
609 : : {
610 : : int _lwork = -1;
611 : : int _liwork = -1;
612 : : double query_work_size = 0;
613 : : int query_iwork_size = 0;
614 : : // Make an empty call to find the optimal work vector size
615 : : MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, NULL, &N, NULL, &query_work_size, &_lwork,
616 : : &query_iwork_size, &_liwork, &info );
617 : : lwork = (int)query_work_size;
618 : : dwork.resize( lwork );
619 : : liwork = query_iwork_size;
620 : : iwork.resize( liwork );
621 : : std::cout << "DSYEVD: Optimal work vector: dsize = " << lwork << ", and isize = " << liwork << ".\n";
622 : : }
623 : :
624 [ + - ]: 22703 : MOAB_dsyevd( &dgeev_opts[0], &dgeev_opts[1], &N, &devmat[0], &N, devreal, &dwork[0], &lwork, &iwork[0],
625 [ + - ][ + - ]: 22703 : &liwork, &info );
[ + - ]
626 [ + + ]: 227030 : for( int i = 0; i < 9; ++i )
627 [ + - ]: 204327 : drevecs[i] = devmat[i];
628 : : // The result eigenvalues are ordered as low-->high, but vectors are in rows of A.
629 [ + - ]: 22703 : evals[0] = devreal[0];
630 [ + - ]: 22703 : evals[1] = devreal[1];
631 [ + - ]: 22703 : evals[2] = devreal[2];
632 : 22703 : evecs._mat[0] = drevecs[0];
633 : 22703 : evecs._mat[3] = drevecs[1];
634 : 22703 : evecs._mat[6] = drevecs[2];
635 : 22703 : evecs._mat[1] = drevecs[3];
636 : 22703 : evecs._mat[4] = drevecs[4];
637 : 22703 : evecs._mat[7] = drevecs[5];
638 : 22703 : evecs._mat[2] = drevecs[6];
639 : 22703 : evecs._mat[5] = drevecs[7];
640 : 22703 : evecs._mat[8] = drevecs[8];
641 : : }
642 : :
643 [ + - ]: 22703 : if( !info ) { return MB_SUCCESS; }
644 : : else
645 : : {
646 [ # # ][ # # ]: 0 : std::cout << "Failure in LAPACK_" << ( bisSymmetric ? "DSYEVD" : "DGEEV" )
[ # # ][ # # ]
647 : : << " call for eigen decomposition.\n";
648 [ # # ][ # # ]: 0 : std::cout << "Failed with error = " << info << ".\n";
[ # # ]
649 : 22703 : return MB_FAILURE;
650 : : }
651 : : #endif
652 : : }
653 : :
654 : : inline void transpose_inplace()
655 : : {
656 : : #ifndef MOAB_HAVE_LAPACK
657 : : _mat.transposeInPlace();
658 : : #else
659 : : Matrix3 mtmp( *this );
660 : : _mat[1] = mtmp._mat[3];
661 : : _mat[3] = mtmp._mat[1];
662 : : _mat[2] = mtmp._mat[6];
663 : : _mat[6] = mtmp._mat[2];
664 : : _mat[5] = mtmp._mat[7];
665 : : _mat[7] = mtmp._mat[5];
666 : : #endif
667 : : }
668 : :
669 : 19681 : inline Matrix3 transpose() const
670 : : {
671 : : #ifndef MOAB_HAVE_LAPACK
672 : : return Matrix3( _mat.transpose() );
673 : : #else
674 : 19681 : Matrix3 mtmp( *this );
675 : 19681 : mtmp._mat[1] = _mat[3];
676 : 19681 : mtmp._mat[3] = _mat[1];
677 : 19681 : mtmp._mat[2] = _mat[6];
678 : 19681 : mtmp._mat[6] = _mat[2];
679 : 19681 : mtmp._mat[5] = _mat[7];
680 : 19681 : mtmp._mat[7] = _mat[5];
681 : 19681 : return mtmp;
682 : : #endif
683 : : }
684 : :
685 : : template < typename Vector >
686 : : inline void copycol( int index, Vector& vol )
687 : : {
688 : : #ifndef MOAB_HAVE_LAPACK
689 : : _mat.col( index ).swap( vol );
690 : : #else
691 : : switch( index )
692 : : {
693 : : case 0:
694 : : _mat[0] = vol[0];
695 : : _mat[3] = vol[1];
696 : : _mat[6] = vol[2];
697 : : break;
698 : : case 1:
699 : : _mat[1] = vol[0];
700 : : _mat[4] = vol[1];
701 : : _mat[7] = vol[2];
702 : : break;
703 : : case 2:
704 : : _mat[2] = vol[0];
705 : : _mat[5] = vol[1];
706 : : _mat[8] = vol[2];
707 : : break;
708 : : }
709 : : #endif
710 : : }
711 : :
712 : 1047 : inline void swapcol( int srcindex, int destindex )
713 : : {
714 [ - + ]: 1047 : assert( srcindex < Matrix3::size );
715 [ - + ]: 1047 : assert( destindex < Matrix3::size );
716 : : #ifndef MOAB_HAVE_LAPACK
717 : : _mat.col( srcindex ).swap( _mat.col( destindex ) );
718 : : #else
719 [ + - ]: 1047 : CartVect svol = this->vcol< CartVect >( srcindex );
720 [ + - ]: 1047 : CartVect dvol = this->vcol< CartVect >( destindex );
721 [ + + - - ]: 1047 : switch( srcindex )
722 : : {
723 : : case 0:
724 [ + - ]: 116 : _mat[0] = dvol[0];
725 [ + - ]: 116 : _mat[3] = dvol[1];
726 [ + - ]: 116 : _mat[6] = dvol[2];
727 : 116 : break;
728 : : case 1:
729 [ + - ]: 931 : _mat[1] = dvol[0];
730 [ + - ]: 931 : _mat[4] = dvol[1];
731 [ + - ]: 931 : _mat[7] = dvol[2];
732 : 931 : break;
733 : : case 2:
734 [ # # ]: 0 : _mat[2] = dvol[0];
735 [ # # ]: 0 : _mat[5] = dvol[1];
736 [ # # ]: 0 : _mat[8] = dvol[2];
737 : 0 : break;
738 : : }
739 [ - + + - ]: 1047 : switch( destindex )
740 : : {
741 : : case 0:
742 [ # # ]: 0 : _mat[0] = svol[0];
743 [ # # ]: 0 : _mat[3] = svol[1];
744 [ # # ]: 0 : _mat[6] = svol[2];
745 : 0 : break;
746 : : case 1:
747 [ + - ]: 80 : _mat[1] = svol[0];
748 [ + - ]: 80 : _mat[4] = svol[1];
749 [ + - ]: 80 : _mat[7] = svol[2];
750 : 80 : break;
751 : : case 2:
752 [ + - ]: 967 : _mat[2] = svol[0];
753 [ + - ]: 967 : _mat[5] = svol[1];
754 [ + - ]: 967 : _mat[8] = svol[2];
755 : 967 : break;
756 : : }
757 : : #endif
758 : 1047 : }
759 : :
760 : : template < typename Vector >
761 : 2094 : inline Vector vcol( int index ) const
762 : : {
763 [ - + ]: 2094 : assert( index < Matrix3::size );
764 : : #ifndef MOAB_HAVE_LAPACK
765 : : return _mat.col( index );
766 : : #else
767 [ + + + - ]: 2094 : switch( index )
768 : : {
769 : : case 0:
770 : 116 : return Vector( _mat[0], _mat[3], _mat[6] );
771 : : case 1:
772 : 1011 : return Vector( _mat[1], _mat[4], _mat[7] );
773 : : case 2:
774 : 967 : return Vector( _mat[2], _mat[5], _mat[8] );
775 : : }
776 : 0 : return Vector( 0.0 );
777 : : #endif
778 : : }
779 : :
780 : 15 : inline void colscale( int index, double scale )
781 : : {
782 [ - + ]: 15 : assert( index < Matrix3::size );
783 : : #ifndef MOAB_HAVE_LAPACK
784 : : _mat.col( index ) *= scale;
785 : : #else
786 [ + + + - ]: 15 : switch( index )
787 : : {
788 : : case 0:
789 : 5 : _mat[0] *= scale;
790 : 5 : _mat[3] *= scale;
791 : 5 : _mat[6] *= scale;
792 : 5 : break;
793 : : case 1:
794 : 5 : _mat[1] *= scale;
795 : 5 : _mat[4] *= scale;
796 : 5 : _mat[7] *= scale;
797 : 5 : break;
798 : : case 2:
799 : 5 : _mat[2] *= scale;
800 : 5 : _mat[5] *= scale;
801 : 5 : _mat[8] *= scale;
802 : 5 : break;
803 : : }
804 : : #endif
805 : 15 : }
806 : :
807 : : inline void rowscale( int index, double scale )
808 : : {
809 : : assert( index < Matrix3::size );
810 : : #ifndef MOAB_HAVE_LAPACK
811 : : _mat.row( index ) *= scale;
812 : : #else
813 : : switch( index )
814 : : {
815 : : case 0:
816 : : _mat[0] *= scale;
817 : : _mat[1] *= scale;
818 : : _mat[2] *= scale;
819 : : break;
820 : : case 1:
821 : : _mat[3] *= scale;
822 : : _mat[4] *= scale;
823 : : _mat[5] *= scale;
824 : : break;
825 : : case 2:
826 : : _mat[6] *= scale;
827 : : _mat[7] *= scale;
828 : : _mat[8] *= scale;
829 : : break;
830 : : }
831 : : #endif
832 : : }
833 : :
834 : 3566546 : inline CartVect col( int index ) const
835 : : {
836 [ - + ]: 3566546 : assert( index < Matrix3::size );
837 : : #ifndef MOAB_HAVE_LAPACK
838 : : Eigen::Vector3d mvec = _mat.col( index );
839 : : return CartVect( mvec[0], mvec[1], mvec[2] );
840 : : #else
841 [ + + + - ]: 3566546 : switch( index )
842 : : {
843 : : case 0:
844 : 949466 : return CartVect( _mat[0], _mat[3], _mat[6] );
845 : : case 1:
846 : 1080033 : return CartVect( _mat[1], _mat[4], _mat[7] );
847 : : case 2:
848 : 1537047 : return CartVect( _mat[2], _mat[5], _mat[8] );
849 : : }
850 : 0 : return CartVect( 0.0 );
851 : : #endif
852 : : }
853 : :
854 : : inline CartVect row( int index ) const
855 : : {
856 : : assert( index < Matrix3::size );
857 : : #ifndef MOAB_HAVE_LAPACK
858 : : Eigen::Vector3d mvec = _mat.row( index );
859 : : return CartVect( mvec[0], mvec[1], mvec[2] );
860 : : #else
861 : : switch( index )
862 : : {
863 : : case 0:
864 : : return CartVect( _mat[0], _mat[1], _mat[2] );
865 : : case 1:
866 : : return CartVect( _mat[3], _mat[4], _mat[5] );
867 : : case 2:
868 : : return CartVect( _mat[6], _mat[7], _mat[8] );
869 : : }
870 : : return CartVect( 0.0 );
871 : : #endif
872 : : }
873 : :
874 : : friend Matrix3 operator+( const Matrix3& a, const Matrix3& b );
875 : : friend Matrix3 operator-( const Matrix3& a, const Matrix3& b );
876 : : friend Matrix3 operator*( const Matrix3& a, const Matrix3& b );
877 : :
878 : 46018 : inline double determinant() const
879 : : {
880 : : #ifndef MOAB_HAVE_LAPACK
881 : : return _mat.determinant();
882 : : #else
883 : 92036 : return ( _mat[0] * _mat[4] * _mat[8] + _mat[1] * _mat[5] * _mat[6] + _mat[2] * _mat[3] * _mat[7] -
884 : 46018 : _mat[0] * _mat[5] * _mat[7] - _mat[1] * _mat[3] * _mat[8] - _mat[2] * _mat[4] * _mat[6] );
885 : : #endif
886 : : }
887 : :
888 : 17642 : inline Matrix3 inverse() const
889 : : {
890 : : #ifndef MOAB_HAVE_LAPACK
891 : : return Matrix3( _mat.inverse() );
892 : : #else
893 : : // return Matrix::compute_inverse( *this, this->determinant() );
894 : 17642 : Matrix3 m( 0.0 );
895 : 17642 : const double d_determinant = 1.0 / this->determinant();
896 : 17642 : m._mat[0] = d_determinant * ( _mat[4] * _mat[8] - _mat[5] * _mat[7] );
897 : 17642 : m._mat[1] = d_determinant * ( _mat[2] * _mat[7] - _mat[8] * _mat[1] );
898 : 17642 : m._mat[2] = d_determinant * ( _mat[1] * _mat[5] - _mat[4] * _mat[2] );
899 : 17642 : m._mat[3] = d_determinant * ( _mat[5] * _mat[6] - _mat[8] * _mat[3] );
900 : 17642 : m._mat[4] = d_determinant * ( _mat[0] * _mat[8] - _mat[6] * _mat[2] );
901 : 17642 : m._mat[5] = d_determinant * ( _mat[2] * _mat[3] - _mat[5] * _mat[0] );
902 : 17642 : m._mat[6] = d_determinant * ( _mat[3] * _mat[7] - _mat[6] * _mat[4] );
903 : 17642 : m._mat[7] = d_determinant * ( _mat[1] * _mat[6] - _mat[7] * _mat[0] );
904 : 17642 : m._mat[8] = d_determinant * ( _mat[0] * _mat[4] - _mat[3] * _mat[1] );
905 : 17642 : return m;
906 : : #endif
907 : : }
908 : :
909 : : inline bool invert()
910 : : {
911 : : bool invertible = false;
912 : : double d_determinant;
913 : : #ifndef MOAB_HAVE_LAPACK
914 : : Eigen::Matrix3d invMat;
915 : : _mat.computeInverseAndDetWithCheck( invMat, d_determinant, invertible );
916 : : if( !Util::is_finite( d_determinant ) ) return false;
917 : : _mat = invMat;
918 : : return invertible;
919 : : #else
920 : : d_determinant = this->determinant();
921 : : if( d_determinant > 1e-13 ) invertible = true;
922 : : d_determinant = 1.0 / d_determinant; // invert the determinant
923 : : std::vector< double > _m;
924 : : _m.assign( _mat, _mat + size );
925 : : _mat[0] = d_determinant * ( _m[4] * _m[8] - _m[5] * _m[7] );
926 : : _mat[1] = d_determinant * ( _m[2] * _m[7] - _m[8] * _m[1] );
927 : : _mat[2] = d_determinant * ( _m[1] * _m[5] - _m[4] * _m[2] );
928 : : _mat[3] = d_determinant * ( _m[5] * _m[6] - _m[8] * _m[3] );
929 : : _mat[4] = d_determinant * ( _m[0] * _m[8] - _m[6] * _m[2] );
930 : : _mat[5] = d_determinant * ( _m[2] * _m[3] - _m[5] * _m[0] );
931 : : _mat[6] = d_determinant * ( _m[3] * _m[7] - _m[6] * _m[4] );
932 : : _mat[7] = d_determinant * ( _m[1] * _m[6] - _m[7] * _m[0] );
933 : : _mat[8] = d_determinant * ( _m[0] * _m[4] - _m[3] * _m[1] );
934 : : #endif
935 : : return invertible;
936 : : }
937 : :
938 : : // Calculate determinant of 2x2 submatrix composed of the
939 : : // elements not in the passed row or column.
940 : : inline double subdet( int r, int c ) const
941 : : {
942 : : assert( r >= 0 && c >= 0 );
943 : : if( r < 0 || c < 0 ) return DBL_MAX;
944 : : #ifndef MOAB_HAVE_LAPACK
945 : : const int r1 = ( r + 1 ) % 3, r2 = ( r + 2 ) % 3;
946 : : const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
947 : : return _mat( r1, c1 ) * _mat( r2, c2 ) - _mat( r1, c2 ) * _mat( r2, c1 );
948 : : #else
949 : : const int r1 = Matrix3::size * ( ( r + 1 ) % 3 ), r2 = Matrix3::size * ( ( r + 2 ) % 3 );
950 : : const int c1 = ( c + 1 ) % 3, c2 = ( c + 2 ) % 3;
951 : : return _mat[r1 + c1] * _mat[r2 + c2] - _mat[r1 + c2] * _mat[r2 + c1];
952 : : #endif
953 : : }
954 : :
955 : : inline void print( std::ostream& s ) const
956 : : {
957 : : #ifndef MOAB_HAVE_LAPACK
958 : : s << "| " << _mat( 0 ) << " " << _mat( 1 ) << " " << _mat( 2 ) << " | " << _mat( 3 ) << " " << _mat( 4 ) << " "
959 : : << _mat( 5 ) << " | " << _mat( 6 ) << " " << _mat( 7 ) << " " << _mat( 8 ) << " |";
960 : : #else
961 : : s << "| " << _mat[0] << " " << _mat[1] << " " << _mat[2] << " | " << _mat[3] << " " << _mat[4] << " " << _mat[5]
962 : : << " | " << _mat[6] << " " << _mat[7] << " " << _mat[8] << " |";
963 : : #endif
964 : : }
965 : :
966 : : }; // class Matrix3
967 : :
968 : : template < typename Vector >
969 : 2872564 : inline Matrix3 outer_product( const Vector& u, const Vector& v )
970 : : {
971 : 2872564 : 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],
972 : 5745128 : u[2] * v[1], u[2] * v[2] );
973 : : }
974 : :
975 : 2137394 : inline Matrix3 operator+( const Matrix3& a, const Matrix3& b )
976 : : {
977 : : #ifndef MOAB_HAVE_LAPACK
978 : : return Matrix3( a._mat + b._mat );
979 : : #else
980 : 2137394 : Matrix3 s( a );
981 [ + + ]: 21373940 : for( int i = 0; i < Matrix3::size; ++i )
982 : 19236546 : s( i ) += b._mat[i];
983 : 2137394 : return s;
984 : : #endif
985 : : }
986 : :
987 : : inline Matrix3 operator-( const Matrix3& a, const Matrix3& b )
988 : : {
989 : : #ifndef MOAB_HAVE_LAPACK
990 : : return Matrix3( a._mat - b._mat );
991 : : #else
992 : : Matrix3 s( a );
993 : : for( int i = 0; i < Matrix3::size; ++i )
994 : : s( i ) -= b._mat[i];
995 : : return s;
996 : : #endif
997 : : }
998 : :
999 : 1424945 : inline Matrix3 operator*( const Matrix3& a, const Matrix3& b )
1000 : : {
1001 : : #ifndef MOAB_HAVE_LAPACK
1002 : : return Matrix3( a._mat * b._mat );
1003 : : #else
1004 : 1424945 : return Matrix::mmult3( a, b );
1005 : : #endif
1006 : : }
1007 : :
1008 : : template < typename T >
1009 : : inline std::vector< T > operator*( const Matrix3& m, const std::vector< T >& v )
1010 : : {
1011 : : return moab::Matrix::matrix_vector( m, v );
1012 : : }
1013 : :
1014 : : template < typename T >
1015 : : inline std::vector< T > operator*( const std::vector< T >& v, const Matrix3& m )
1016 : : {
1017 : : return moab::Matrix::vector_matrix( v, m );
1018 : : }
1019 : :
1020 : 165854 : inline CartVect operator*( const Matrix3& m, const CartVect& v )
1021 : : {
1022 : 165854 : return moab::Matrix::matrix_vector( m, v );
1023 : : }
1024 : :
1025 : : inline CartVect operator*( const CartVect& v, const Matrix3& m )
1026 : : {
1027 : : return moab::Matrix::vector_matrix( v, m );
1028 : : }
1029 : :
1030 : : } // namespace moab
1031 : :
1032 : : #ifdef MOAB_HAVE_LAPACK
1033 : : #undef MOAB_DMEMZERO
1034 : : #endif
1035 : :
1036 : : #ifndef MOAB_MATRIX3_OPERATORLESS
1037 : : #define MOAB_MATRIX3_OPERATORLESS
1038 : 0 : inline std::ostream& operator<<( std::ostream& s, const moab::Matrix3& m )
1039 : : {
1040 : 0 : return s << "| " << m( 0, 0 ) << " " << m( 0, 1 ) << " " << m( 0, 2 ) << " | " << m( 1, 0 ) << " " << m( 1, 1 )
1041 : 0 : << " " << m( 1, 2 ) << " | " << m( 2, 0 ) << " " << m( 2, 1 ) << " " << m( 2, 2 ) << " |";
1042 : : }
1043 : : #endif // MOAB_MATRIX3_OPERATORLESS
1044 : : #endif // MOAB_MATRIX3_HPP
|
