| /* ----------------------------------------------------------------------- *//** |
| * |
| * @file sparse_linear_systems.cpp |
| * |
| * @brief sparse linear systems |
| * |
| * We implement sparse linear systems using the direct. |
| * |
| *//* ----------------------------------------------------------------------- */ |
| #include <limits> |
| #include <dbconnector/dbconnector.hpp> |
| #include <modules/shared/HandleTraits.hpp> |
| #include <modules/prob/boost.hpp> |
| #include "sparse_linear_systems.hpp" |
| #include <Eigen/Sparse> |
| |
| // Common SQL functions used by all modules |
| //#include "sparse_linear_systems_states.hpp" |
| |
| namespace madlib { |
| |
| // Use Eigen |
| using namespace dbal::eigen_integration; |
| using namespace Eigen; |
| |
| namespace modules { |
| |
| // Import names from other MADlib modules |
| using dbal::NoSolutionFoundException; |
| |
| namespace linear_systems{ |
| |
| |
| // --------------------------------------------------------------------------- |
| // Direct sparse Linear System States |
| // --------------------------------------------------------------------------- |
| |
| // Function protofype of Internal functions |
| AnyType direct_sparse_stateToResult( |
| const Allocator &inAllocator, |
| const ColumnVector &x, |
| const double residual); |
| |
| |
| /** |
| * @brief States for linear systems |
| * |
| * TransitionState encapsualtes the transition state during the |
| * logistic-regression aggregate function. To the database, the state is |
| * exposed as a single DOUBLE PRECISION array, to the C++ code it is a proper |
| * object containing scalars and vectors. |
| * |
| * Note: We assume that the DOUBLE PRECISION array is initialized by the |
| * database with length at least 5, and all elemenets are 0. |
| * |
| */ |
| template <class Handle> |
| class SparseDirectLinearSystemTransitionState { |
| template <class OtherHandle> |
| friend class SparseDirectLinearSystemTransitionState; |
| |
| public: |
| SparseDirectLinearSystemTransitionState(const AnyType &inArray) |
| : mStorage(inArray.getAs<Handle>()) { |
| |
| rebind( static_cast<uint32_t>(mStorage[1]), |
| static_cast<uint32_t>(mStorage[2])); |
| } |
| |
| /** |
| * @brief Convert to backend representation |
| * |
| * We define this function so that we can use State in the |
| * argument list and as a return type. |
| */ |
| inline operator AnyType() const { |
| return mStorage; |
| } |
| |
| /** |
| * @brief Initialize the conjugate-gradient state. |
| * |
| * This function is only called for the first iteration, for the first row. |
| */ |
| inline void initialize(const Allocator &inAllocator, |
| uint32_t innumVars, |
| uint32_t innumEquations, |
| uint32_t inNNZA |
| ) { |
| // Array size does not depend in numVars |
| mStorage = inAllocator.allocateArray<double, dbal::AggregateContext, |
| dbal::DoZero, dbal::ThrowBadAlloc>(arraySize(innumEquations, inNNZA)); |
| rebind(innumEquations, inNNZA); |
| numVars = innumVars; |
| numEquations = innumEquations; |
| NNZA = inNNZA; |
| } |
| |
| /** |
| * @brief We need to support assigning the previous state |
| */ |
| template <class OtherHandle> |
| SparseDirectLinearSystemTransitionState &operator=( |
| const SparseDirectLinearSystemTransitionState<OtherHandle> &inOtherState) { |
| |
| for (size_t i = 0; i < mStorage.size(); i++) |
| mStorage[i] = inOtherState.mStorage[i]; |
| return *this; |
| } |
| |
| /** |
| * @brief Merge with another State object by copying the intra-iteration |
| * fields |
| */ |
| template <class OtherHandle> |
| SparseDirectLinearSystemTransitionState &operator+=( |
| const SparseDirectLinearSystemTransitionState<OtherHandle> &inOtherState) { |
| |
| if (mStorage.size() != inOtherState.mStorage.size() || |
| numVars != inOtherState.numVars || |
| NNZA != inOtherState.NNZA|| |
| numEquations != inOtherState.numEquations) |
| throw std::logic_error("Internal error: Incompatible transition " |
| "states"); |
| |
| b += inOtherState.b; |
| b_stored += inOtherState.b_stored; |
| |
| // Merge state for sparse loading is not an array add operation |
| // but it is an array-append operation |
| for (uint32_t i=0; i < inOtherState.nnz_processed; i++){ |
| r(nnz_processed + i) += inOtherState.r(i); |
| c(nnz_processed + i) += inOtherState.c(i); |
| v(nnz_processed + i) += inOtherState.v(i); |
| } |
| nnz_processed += inOtherState.nnz_processed; |
| return *this; |
| } |
| |
| /** |
| * @brief Reset the inter-iteration fields. |
| */ |
| inline void reset() { |
| nnz_processed = 0; |
| r.fill(0); |
| c.fill(0); |
| v.fill(0); |
| b.fill(0); |
| b_stored.fill(0); |
| } |
| |
| private: |
| static inline size_t arraySize(const uint32_t innumEquations, |
| const uint32_t inNNZA) { |
| return 5 + 3 * inNNZA + 2 * innumEquations; |
| } |
| |
| /** |
| * @brief Rebind to a new storage array |
| * |
| * @param innumVars The number of independent variables. |
| * |
| * Array layout (iteration refers to one aggregate-function call): |
| * Inter-iteration components (updated in final function): |
| * - 0: numVars (Total number of variables) |
| * - 1: numEquations (Total number of equations) |
| * - 2: nnZ (Total number of non-zeros) |
| * - 3: algorithm |
| * - 4: nnz_processed Number of non-zeros processed by a node |
| * - 4: b (RHS vector) |
| * - 4 + 1 * numEquations: r (LHS matrix row) |
| * - 4 + 1 * numEquations + 1 * nnzA: c (LHS matrix column) |
| * - 4 + 1 * numEquations + 2 * nnzA: v (LHS matrix value) |
| */ |
| void rebind(uint32_t innumEquations, uint32_t inNNZA) { |
| numVars.rebind(&mStorage[0]); |
| numEquations.rebind(&mStorage[1]); |
| NNZA.rebind(&mStorage[2]); |
| algorithm.rebind(&mStorage[3]); |
| nnz_processed.rebind(&mStorage[4]); |
| b.rebind(&mStorage[5], innumEquations); |
| b_stored.rebind(&mStorage[5 + innumEquations], innumEquations); |
| r.rebind(&mStorage[5 + 2*innumEquations], inNNZA); |
| c.rebind(&mStorage[5 + 2*innumEquations + inNNZA], inNNZA); |
| v.rebind(&mStorage[5 + 2*innumEquations + 2 * inNNZA], inNNZA); |
| } |
| |
| Handle mStorage; |
| |
| public: |
| typename HandleTraits<Handle>::ReferenceToUInt32 numVars; |
| typename HandleTraits<Handle>::ReferenceToUInt32 numEquations; |
| typename HandleTraits<Handle>::ReferenceToUInt32 NNZA; |
| typename HandleTraits<Handle>::ReferenceToUInt32 nnz_processed; |
| typename HandleTraits<Handle>::ReferenceToUInt32 algorithm; |
| |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap b_stored; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap b; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap r; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap c; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap v; |
| }; |
| |
| |
| /** |
| * @brief Perform the sparse linear system transition step |
| */ |
| AnyType |
| sparse_direct_linear_system_transition::run(AnyType &args) { |
| |
| SparseDirectLinearSystemTransitionState<MutableArrayHandle<double> > |
| state = args[0]; |
| int32_t row_id = args[1].getAs<int32_t>(); |
| int32_t col_id = args[2].getAs<int32_t>(); |
| double value = args[3].getAs<double>(); |
| double _b = args[4].getAs<double>(); |
| |
| |
| |
| // When the segment recieves the first non-zero in the sparse matrix |
| // we initialize the state |
| if (state.nnz_processed == 0) { |
| |
| int32_t num_equations = args[5].getAs<int32_t>(); |
| int32_t num_vars = args[6].getAs<int32_t>(); |
| int32_t total_nnz = args[7].getAs<int32_t>(); |
| int algorithm = args[8].getAs<int>(); |
| state.initialize(*this, |
| num_vars, |
| num_equations, |
| total_nnz |
| ); |
| state.algorithm = algorithm; |
| state.NNZA = total_nnz; |
| state.numVars = num_vars; |
| state.numEquations = num_equations; |
| state.b_stored.setZero(); |
| state.b.setZero(); |
| |
| } |
| |
| // Now do the transition step |
| // First we create a block of zeros in memory |
| // and then add the vector in the appropriate position |
| ColumnVector bvec(static_cast<uint32_t>(state.numEquations)); |
| ColumnVector rvec(static_cast<uint32_t>(state.NNZA)); |
| ColumnVector cvec(static_cast<uint32_t>(state.NNZA)); |
| ColumnVector vvec(static_cast<uint32_t>(state.NNZA)); |
| |
| bvec.setZero(); |
| rvec.setZero(); |
| cvec.setZero(); |
| vvec.setZero(); |
| |
| rvec(state.nnz_processed) = row_id; |
| cvec(state.nnz_processed) = col_id; |
| vvec(state.nnz_processed) = value; |
| |
| if (state.b_stored(row_id) == 0) { |
| bvec(row_id) = _b; |
| state.b_stored(row_id) = 1; |
| state.b += bvec; |
| } |
| |
| // Build the vector & matrices based on row_id |
| state.r += rvec; |
| state.c += cvec; |
| state.v += vvec; |
| state.nnz_processed++; |
| |
| return state; |
| } |
| |
| |
| /** |
| * @brief Perform the perliminary aggregation function: Merge transition states |
| */ |
| AnyType |
| sparse_direct_linear_system_merge_states::run(AnyType &args) { |
| SparseDirectLinearSystemTransitionState<MutableArrayHandle<double> > stateLeft = args[0]; |
| SparseDirectLinearSystemTransitionState<ArrayHandle<double> > stateRight = args[1]; |
| |
| // We first handle the trivial case where this function is called with one |
| // of the states being the initial state |
| if (stateLeft.numEquations == 0) |
| return stateRight; |
| else if (stateRight.numEquations == 0) |
| return stateLeft; |
| |
| // Merge states together and return |
| stateLeft += stateRight; |
| return stateLeft; |
| } |
| |
| |
| |
| /** |
| * @brief Perform the linear system final step |
| */ |
| AnyType |
| sparse_direct_linear_system_final::run(AnyType &args) { |
| // We request a mutable object. Depending on the backend, this might perform |
| // a deep copy. |
| SparseDirectLinearSystemTransitionState<MutableArrayHandle<double> > state = args[0]; |
| |
| // Aggregates that haven't seen any data just return Null. |
| if (state.numEquations == 0) |
| return Null(); |
| |
| // Eigen works better when you reserve the number of nnz |
| // Note: When calling reserve(), it is not required that nnz is the |
| // exact number of nonzero elements in the final matrix. However, an |
| // exact estimation will avoid multiple reallocations during the |
| // insertion phase. |
| SparseMatrix A(state.numEquations, state.numVars); |
| A.reserve(state.NNZA); |
| ColumnVector x; |
| |
| for (uint32_t i=0; i < state.NNZA; i++){ |
| A.insert(static_cast<uint32_t>(state.r(i)), |
| static_cast<uint32_t>(state.c(i))) = state.v(i); |
| } |
| |
| // Switch case needs scoping in C++ if you want to declare inside it |
| // Unfortunately, this means that I have to write the code to call the |
| // solver in each switch case separtely |
| switch (state.algorithm) { |
| case 1:{ |
| Eigen::SimplicialLLT<SparseMatrix> solver; |
| solver.compute(A); |
| x = solver.solve(state.b); |
| break; |
| } |
| case 2:{ |
| Eigen::SimplicialLDLT<SparseMatrix> solver; |
| solver.compute(A); |
| x = solver.solve(state.b); |
| break; |
| } |
| } |
| |
| |
| // It is unclear whether I should do the residual computation in-memory or |
| // in-database (as done in the dense linear systems case) |
| ColumnVector bvec = state.b; |
| double residual = (A*x-bvec).norm() / bvec.norm(); |
| |
| // Compute the residual |
| return direct_sparse_stateToResult(*this, x, residual); |
| } |
| |
| /** |
| * @brief Helper function that computes the final statistics |
| */ |
| |
| AnyType direct_sparse_stateToResult( |
| const Allocator &inAllocator, |
| const ColumnVector &x, |
| const double residual) { |
| |
| MutableNativeColumnVector solution( |
| inAllocator.allocateArray<double>(x.size())); |
| |
| for (Index i = 0; i < x.size(); ++i) { |
| solution(i) = x(i); |
| } |
| |
| // Return all coefficients, standard errors, etc. in a tuple |
| AnyType tuple; |
| tuple << solution |
| << residual |
| << Null(); |
| |
| return tuple; |
| } |
| |
| |
| // --------------------------------------------------------------------------- |
| // In-memory iterative sparse Linear System States |
| // --------------------------------------------------------------------------- |
| |
| // Function protofype of Internal functions |
| AnyType inmem_iterative_sparse_stateToResult( |
| const Allocator &inAllocator, |
| const ColumnVector &x, |
| const int iters, |
| const double residual_norm); |
| |
| |
| /** |
| * @brief States for linear systems |
| * |
| * TransitionState encapsualtes the transition state during the |
| * logistic-regression aggregate function. To the database, the state is |
| * exposed as a single DOUBLE PRECISION array, to the C++ code it is a proper |
| * object containing scalars and vectors. |
| * |
| * Note: We assume that the DOUBLE PRECISION array is initialized by the |
| * database with length at least 5, and all elemenets are 0. |
| * |
| */ |
| template <class Handle> |
| class SparseInMemIterativeLinearSystemTransitionState { |
| template <class OtherHandle> |
| friend class SparseInMemIterativeLinearSystemTransitionState; |
| |
| public: |
| SparseInMemIterativeLinearSystemTransitionState(const AnyType &inArray) |
| : mStorage(inArray.getAs<Handle>()) { |
| |
| rebind( static_cast<uint32_t>(mStorage[1]), |
| static_cast<uint32_t>(mStorage[2])); |
| } |
| |
| /** |
| * @brief Convert to backend representation |
| * |
| * We define this function so that we can use State in the |
| * argument list and as a return type. |
| */ |
| inline operator AnyType() const { |
| return mStorage; |
| } |
| |
| /** |
| * @brief Initialize the conjugate-gradient state. |
| * |
| * This function is only called for the first iteration, for the first row. |
| */ |
| inline void initialize(const Allocator &inAllocator, |
| uint32_t innumVars, |
| uint32_t innumEquations, |
| uint32_t inNNZA |
| ) { |
| // Array size does not depend in numVars |
| mStorage = inAllocator.allocateArray<double, dbal::AggregateContext, |
| dbal::DoZero, dbal::ThrowBadAlloc>(arraySize(innumEquations, inNNZA)); |
| rebind(innumEquations, inNNZA); |
| numVars = innumVars; |
| numEquations = innumEquations; |
| NNZA = inNNZA; |
| } |
| |
| /** |
| * @brief We need to support assigning the previous state |
| */ |
| template <class OtherHandle> |
| SparseInMemIterativeLinearSystemTransitionState &operator=( |
| const SparseInMemIterativeLinearSystemTransitionState<OtherHandle> &inOtherState) { |
| |
| for (size_t i = 0; i < mStorage.size(); i++) |
| mStorage[i] = inOtherState.mStorage[i]; |
| return *this; |
| } |
| |
| /** |
| * @brief Merge with another State object by copying the intra-iteration |
| * fields |
| */ |
| template <class OtherHandle> |
| SparseInMemIterativeLinearSystemTransitionState &operator+=( |
| const SparseInMemIterativeLinearSystemTransitionState<OtherHandle> &inOtherState) { |
| |
| if (mStorage.size() != inOtherState.mStorage.size() || |
| numVars != inOtherState.numVars || |
| NNZA != inOtherState.NNZA|| |
| numEquations != inOtherState.numEquations) |
| throw std::logic_error("Internal error: Incompatible transition " |
| "states"); |
| |
| b += inOtherState.b; |
| b_stored += inOtherState.b_stored; |
| |
| // Merge state for sparse loading is not an array add operation |
| // but it is an array-append operation |
| for (uint32_t i=0; i < inOtherState.nnz_processed; i++){ |
| r(nnz_processed + i) += inOtherState.r(i); |
| c(nnz_processed + i) += inOtherState.c(i); |
| v(nnz_processed + i) += inOtherState.v(i); |
| } |
| nnz_processed += inOtherState.nnz_processed; |
| return *this; |
| } |
| |
| /** |
| * @brief Reset the inter-iteration fields. |
| */ |
| inline void reset() { |
| nnz_processed = 0; |
| r.fill(0); |
| c.fill(0); |
| v.fill(0); |
| b.fill(0); |
| b_stored.fill(0); |
| } |
| |
| private: |
| static inline size_t arraySize(const uint32_t innumEquations, |
| const uint32_t inNNZA) { |
| return 7 + 3 * inNNZA + 2 * innumEquations; |
| } |
| |
| /** |
| * @brief Rebind to a new storage array |
| * |
| * @param innumVars The number of independent variables. |
| * |
| * Array layout (iteration refers to one aggregate-function call): |
| * Inter-iteration components (updated in final function): |
| * - 0: numVars (Total number of variables) |
| * - 1: numEquations (Total number of equations) |
| * - 2: nnZ (Total number of non-zeros) |
| * - 3: algorithm |
| * - 4: maxIter |
| * - 5: termToler |
| * - 6: nnz_processed Number of non-zeros processed by a node |
| * - 7: b (RHS vector) |
| * - 7: b_stored (Boolean: Has the b_vector already been loaded?) |
| * - 7 + 2 * numEquations: r (LHS matrix row) |
| * - 7 + 2 * numEquations + 1 * nnzA: c (LHS matrix column) |
| * - 7 + 2 * numEquations + 2 * nnzA: v (LHS matrix value) |
| */ |
| void rebind(uint32_t innumEquations, uint32_t inNNZA) { |
| numVars.rebind(&mStorage[0]); |
| numEquations.rebind(&mStorage[1]); |
| NNZA.rebind(&mStorage[2]); |
| algorithm.rebind(&mStorage[3]); |
| nnz_processed.rebind(&mStorage[4]); |
| maxIter.rebind(&mStorage[5]); |
| termToler.rebind(&mStorage[6]); |
| b.rebind(&mStorage[7], innumEquations); |
| b_stored.rebind(&mStorage[7 + innumEquations], innumEquations); |
| r.rebind(&mStorage[7 + 2*innumEquations], inNNZA); |
| c.rebind(&mStorage[7 + 2*innumEquations + inNNZA], inNNZA); |
| v.rebind(&mStorage[7 + 2*innumEquations + 2 * inNNZA], inNNZA); |
| } |
| |
| Handle mStorage; |
| |
| public: |
| typename HandleTraits<Handle>::ReferenceToUInt32 numVars; |
| typename HandleTraits<Handle>::ReferenceToUInt32 numEquations; |
| typename HandleTraits<Handle>::ReferenceToUInt32 NNZA; |
| typename HandleTraits<Handle>::ReferenceToUInt32 nnz_processed; |
| typename HandleTraits<Handle>::ReferenceToUInt32 algorithm; |
| typename HandleTraits<Handle>::ReferenceToUInt32 maxIter; |
| typename HandleTraits<Handle>::ReferenceToDouble termToler; |
| |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap b_stored; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap b; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap r; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap c; |
| typename HandleTraits<Handle>::ColumnVectorTransparentHandleMap v; |
| }; |
| |
| |
| /** |
| * @brief Perform the sparse linear system transition step |
| */ |
| AnyType |
| sparse_inmem_iterative_linear_system_transition::run(AnyType &args) { |
| |
| SparseInMemIterativeLinearSystemTransitionState<MutableArrayHandle<double> > |
| state = args[0]; |
| int32_t row_id = args[1].getAs<int32_t>(); |
| int32_t col_id = args[2].getAs<int32_t>(); |
| double value = args[3].getAs<double>(); |
| double _b = args[4].getAs<double>(); |
| |
| |
| // When the segment recieves the first non-zero in the sparse matrix |
| // we initialize the state |
| if (state.nnz_processed == 0) { |
| |
| int32_t num_equations = args[5].getAs<int32_t>(); |
| int32_t num_vars = args[6].getAs<int32_t>(); |
| int32_t total_nnz = args[7].getAs<int32_t>(); |
| int algorithm = args[8].getAs<int>(); |
| int32_t maxIter = args[9].getAs<int>(); |
| double termToler = args[10].getAs<double>(); |
| |
| state.initialize(*this, |
| num_vars, |
| num_equations, |
| total_nnz |
| ); |
| |
| state.algorithm = algorithm; |
| state.NNZA = total_nnz; |
| state.numVars = num_vars; |
| state.numEquations = num_equations; |
| state.maxIter = maxIter; |
| state.termToler = termToler; |
| state.b_stored.setZero(); |
| state.b.setZero(); |
| |
| } |
| |
| // Now do the transition step |
| // First we create a block of zeros in memory |
| // and then add the vector in the appropriate position |
| ColumnVector bvec(static_cast<uint32_t>(state.numEquations)); |
| ColumnVector rvec(static_cast<uint32_t>(state.NNZA)); |
| ColumnVector cvec(static_cast<uint32_t>(state.NNZA)); |
| ColumnVector vvec(static_cast<uint32_t>(state.NNZA)); |
| |
| bvec.setZero(); |
| rvec.setZero(); |
| cvec.setZero(); |
| vvec.setZero(); |
| |
| rvec(state.nnz_processed) = row_id; |
| cvec(state.nnz_processed) = col_id; |
| vvec(state.nnz_processed) = value; |
| |
| if (state.b_stored(row_id) == 0) { |
| bvec(row_id) = _b; |
| state.b_stored(row_id) = 1; |
| state.b += bvec; |
| } |
| |
| // Build the vector & matrices based on row_id |
| state.r += rvec; |
| state.c += cvec; |
| state.v += vvec; |
| state.nnz_processed++; |
| |
| return state; |
| } |
| |
| |
| /** |
| * @brief Perform the perliminary aggregation function: Merge transition states |
| */ |
| AnyType |
| sparse_inmem_iterative_linear_system_merge_states::run(AnyType &args) { |
| SparseInMemIterativeLinearSystemTransitionState<MutableArrayHandle<double> > stateLeft = args[0]; |
| SparseInMemIterativeLinearSystemTransitionState<ArrayHandle<double> > stateRight = args[1]; |
| |
| // We first handle the trivial case where this function is called with one |
| // of the states being the initial state |
| if (stateLeft.numEquations == 0) |
| return stateRight; |
| else if (stateRight.numEquations == 0) |
| return stateLeft; |
| |
| // Merge states together and return |
| stateLeft += stateRight; |
| return stateLeft; |
| } |
| |
| |
| |
| /** |
| * @brief Perform the linear system final step |
| */ |
| AnyType |
| sparse_inmem_iterative_linear_system_final::run(AnyType &args) { |
| // We request a mutable object. Depending on the backend, this might perform |
| // a deep copy. |
| SparseInMemIterativeLinearSystemTransitionState<MutableArrayHandle<double> > state = args[0]; |
| |
| // Aggregates that haven't seen any data just return Null. |
| if (state.numEquations == 0) |
| return Null(); |
| |
| // Eigen works better when you reserve the number of nnz |
| // Note: When calling reserve(), it is not required that nnz is the |
| // exact number of nonzero elements in the final matrix. However, an |
| // exact estimation will avoid multiple reallocations during the |
| // insertion phase. |
| SparseMatrix A(state.numEquations, state.numVars); |
| A.reserve(state.NNZA); |
| ColumnVector x; |
| |
| for (uint32_t i=0; i < state.NNZA; i++){ |
| A.insert(static_cast<uint32_t>(state.r(i)), |
| static_cast<uint32_t>(state.c(i))) = state.v(i); |
| } |
| |
| |
| int iters = 0; |
| double error = 0.; |
| |
| |
| // Switch case needs scoping in C++ if you want to declare inside it |
| // Unfortunately, this means that I have to write the code to call the |
| // solver in each switch case separtely |
| // Note: CG uses a diagonal preconditioner by default |
| switch (state.algorithm) { |
| case 1:{ |
| Eigen::ConjugateGradient<SparseMatrix> solver; |
| solver.setTolerance(state.termToler); |
| solver.setMaxIterations(static_cast<int>(state.maxIter)); |
| x = solver.compute(A).solve(state.b); |
| iters = solver.iterations(); |
| error = solver.error(); |
| break; |
| } |
| case 2:{ |
| Eigen::BiCGSTAB<SparseMatrix> solver; |
| solver.setTolerance(state.termToler); |
| solver.setMaxIterations(static_cast<int>(state.maxIter)); |
| x = solver.compute(A).solve(state.b); |
| iters = solver.iterations(); |
| error = solver.error(); |
| break; |
| } |
| // Preconditioned CG: Uses an incomplete LUT preconditioner |
| // For the incomplete LUT preconditioner. You might want to see Eigen docs |
| // for factors like fill-in which will make the algorithm more suitable |
| // for handling tougher linear systgems |
| case 3:{ |
| // 3 Arguments in template |
| // ConjugateGradient< _MatrixType, _UpLo, _Preconditioner >: |
| Eigen::ConjugateGradient<SparseMatrix, 1, |
| Eigen::IncompleteLUT<double> > solver; |
| solver.setTolerance(state.termToler); |
| solver.setMaxIterations(static_cast<int>(state.maxIter)); |
| x = solver.compute(A).solve(state.b); |
| iters = solver.iterations(); |
| error = solver.error(); |
| break; |
| } |
| case 4:{ |
| // 2 Arguments in template: |
| // ConjugateGradient< _MatrixType, _UpLo, _Preconditioner >: |
| // NOTE: BiCGSTAB is a bi-conjugate gradient that does not |
| // require a lower or upper triangular option |
| Eigen::BiCGSTAB<SparseMatrix, Eigen::IncompleteLUT<double> > solver; |
| solver.setTolerance(state.termToler); |
| solver.setMaxIterations(static_cast<int>(state.maxIter)); |
| x = solver.compute(A).solve(state.b); |
| iters = solver.iterations(); |
| error = solver.error(); |
| break; |
| } |
| } |
| |
| |
| // Compute the residual |
| return inmem_iterative_sparse_stateToResult(*this, x, iters, error); |
| } |
| |
| /** |
| * @brief Helper function that computes the final statistics |
| */ |
| |
| AnyType inmem_iterative_sparse_stateToResult( |
| const Allocator &inAllocator, |
| const ColumnVector &x, |
| const int iters, |
| const double residual_norm) { |
| |
| MutableNativeColumnVector solution( |
| inAllocator.allocateArray<double>(x.size())); |
| |
| for (Index i = 0; i < x.size(); ++i) { |
| solution(i) = x(i); |
| } |
| |
| // Return all coefficients, standard errors, etc. in a tuple |
| AnyType tuple; |
| tuple << solution |
| << residual_norm |
| << iters; |
| |
| return tuple; |
| } |
| |
| |
| |
| |
| |
| } // namespace linear_systems |
| |
| } // namespace modules |
| |
| } // namespace madlib |
