LA: add lapack-based solveLU functions

src/LA/system_solvers_small.cpp

@@ -43,10 +43,20 @@
 # include "multiplication.hpp"
 # include "system_solvers_small.hpp"
 
+
+# include "bitpit_private_lapacke.hpp"
+
 namespace bitpit{
 
 namespace linearalgebra{
 
+namespace constants{
+
+const int ROW_MAJOR = LAPACK_ROW_MAJOR;
+const int COL_MAJOR = LAPACK_COL_MAJOR;
+
+}
+
 /*!
  * @ingroup system_solver_small
  * @{
@@ -353,6 +363,62 @@ backwardSubstitution(U, z, x);
 
 return; };
 
+// -------------------------------------------------------------------------- //
+/*!
+    Solve a linear system using Lapacke_dgesv (partial pivoting
+    and LU factorization). Containers must be correctly sized.
+    This routine works on the entire matrix, the number of considered equations
+    is equal to the leading dimension of the matrix, the same is for the r.h.s.
+    Multiple r.h.s. are not allowed.
+
+    \param[in] layout matrix layout (possible values are ROW_MAJOR and COL_MAJOR)
+    \param[in] A coeffs of square matrix (linear)
+    \param[in] B r.h.s. of the linear system
+    \param[in,out] x on output stores the solution of the linear system
+    \return information on the execution of LAPACKE_dgesv, see http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga5ee879032a8365897c3ba91e3dc8d512.html
+*/
+int solveLU(
+    int                              layout,
+    std::vector<double>                  &A,
+    std::vector<double>                  &B,
+    std::vector<double>                  &x
+) {
+
+    int info = solveLU(layout, B.size(), A.data(), B.data(), x.data()); 
+    return info;
+
+};
+
+// -------------------------------------------------------------------------- //
+/*!
+    Solve a linear system using Lapacke_degesv (partial pivoting
+    and LU factorization).
+    This routine works on the entire matrix, the number of considered equations
+    is equal to the leading dimension of the matrix, the same is for the r.h.s.
+    Multiple r.h.s. are not allowed.
+
+    \param[in] layout matrix layout (possible values are ROW_MAJOR and COL_MAJOR)
+    \param[in] matrixOrder number of equations in linear system
+    \param[in] A coeffs of square matrix (linear)
+    \param[in] B r.h.s. of the linear system
+    \param[in,out] x on output stores the solution of the linear system
+    \return information on the execution of LAPACKE_dgesv, see http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga5ee879032a8365897c3ba91e3dc8d512.html
+*/
+int solveLU(
+    int                              layout,
+    int                         matrixOrder,
+    double                               *A,
+    double                               *B,
+    double                               *x
+) {
+
+    std::vector<int> ipiv(matrixOrder);
+    int ldb = (layout == constants::ROW_MAJOR ? 1 : matrixOrder);
+    int info = LAPACKE_dgesv(layout, matrixOrder, 1, A, matrixOrder, ipiv.data(), B, ldb);
+    return info;
+
+};
+
 /*!
  * @}
  */

src/LA/system_solvers_small.hpp

@@ -57,6 +57,13 @@ namespace bitpit{
  */
 namespace linearalgebra{
 
+namespace constants{
+
+extern const int ROW_MAJOR;
+extern const int COL_MAJOR;
+
+}
+
 template <class T>
 void cramer(                                                                  // Solve linear system using Cramer's rule
     std::vector< std::vector < T > >            &,                            // (input) coeff. matrix
@@ -125,6 +132,21 @@ void solveLU(                                                                 //
     std::array<double, m>                       &                             // (input/output) Solution
 );
 
+int solveLU(                                                                 // Solve linear system using Lapack DGESV
+    int                                          ,                            // (input) Matrix layout
+    std::vector<double>                         &,                            // (input) Input coeffs. matrix (linear)
+    std::vector<double>                         &,                            // (input) Source term
+    std::vector<double>                         &                             // (input/output) Solution
+);
+
+int solveLU(                                                                 // Solve linear system using Lapack DGESV
+    int                                          ,                            // (input) Matrix layout
+    int                                          ,                            // (input) Matrix order
+    double                                      *,                            // (input) Pointer to input coeffs. matrix (linear)
+    double                                      *,                            // (input) Pointer to source term
+    double                                      *                             // (input/output) Pointer to solution
+);
+
 }
 
 }