LinearAlgebra.h

Go to the documentation of this file.

#ifndef _FdEnergyDepositFinderKG_LinearAlgebra_h_
#define _FdEnergyDepositFinderKG_LinearAlgebra_h_

#include <utl/AugerException.h>
#include <vector>


namespace FdEnergyDepositFinderKG {

  class DiagonalMatrix;
  class LowerTriangularMatrix;
  class ColumnVector;


  class ColumnVector {
  public:
    ColumnVector(const int n) : fElements(n) { }

    ColumnVector(const int n, const double x)
      : fElements(n, x) { }

    ColumnVector() { }

    const double& operator()(const int row) const
    { return fElements[row]; }
    double& operator()(const int row) { return fElements[row]; }

    void Clear() { fElements.clear(); }
    void Zero();

    void SetSize(const int n) { fElements.resize(n); }
    unsigned int GetSize() const { return fElements.size(); }

    ColumnVector operator+(const ColumnVector& rhs) const;

  private:
    std::vector<double> fElements;

    friend class LowerTriangularMatrix;
    friend class DiagonalMatrix;
  };


  class LowerTriangularMatrix {
  public:
    LowerTriangularMatrix(const int n) : fElements(n*(n+1)/2) { }

    LowerTriangularMatrix(const int n, const double x)
      : fElements(n*(n+1)/2, x) { }

    LowerTriangularMatrix() { }

    // Note: for proper access asserts, we'd have to make this non-inlined...
    const double& operator()(const int row, const int col) const
    { return fElements[(row*(row+1))/2 + col]; }
    double& operator()(const int row, const int col)
    { return fElements[(row*(row+1))/2 + col]; }

    void Clear() { fElements.clear(); }
    void Zero();

    void SetSize(const int n) { fElements.resize(n*(n+1)/2); }
    unsigned int GetSize() const;

    LowerTriangularMatrix GetInverse() const;

    void Print() const;

    LowerTriangularMatrix operator+(const LowerTriangularMatrix& rhs) const;
    LowerTriangularMatrix operator+(const DiagonalMatrix& rhs) const;
    ColumnVector operator*(const ColumnVector& rhs) const;

  private:
    std::vector<double> fElements;

    friend class DiagonalMatrix;
  };


  class DiagonalMatrix {
  public:
    DiagonalMatrix(const int n) : fElements(n) { }

    DiagonalMatrix(const int n, const double x)
      : fElements(n, x) { }

    DiagonalMatrix() { }

    const double& operator()(const int row, const int col) const;
    double& operator()(const int row, const int col);

    const double& operator()(const int iDiagonal) const
    { return fElements[iDiagonal]; }
    double& operator()(const int iDiagonal)
    { return fElements[iDiagonal]; }

    void Clear() { fElements.clear(); }
    void Zero();

    void SetSize(const int n) { fElements.resize(n); }
    unsigned int GetSize() const { return fElements.size(); }

    void Print() const;

    DiagonalMatrix operator+(const DiagonalMatrix& rhs) const;
    ColumnVector operator*(const ColumnVector& rhs) const;

    DiagonalMatrix GetInverse() const;

  private:
    std::vector<double> fElements;

    friend class LowerTriangularMatrix;
  };


  class SingularMatrixException : public utl::AugerException {
  public:
    SingularMatrixException(const std::string& message = std::string())
      : AugerException(message) { }

    virtual std::string GetExceptionName() const
    { return "Singular Matrix Exception"; }
  };

}


#endif