Add lapacke-based solveLU

src/LA/matrix_utilities.tpp

@@ -261,7 +261,6 @@ m = A.size();
 if (m == 0) { return; }
 n = A[0].size();
 if (n == 0) { return; }
-i--; j--;
 if ((i >= m) || (i < 0)) {
     return;
 }
@@ -275,9 +274,9 @@ if ((j >= n) || (j < 0)) {
 
 // Resize output variables
 B.resize(m-1);
-for (i = 0; i < m-1; ++i) {
-    B[i].resize(n-1, 0.0);
-} //next i
+for (int p = 0; p < m-1; ++p) {
+    B[p].resize(n-1, 0.0);
+} //next p
 
 // Extract complement
 for (l = 0; l < i; l++) {
@@ -332,7 +331,6 @@ int            l, k;
 // ========================================================================== //
 if (m == 0) { return; }
 if (n == 0) { return; }
-i--; j--;
 if ((i >= (long) m) || (i < 0)) {
     return;
 }
@@ -692,9 +690,6 @@ int                     m, n;
 T                       d = (T) 1.0e+18;
 std::vector< std::vector < T > >  C;
 
-// Counters
-int                     i;
-
 // ========================================================================== //
 // CHECK INPUT                                                                //
 // ========================================================================== //
@@ -722,9 +717,9 @@ else if (m == 3) {
 }
 else {
     d = (T) 0.0;
-    for (i = 0; i < m; i++) {
-        complement(1, i+1, A, C);
-        d += pow(-1.0, i+2) * A[0][i] * det(C);
+    for (int i = 0; i < m; i++) {
+        complement(0, i, A, C);
+        d += pow(-1.0, i) * A[0][i] * det(C);
     } //next i
 }
 
@@ -781,8 +776,8 @@ else if (m == 3) {
 else {
     std::array< std::array < double, m-1 >, m-1 >     C;
     for (i = 0; i < m; i++) {
-        complement(1, i+1, A, C);
-        d += pow(-1.0, i+2) * A[0][i] * det(C);
+        complement(0, i, A, C);
+        d += pow(-1.0, i) * A[0][i] * det(C);
     } //next i
 }
 

src/LA/system_solvers_small.cpp

@@ -43,10 +43,19 @@
 # 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 +362,72 @@ backwardSubstitution(U, z, x);
 
 return; };
 
+// -------------------------------------------------------------------------- //
+/*!
+    Solve a linear system using partial pivoting and LU factorization (the
+    LAPACK function Lapacke_dgesv is used).
+
+    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. Containers MUST be correctly sized.
+    Permutation matrix not provided.
+
+    This function is intended for solution computation ONLY.
+
+    \param[in] layout matrix layout (possible values are ROW_MAJOR and COL_MAJOR)
+    \param[in,out] A in input the coefficients of square matrix (linear), the
+    solver will the use the matrix as a storage for temporary data needed during
+    the solution of the system, hence on output the matrix coefficients will be
+    overwritten and the original matrix coefficients cannot be recovered
+    \param[in,out] B in input r.h.s. of the linear system, in output the solution.
+    \return information on the execution of LAPACKE_dgesv, see LAPACK  documentation
+*/
+int solveLU(
+    int                              layout,
+    std::vector<double>                  &A,
+    std::vector<double>                  &B
+) {
+
+    std::vector<int> ipiv(B.size());
+    int info = solveLU(layout, B.size(), A.data(), B.data(),ipiv.data());
+    return info;
+
+};
+
+// -------------------------------------------------------------------------- //
+/*!
+    Solve a linear system using partial pivoting and LU factorization (the
+    LAPACK function Lapacke_dgesv is used).
+
+    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. Containers MUST be correctly sized.
+
+    This function is intended for both solution and factorization computation.
+
+    \param[in] layout matrix layout (possible values are ROW_MAJOR and COL_MAJOR).
+    \param[in] matrixOrder number of equations in linear system.
+    \param[in,out] A in input the coefficients of square matrix (linear), in output the
+    factors L and U, A = P * L * U.
+    \param[in,out] B in input r.h.s. of the linear system, in output the solution.
+    \param[out] ipiv the pivot indices that define the permutation matrix P, see LAPACK
+    documentation.
+    \return information on the execution of LAPACKE_dgesv, see LAPACK documentation
+*/
+int solveLU(
+    int                              layout,
+    int                         matrixOrder,
+    double                               *A,
+    double                               *B,
+    int                               *ipiv
+) {
+
+    int ldb = (layout == constants::ROW_MAJOR ? 1 : matrixOrder);
+    int info = LAPACKE_dgesv(layout, matrixOrder, 1, A, matrixOrder, ipiv, 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,20 @@ 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/output) Source term / 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/output) Pointer to source term / solution
+    int                                         *                             // (output) Pivot indeces of the permutation matrix
+);
+
 }
 
 }