DianaLSQRSolver.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_LSQR_SOLVER_HPP
#define DIANA_LSQR_SOLVER_HPP

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

namespace Diana {

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

    DianaLSQRSolver& operator=(const DianaSparseArray& spar)
      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);

    double norm2(int n, const double* vec);

  private:
    DianaSparseArray::crc_type crc32sa;

    long                nRows;
    long                nCols;
    long                nNonZeros;
    std::vector<long>   irow;
    std::vector<long>   jcol;
    std::vector<double> aval;
    Common::Types::CapeDouble  dblUround;
  };
};

#endif // DIANA_LSQR_SOLVER_HPP