From 475b73d79b204270181c9ce0e5ddb9a24ed60aa3 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 11 Jun 2018 14:38:54 +0200
Subject: [PATCH] Few clean-up in TinyMatrix
Now det for arbitrary dimension is only available to floating point
types.
---
src/algebra/TinyMatrix.hpp | 96 +++++++++++++++++++++-----------------
1 file changed, 53 insertions(+), 43 deletions(-)
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 577bd4fa0..9a9e9aa22 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -52,11 +52,11 @@ public:
TinyMatrix AB;
for (size_t i=0; i<N; ++i) {
for (size_t j=0; j<N; ++j) {
- T sum = A(i,0)*B(0,j);
- for (size_t k=1; k<N; ++k) {
- sum += A(i,k)*B(k,j);
- }
- AB(i,j) = sum;
+ T sum = A(i,0)*B(0,j);
+ for (size_t k=1; k<N; ++k) {
+ sum += A(i,k)*B(k,j);
+ }
+ AB(i,j) = sum;
}
}
return std::move(AB);
@@ -71,7 +71,7 @@ public:
for (size_t i=0; i<N; ++i) {
T sum = A(i,0)*x[0];
for (size_t j=1; j<N; ++j) {
- sum += A(i,j)*x[j];
+ sum += A(i,j)*x[j];
}
Ax[i] = sum;
}
@@ -83,22 +83,37 @@ public:
friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
{
if constexpr(N==1) {
- os << A(0,0);
+ os << A(0,0);
}
else {
os << '[';
for (size_t i=0; i<N; ++i) {
- os << '(' << A(i,0);
- for (size_t j=1; j<N; ++j) {
- os << ',' << A(i,j);
- }
- os << ')';
+ os << '(' << A(i,0);
+ for (size_t j=1; j<N; ++j) {
+ os << ',' << A(i,j);
+ }
+ os << ')';
}
os << ']';
}
return os;
}
+ KOKKOS_INLINE_FUNCTION
+ bool operator==(const TinyMatrix& A) const
+ {
+ for (size_t i=0; i<N*N; ++i) {
+ if (m_values[i] != A.m_values[i]) return false;
+ }
+ return true;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ bool operator!=(const TinyMatrix& A) const
+ {
+ return not this->operator==(A);
+ }
+
KOKKOS_INLINE_FUNCTION
TinyMatrix operator+(const TinyMatrix& A) const
{
@@ -167,14 +182,14 @@ public:
static_assert(std::is_arithmetic<T>(),"Cannot assign 'identity' value for non-arithmetic types");
for (size_t i=0; i<N; ++i) {
for (size_t j=0; j<N; ++j) {
- m_values[_index(i,j)] = (i==j) ? 1 : 0;
+ m_values[_index(i,j)] = (i==j) ? 1 : 0;
}
}
return *this;
}
KOKKOS_INLINE_FUNCTION
- const TinyMatrix& operator=(const TinyMatrix& A)
+ TinyMatrix& operator=(const TinyMatrix& A)
{
for (size_t i=0; i<N*N; ++i) {
m_values[i] = A.m_values[i];
@@ -194,10 +209,7 @@ public:
}
KOKKOS_INLINE_FUNCTION
- TinyMatrix()
- {
- ;
- }
+ TinyMatrix()=default;
KOKKOS_INLINE_FUNCTION
TinyMatrix(const ZeroType& z)
@@ -214,7 +226,7 @@ public:
static_assert(std::is_arithmetic<T>(),"Cannot construct from 'identity' value for non-arithmetic types");
for (size_t i=0; i<N; ++i) {
for (size_t j=0; j<N; ++j) {
- m_values[_index(i,j)] = (i==j) ? 1 : 0;
+ m_values[_index(i,j)] = (i==j) ? 1 : 0;
}
}
}
@@ -231,16 +243,13 @@ public:
TinyMatrix(TinyMatrix&& A) = default;
KOKKOS_INLINE_FUNCTION
- ~TinyMatrix()
- {
- ;
- }
+ ~TinyMatrix()=default;
};
template <size_t N, typename T>
KOKKOS_INLINE_FUNCTION
TinyMatrix<N,T> tensorProduct(const TinyVector<N,T>& x,
- const TinyVector<N,T>& y)
+ const TinyVector<N,T>& y)
{
TinyMatrix<N,T> A;
for (size_t i=0; i<N; ++i) {
@@ -256,39 +265,40 @@ KOKKOS_INLINE_FUNCTION
T det(const TinyMatrix<N,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+ static_assert(std::is_floating_point<T>::value, "determinent for arbitrary dimension N is defined for floating point types only");
TinyMatrix<N,T> M = A;
- TinyVector<N, int> index;
- for (int i=0; i<N; ++i) index[i]=i;
+ TinyVector<N, size_t> index;
+ for (size_t i=0; i<N; ++i) index[i]=i;
T determinent = 1;
- for (int i=0; i<N; ++i) {
- for (int j=i; j<N; ++j) {
- int l = j;
- const int J = index[j];
- for (int k=j+1; k<N; ++k) {
- if (std::abs(M(index[k],i)) > std::abs(M(J,i))) {
- l=k;
- }
+ for (size_t i=0; i<N; ++i) {
+ for (size_t j=i; j<N; ++j) {
+ size_t l = j;
+ const size_t J = index[j];
+ for (size_t k=j+1; k<N; ++k) {
+ if (std::abs(M(index[k],i)) > std::abs(M(J,i))) {
+ l=k;
+ }
}
if (l != j) {
- std::swap(index[l], index[j]);
- determinent *= -1;
+ std::swap(index[l], index[j]);
+ determinent *= -1;
}
}
- const int I = index[i];
+ const size_t I = index[i];
if (M(I,i)==0) return 0;
- double inv_Mii = 1./M(I,i);
- for (int k=i+1; k<N; ++k) {
- const int K = index[k];
+ const T inv_Mii = 1./M(I,i);
+ for (size_t k=i+1; k<N; ++k) {
+ const size_t K = index[k];
const T factor = M(K,i)*inv_Mii;
- for (int l=i+1; l<N; ++l) {
- M(K,l) -= factor*M(I,l);
+ for (size_t l=i+1; l<N; ++l) {
+ M(K,l) -= factor*M(I,l);
}
}
}
- for (int i=0; i<N; ++i) {
+ for (size_t i=0; i<N; ++i) {
determinent *= M(index[i],i);
}
return determinent;
--
GitLab