diff --git a/src/LA/matrix_utilities.tpp b/src/LA/matrix_utilities.tpp index 635913e4eb..43dd7561cc 100644 --- a/src/LA/matrix_utilities.tpp +++ b/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 } diff --git a/src/LA/system_solvers_small.cpp b/src/LA/system_solvers_small.cpp index ddb007457a..9cafcb73c7 100644 --- a/src/LA/system_solvers_small.cpp +++ b/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 &A, + std::vector &B +) { + + std::vector 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; + +}; + /*! * @} */ diff --git a/src/LA/system_solvers_small.hpp b/src/LA/system_solvers_small.hpp index af01b78973..d959c89ce2 100644 --- a/src/LA/system_solvers_small.hpp +++ b/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 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 & // (input/output) Solution ); +int solveLU( // Solve linear system using Lapack DGESV + int , // (input) Matrix layout + std::vector &, // (input) Input coeffs. matrix (linear) + std::vector & // (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 +); + } }