DianaSparseSolver.hpp

/* -----------------------------------------------------------------------------
 *  Diana process modelling, simulation and analysis software
 *  Copyright (c) 2005, Michael Krasnyk
 *  All rights reserved.
 *
 *  This file is a part of Diana process modelling, simulation and analysis software
 *
 *  Diana is free software; you can redistribute it and/or modify it
 *  under the terms of the GNU General Public License as published
 *  by the Free Software Foundation (see accompanying file LICENSE)
 * -----------------------------------------------------------------------------
 *  $Id: DianaSparseSolver.hpp 8160 2009-03-16 18:36:30Z miha $
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef DIANA_SPARSE_SOLVER_HPP
#define DIANA_SPARSE_SOLVER_HPP

#include <vector>
#include <Basic/DianaSparseArray.hpp>
#include "IDianaLinearAlgebraSolver.hpp"

namespace Diana {

  class DianaSparseSolver : virtual public IDianaLinearAlgebraSolver {
  public:
    DianaSparseSolver();
    virtual ~DianaSparseSolver();
    
    DianaSparseSolver(const DianaSparseArray& spar)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    DianaSparseSolver& operator=(const DianaSparseArray& spar)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    void decompose()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    void solve(const Common::Types::CapeArrayDouble& rhs,
               Common::Types::CapeArrayDouble& sol,
               Common::Types::CapeBoolean transpose=false)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    void solve(int n, const double* rhs, double* sol,
               Common::Types::CapeBoolean transpose=false)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    void mulvec(const Common::Types::CapeArrayDouble& vec,
                Common::Types::CapeArrayDouble& res,
                Common::Types::CapeBoolean transpose=false)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    void mulvec(int nvec, const double* vec, int nres, double* res,
                Common::Types::CapeBoolean transpose=false)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);
  private:
    void init_matrix(const DianaSparseArray& spar)
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds, Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

    void lusymb()
      throw (Common::Error::ECapeUnknown, Common::Error::ECapeBadArgument,
             Common::Error::ECapeOutOfBounds,Common::Error::ECapeSolvingError,
             Common::Error::ECapeNoMemory);

  private:
    Common::Types::CapeBoolean                       bNewStructure;
    Common::Types::CapeBoolean                       bNewValues;
    DianaSparseArray::crc_type crc32sa;

    void*               Symbolic;
    void*               Numeric;
    long                nRows;
    long                nCols;
    long                NZ;
    std::vector<double> Info;
    std::vector<double> Control;
    std::vector<double> WorkD;
    std::vector<double> Ax;
    std::vector<long>   WorkI;
    std::vector<long>   Ti;
    std::vector<long>   Tj;
    std::vector<long>   Ap;
    std::vector<long>   Ai;
    std::vector<long>   Mapping;
  };
};

#endif // DIANA_SPARSE_SOLVER_HPP