1 // Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
2 // other Axom Project Developers. See the top-level LICENSE file for details.
3 //
4 // SPDX-License-Identifier: (BSD-3-Clause)
5
6 #ifndef AXOM_NUMERICS_DETERMINANTS_HPP_
7 #define AXOM_NUMERICS_DETERMINANTS_HPP_
8
9 #include "axom/core/numerics/LU.hpp" // for lu_decompose()
10 #include "axom/core/numerics/Matrix.hpp" // for Matrix
11
12 #include "axom/core/Macros.hpp"
13
14 namespace axom
15 {
16 namespace numerics
17 {
18 // clang-format off
19 /// \name Matrix Operators
20 /// @{
21
22 /*!
23 * \brief Computes a 2X2 determinant of the given matrix.
24 * \accelerated
25 * \param [in] a00 matrix element
26 * \param [in] a01 matrix element
27 * \param [in] a10 matrix element
28 * \param [in] a11 matrix element
29 * \return det the determinant of the 2X2 matrix
30 */
31 template <typename real>
32 inline AXOM_HOST_DEVICE
determinant(const real & a00,const real & a01,const real & a10,const real & a11)33 real determinant(const real& a00, const real& a01,
34 const real& a10, const real& a11)
35 {
36 const real det = a00 * a11 - a10 * a01;
37 return det;
38 }
39
40 /*!
41 * \brief Computes the 3x3 determinant of the given matrix.
42 * \param [in] a00 matrix element
43 * \param [in] a01 matrix element
44 * \param [in] a02 matrix element
45 * \param [in] a10 matrix element
46 * \param [in] a11 matrix element
47 * \param [in] a12 matrix element
48 * \param [in] a20 matrix element
49 * \param [in] a21 matrix element
50 * \param [in] a22 matrix element
51 * \return det the determinant of the 3x3 matrix.
52 */
53 template <typename real>
determinant(const real & a00,const real & a01,const real & a02,const real & a10,const real & a11,const real & a12,const real & a20,const real & a21,const real & a22)54 inline real determinant(const real& a00, const real& a01, const real& a02,
55 const real& a10, const real& a11, const real& a12,
56 const real& a20, const real& a21, const real& a22)
57 {
58 const real m01 = a00 * a11 - a10 * a01;
59 const real m02 = a00 * a21 - a20 * a01;
60 const real m12 = a10 * a21 - a20 * a11;
61
62 const real det = m01 * a22 - m02 * a12 + m12 * a02;
63 return det;
64 }
65
66 /*!
67 * \brief Computes the 4x4 determinant of the given matrix.
68 *
69 * \accelerated
70 * \param [in] a00 matrix element
71 * \param [in] a01 matrix element
72 * \param [in] a02 matrix element
73 * \param [in] a03 matrix element
74 * \param [in] a10 matrix element
75 * \param [in] a11 matrix element
76 * \param [in] a12 matrix element
77 * \param [in] a13 matrix element
78 * \param [in] a20 matrix element
79 * \param [in] a21 matrix element
80 * \param [in] a22 matrix element
81 * \param [in] a23 matrix element
82 * \param [in] a30 matrix element
83 * \param [in] a31 matrix element
84 * \param [in] a32 matrix element
85 * \param [in] a33 matrix element
86 * \return det the determinant of the 4x4 matrix
87 */
88 template <typename real>
89 AXOM_HOST_DEVICE
determinant(const real & a00,const real & a01,const real & a02,const real & a03,const real & a10,const real & a11,const real & a12,const real & a13,const real & a20,const real & a21,const real & a22,const real & a23,const real & a30,const real & a31,const real & a32,const real & a33)90 inline real determinant(
91 const real& a00, const real& a01, const real& a02, const real& a03,
92 const real& a10, const real& a11, const real& a12, const real& a13,
93 const real& a20, const real& a21, const real& a22, const real& a23,
94 const real& a30, const real& a31, const real& a32, const real& a33)
95 {
96 const real m01 = a10 * a01 - a00 * a11;
97 const real m02 = a20 * a01 - a00 * a21;
98 const real m03 = a30 * a01 - a00 * a31;
99 const real m12 = a20 * a11 - a10 * a21;
100 const real m13 = a30 * a11 - a10 * a31;
101 const real m23 = a30 * a21 - a20 * a31;
102
103 const real m012 = m12 * a02 - m02 * a12 + m01 * a22;
104 const real m013 = m13 * a02 - m03 * a12 + m01 * a32;
105 const real m023 = m23 * a02 - m03 * a22 + m02 * a32;
106 const real m123 = m23 * a12 - m13 * a22 + m12 * a32;
107
108 const real det = m123 * a03 - m023 * a13 + m013 * a23 - m012 * a33;
109 return det;
110 }
111
112 /*!
113 * \brief Computes the determinant of the given square matrix.
114 * \param [in] A an \f$ N \times N \f$ input matrix
115 * \return det the computed determinant.
116 * \note if \f$ A \f$ is not square or empty, this function will return 0.0
117 */
118 template <typename real>
determinant(const Matrix<real> & A)119 real determinant(const Matrix<real>& A)
120 {
121 real det = 0.0;
122
123 if(!A.isSquare() || A.empty())
124 {
125 /* short-circuit */
126 return det;
127 }
128
129 const int N = A.getNumColumns();
130 if(N == 1)
131 {
132 det = A(0, 0);
133 }
134 else if(N == 2)
135 {
136
137 det = determinant(A(0,0), A(0,1),
138 A(1,0), A(1,1));
139
140 }
141 else if(N == 3)
142 {
143
144 det = determinant(A(0,0), A(0,1), A(0,2),
145 A(1,0), A(1,1), A(1,2),
146 A(2,0), A(2,1), A(2,2));
147
148 }
149 else if(N == 4)
150 {
151
152 det = determinant(A(0,0), A(0,1), A(0,2), A(0,3),
153 A(1,0), A(1,1), A(1,2), A(1,3),
154 A(2,0), A(2,1), A(2,2), A(2,3),
155 A(3,0), A(3,1), A(3,2), A(3,3));
156
157 }
158 else
159 {
160 Matrix<real> lu = A;
161 int* pivots = new int[N];
162
163 int rc = lu_decompose(lu, pivots);
164 if(rc == LU_SUCCESS)
165 {
166 // count number of row interchanges
167 int row_interchanges = 0;
168 for(int i = 0; i < N; ++i)
169 {
170 row_interchanges += (pivots[i] != i) ? 1 : 0;
171 } // END for all rows
172
173 det = ((row_interchanges & 1) == 0) ? 1.0 : -1.0;
174 for(int i = 0; i < N; ++i)
175 {
176 det *= lu(i, i);
177 }
178 }
179
180 delete[] pivots;
181 }
182
183 return (det);
184 }
185 /// @}
186 // clang-format on
187
188 } /* namespace numerics */
189 } /* namespace axom */
190
191 #endif
192