DianaNormalEquationsSolver.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$
 * -----------------------------------------------------------------------------
 *  Description:
 */

#ifndef DIANA_NORMAL_EQUATIONS_SOLVER_HPP
#define DIANA_NORMAL_EQUATIONS_SOLVER_HPP

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

namespace Diana {

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

    DianaNormalEquationsSolver& 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);

    void free_cholmod();
    
  private:
    Common::Types::CapeBoolean bNewStructure;
    Common::Types::CapeBoolean bNewValues;
    DianaSparseArray::crc_type crc32sa;

    long                       nRows;
    long                       nCols;
    long                       nNoneZeros;

    void*                      pCM;    // cholmod_common
    void*                      cspA;   // cholmod_sparse
    void*                      cdeB;   // cholmod_dense
    void*                      cdeX;   // cholmod_dense
    void*                      cdeAtB; // cholmod_dense
    void*                      csfL;   // cholmod_factor
    Common::Types::CapeDouble  dblUround;
  };
};

#endif // DIANA_NORMAL_EQUATIONS_SOLVER_HPP