Merge pull request #1068 from su2code/feature_CGeneralSquareMatrixCM

pcarruscag · web-flow · commit 7cda79e44c64 · 2020-12-12T12:38:56.000Z
Added CSquareMatrixCM to the toolboxes
diff --git a/Common/include/blas_structure.hpp b/Common/include/blas_structure.hpp
@@ -94,6 +94,55 @@ class CBlasStructure {
   void axpy(const int n,    const su2double a,  const su2double *x,
             const int incx, su2double *y,       const int incy);
 
+  /*!
+   * \brief Invert a square matrix.
+   * \param[in] M - Size.
+   * \param[in,out] mat - Matrix, and inverse on exit.
+   */
+  template<class Mat>
+  static void inverse(const int M, Mat& mat) {
+    using Scalar = typename Mat::Scalar;
+
+    /*--- Copy the data from A into the augmented matrix and initialize mat with the identity. ---*/
+    Mat aug = mat;
+    mat = Scalar(0);
+    for(int j=0; j<M; ++j) mat(j,j) = 1;
+
+    /*--- Outer loop of the Gauss-Jordan elimination. ---*/
+    for(int j=0; j<M; ++j) {
+
+      /*--- Find the pivot in the current column. ---*/
+      int jj = j;
+      Scalar valMax = fabs(aug(j,j));
+      for(int i=j+1; i<M; ++i) {
+        Scalar val = fabs(aug(i,j));
+        if(val > valMax){
+          jj = i;
+          valMax = val;
+        }
+      }
+
+      /*--- Swap the rows j and jj, if needed. ---*/
+      if(jj > j) {
+        for(int k=j; k<M; ++k) std::swap(aug(j,k), aug(jj,k));
+        for(int k=0; k<M; ++k) std::swap(mat(j,k), mat(jj,k));
+      }
+
+      /*--- Performing row operations to form required identity
+            matrix out of the input matrix.  ---*/
+      for(int i=0; i<M; ++i) {
+        if(i == j) continue;
+        valMax = aug(i,j)/aug(j,j);
+        for(int k=j; k<M; ++k) aug(i,k) -= valMax*aug(j,k);
+        for(int k=0; k<M; ++k) mat(i,k) -= valMax*mat(j,k);
+      }
+
+      valMax = 1.0/aug(j,j);
+      for(int k=j; k<M; ++k) aug(j,k) *= valMax;
+      for(int k=0; k<M; ++k) mat(j,k) *= valMax;
+    }
+  }
+
 private:
 
 #if !(defined(HAVE_LIBXSMM) || defined(HAVE_BLAS) || defined(HAVE_MKL)) || (defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE))
diff --git a/Common/include/containers/C2DContainer.hpp b/Common/include/containers/C2DContainer.hpp
@@ -613,7 +613,8 @@ class C2DContainer :
  * \brief Useful typedefs with default template parameters
  */
 template<class T> using su2vector = C2DContainer<unsigned long, T, StorageType::ColumnMajor, 64, DynamicSize, 1>;
-template<class T> using su2matrix = C2DContainer<unsigned long, T, StorageType::RowMajor, 64, DynamicSize, DynamicSize>;
+template<class T> using su2matrix = C2DContainer<unsigned long, T, StorageType::RowMajor,    64, DynamicSize, DynamicSize>;
+template<class T> using ColMajorMatrix = C2DContainer<unsigned long, T, StorageType::ColumnMajor, 64, DynamicSize, DynamicSize>;
 
 using su2activevector = su2vector<su2double>;
 using su2activematrix = su2matrix<su2double>;
diff --git a/Common/include/toolboxes/CSquareMatrixCM.hpp b/Common/include/toolboxes/CSquareMatrixCM.hpp
@@ -0,0 +1,135 @@
+﻿/*!
+ * \file CSquareMatrixCM.hpp
+ * \brief Dense general square matrix, used for example in DG standard elements
+ *        in Column Major order storage.
+ * \author Edwin van der Weide, Pedro Gomes.
+ * \version 7.0.8 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+
+#include <vector>
+#include "../containers/C2DContainer.hpp"
+
+/*!
+ * \brief Class to store a dense general square matrix that uses the Column
+ *        Major order storage format. The code should be compiled with
+ *        LAPACK to use optimized matrix inversion and multiplication routines.
+ */
+class CSquareMatrixCM {
+  static_assert(ColMajorMatrix<passivedouble>::Storage == StorageType::ColumnMajor,
+                "Column major storage is assumed for LAPACK.");
+private:
+  ColMajorMatrix<passivedouble> mat;  /*!< \brief Storage of the actual matrix. */
+
+public:
+
+  /*!
+   * \brief Default constructor. Nothing to be done.
+   */
+  CSquareMatrixCM() = default;
+
+  /*!
+   * \overload
+   * \brief Overloaded constructor, which allocates the memory to store
+   *        the matrix.
+   * \param[in] N - Number of rows and colums of the matrix.
+   */
+  CSquareMatrixCM(int N) {Initialize(N);}
+
+  /*!
+   * \brief Operator, which makes available the given matrix element as a reference.
+   * \param[in] i   - Row index of the matrix element.
+   * \param[in] j   - Column index of the matrix element.
+   * \return          Reference to element (i,j).
+   */
+  inline passivedouble& operator() (int i, int j) {return mat(i,j);}
+
+  /*!
+   * \brief Operator, which makes available the given matrix element as a const reference.
+   * \param[in] i   - Row index of the matrix element.
+   * \param[in] j   - Column index of the matrix element.
+   * \return          Constant reference to element (i,j).
+   */
+  inline const passivedouble& operator() (int i, int j) const {return mat(i,j);}
+
+  /*!
+   * \brief Function, which makes available a reference to the actual matrix.
+   * \return A reference to mat.
+   */
+  inline ColMajorMatrix<passivedouble>& GetMat() {return mat;}
+
+  /*!
+   * \brief Function, which makes available a const reference to the actual matrix.
+   * \return A const reference to mat.
+   */
+  inline const ColMajorMatrix<passivedouble>& GetMat() const {return mat;}
+
+  /*!
+   * \brief Function, which allocates the memory for the matrix.
+   * \param[in] N - Number of rows and colums of the matrix.
+   */
+  inline void Initialize(int N) {mat.resize(N,N);}
+
+  /*!
+   * \brief Function, which makes available the size of the matrix.
+   * \return The number of rows, columns of the matrix.
+   */
+  inline int Size() const {return mat.rows();}
+
+  /*!
+   * \brief Function, which carries out the matrix produc of the current matrix
+   *        with mat_in and stores the result in mat_out.
+   * \param[in]  side    - left: mat_out = this * mat_in, right: mat_out = mat_in * this
+   * \param[in]  mat_in  - Matrix to be multiplied by the current matrix.
+   * \param[out] mat_out - Matrix to store the result of the multiplication.
+   */
+  void MatMatMult(const char                          side,
+                  const ColMajorMatrix<passivedouble> &mat_in,
+                  ColMajorMatrix<passivedouble>       &mat_out) const;
+
+  /*!
+   * \brief Naive matrix-vector multiplication with general type.
+   */
+  template<class ForwardIt>
+  void MatVecMult(ForwardIt vec_in, ForwardIt vec_out) const
+  {
+    for (int i = 0; i < Size(); ++i) {
+      *vec_out = 0.0;
+      auto vec = vec_in;
+      for (int k = 0; k < Size(); ++k)
+        *vec_out += *(vec++) * mat(i,k);
+      ++vec_out;
+    }
+  }
+
+  /*!
+   * \brief Function, which inverts the matrix in-place.
+   */
+  void Invert();
+
+  /*!
+   * \brief Function, which transposes the matrix in-place.
+   */
+  void Transpose();
+
+};
diff --git a/Common/include/toolboxes/CSymmetricMatrix.hpp b/Common/include/toolboxes/CSymmetricMatrix.hpp
@@ -42,7 +42,6 @@ class CSymmetricMatrix {
   // Not optimized dense matrix factorization and inversion for portability.
   void CalcInv(bool is_spd);
   void CholeskyDecompose();
-  void LUDecompose(su2passivematrix& decomp, std::vector<int>& perm) const;
   // Matrix inversion using LAPACK routines (LDLT and LLT factorization).
   void CalcInv_sytri();
   void CalcInv_potri();
diff --git a/Common/lib/Makefile.am b/Common/lib/Makefile.am
@@ -113,6 +113,7 @@ lib_sources = \
   ../src/toolboxes/CLinearPartitioner.cpp \
   ../src/toolboxes/C1DInterpolation.cpp \
   ../src/toolboxes/CSymmetricMatrix.cpp \
+  ../src/toolboxes/CSquareMatrixCM.cpp \
   ../src/toolboxes/MMS/CVerificationSolution.cpp \
   ../src/toolboxes/MMS/CIncTGVSolution.cpp \
   ../src/toolboxes/MMS/CInviscidVortexSolution.cpp \
diff --git a/Common/src/toolboxes/CSquareMatrixCM.cpp b/Common/src/toolboxes/CSquareMatrixCM.cpp
@@ -0,0 +1,144 @@
+/*!
+ * \file CSquareMatrixCM.cpp
+ * \brief Implementation of dense matrix helper class in Column Major order (see hpp).
+ * \author Edwin van der Weide, Pedro Gomes.
+ * \version 7.0.8 "Blackbird"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "../../include/toolboxes/CSquareMatrixCM.hpp"
+#include "../../include/mpi_structure.hpp"
+#include "../../include/blas_structure.hpp"
+
+using namespace std;
+
+#if defined(HAVE_MKL)
+#include "mkl.h"
+#ifndef HAVE_LAPACK
+#define HAVE_LAPACK
+#endif
+#elif defined(HAVE_LAPACK)
+/*--- Lapack / Blas routines used in CSquareMatrixCM. ---*/
+extern "C" void dgetrf_(const int*, const int*, passivedouble*, const int*,
+                        int*, int*);
+extern "C" void dgetri_(const int*, passivedouble*, const int*, int*,
+                        passivedouble*, const int*, int*);
+extern "C" void dgemm_(char*, char*, const int*, const int*, const int*,
+                       const passivedouble*, const passivedouble*,
+                       const int *, const passivedouble*, const int*,
+                       const passivedouble*, passivedouble*, const int*);
+#define DGEMM dgemm_
+#endif
+
+void CSquareMatrixCM::Transpose() {
+
+  for(int j=1; j<Size(); ++j)
+    for(int i=0; i<j; ++i)
+      swap(mat(i,j), mat(j,i));
+}
+
+void CSquareMatrixCM::Invert() {
+
+#ifdef HAVE_LAPACK
+
+  /*--- Computation of the inverse using the Lapack routines. ---*/
+  int sz = Size();
+  int info;
+  vector<int> ipiv(sz);
+  vector<passivedouble> work(sz);
+
+  dgetrf_(&sz, &sz, mat.data(), &sz, ipiv.data(), &info);
+  if(info != 0) SU2_MPI::Error(string("Matrix is singular"), CURRENT_FUNCTION);
+
+  dgetri_(&sz, mat.data(), &sz, ipiv.data(), work.data(), &sz, &info);
+  if(info != 0) SU2_MPI::Error(string("Matrix inversion failed"), CURRENT_FUNCTION);
+
+#else
+  CBlasStructure::inverse(Size(), mat);
+#endif
+}
+
+void CSquareMatrixCM::MatMatMult(const char                          side,
+                                        const ColMajorMatrix<passivedouble> &mat_in,
+                                        ColMajorMatrix<passivedouble>       &mat_out) const {
+
+  /*--- Check the type of multiplication to be carried out. ---*/
+  if (side == 'L' || side == 'l') {
+
+    /*--- Left side: mat_out = this * mat_in. Set some sizes
+          and allocate the memory for mat_out. ---*/
+    const int M = Size(), N = mat_in.cols();
+    assert(M == mat_in.rows());
+
+    mat_out.resize(M,N);
+
+#ifdef HAVE_LAPACK
+
+    /*--- The Lapack/blas function dgemm is used to carry out
+          the matrix matrix multiplication. ---*/
+    passivedouble alpha = 1.0, beta = 0.0;
+    char trans = 'N';
+
+    DGEMM(&trans, &trans, &M, &N, &M, &alpha, mat.data(), &M,
+          mat_in.data(), &M, &beta, mat_out.data(), &M);
+#else
+    /*--- Naive product. ---*/
+    for (int i = 0; i < M; ++i) {
+      for (int j = 0; j < N; ++j) {
+        mat_out(i,j) = 0.0;
+        for (int k = 0; k < M; ++k)
+          mat_out(i,j) += mat(i,k) * mat_in(k,j);
+      }
+    }
+#endif
+
+  }
+  else {
+
+    /*--- Right_side: mat_out = mat_in * this. Set some sizes
+          and allocate the memory for mat_out. ---*/
+    const int M = mat_in.rows(), N = Size();
+    assert(N == mat_in.cols());
+
+    mat_out.resize(M,N);
+
+#ifdef HAVE_LAPACK
+
+    /*--- The Lapack/blas function dgemm is used to carry out
+          the matrix matrix multiplication. ---*/
+    passivedouble alpha = 1.0, beta = 0.0;
+    char trans = 'N';
+
+    DGEMM(&trans, &trans, &M, &N, &N, &alpha, mat_in.data(), &M,
+          mat.data(), &N, &beta, mat_out.data(), &M);
+#else
+    /*--- Naive product. ---*/
+    for (int i = 0; i < M; ++i) {
+      for (int j = 0; j < N; ++j) {
+        mat_out(i,j) = 0.0;
+        for (int k = 0; k < N; ++k)
+          mat_out(i,j) += mat_in(i,k) * mat(k,j);
+      }
+    }
+#endif
+  }
+}
diff --git a/Common/src/toolboxes/CSymmetricMatrix.cpp b/Common/src/toolboxes/CSymmetricMatrix.cpp
diff --git a/Common/src/toolboxes/meson.build b/Common/src/toolboxes/meson.build
