DianaHslSolver.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: DianaHslSolver.hpp 8160 2009-03-16 18:36:30Z miha $
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef DIANA_HSL_SOLVER_HPP
#define DIANA_HSL_SOLVER_HPP

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

namespace Diana {

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

    DianaHslSolver& 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;

    fortrani               nRows;
    fortrani               nCols;
    fortrani               Ne;
    fortrani               JobA;
    fortrani               JobB;
    fortrani               JobC;
    fortrani               La;
    fortrand               Rinfo;
    std::vector<fortrand>  A;
    std::vector<fortrani>  Irn;
    std::vector<fortrani>  Jcn;
    std::vector<fortrand>  Ta;
    std::vector<fortrani>  Tirn;
    std::vector<fortrani>  Tjcn;
    std::vector<fortrani>  Keep;
    std::vector<fortrani>  Iw;
    std::vector<fortrani>  Info; //[MA48Info_Size];
    std::vector<fortrand>  W;
    std::vector<fortrand>  Error; //[MA48Error_Size];
    std::vector<fortrand>  Cntl; // [MA48Cntl_Size];
    std::vector<fortrani>  Icntl; // [MA48Icntl_Size]; 
    std::vector<fortrand>  B;
  };
};

#endif // DIANA_HSL_SOLVER_HPP