text: /* ----------------------------------------------------------------------- *//**
   *
   * @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
  