diagonalMatrix.cc

Go to the documentation of this file.

#include "diagonalMatrix.h"
#include "matrixInversionExceptions.h"

using namespace oBLAS;


/*****************************************************
 *
 *  constructor of diagonalMatrix
 *
 *****************************************************/
diagonalMatrix::diagonalMatrix(int i) :
  ublas::banded_matrix<double>(i, i, 0, 0)
{
  Init();
  //std::cout << " diag1 " << this << " " << fInverse << std::endl;
}


/*****************************************************
 *
 *  copy constructor from ublas::matrix
 *
 *****************************************************/
diagonalMatrix::diagonalMatrix(const ublas::matrix<double>& r) :
  ublas::banded_matrix<double>(r.size1(), r.size1(), 0, 0)
{
  const unsigned int size = this->size1();
  for (unsigned int i = 0; i < size; ++i)
      (*this)(i, i) = r(i, i);
  Init();
  //std::cout << " diag2 " << this << " " << fInverse << std::endl;
}


/*****************************************************
 *
 *  initialize
 *
 *****************************************************/
void
diagonalMatrix::Init()
{
  fInverse = NULL;
}


/*****************************************************
 *
 *  destructor of diagonalMatrix
 *
 *****************************************************/
diagonalMatrix::~diagonalMatrix()
{
  //std::cout << " ~diag " << this << " " << fInverse << std::endl;
  delete fInverse;
}


/*****************************************************
 *
 *  assignment to boost matrix
 *
 *****************************************************/
void
diagonalMatrix::operator=(const ublas::matrix<double>& r)
{
  const unsigned int size = this->size1();
  for (unsigned int i = 0; i < size; ++i)
      (*this)(i, i) = r(i, i);
}


/*****************************************************
 *
 *  calculate and return inverse of matrix
 *
 *****************************************************/
const diagonalMatrix*
diagonalMatrix::GetInverse()
{
  //std::cout << "dinv " << this << " " << fInverse << std::endl;
  const unsigned int size = this->size1();
  if (fInverse == NULL)
    fInverse = new diagonalMatrix(size);

  const double epsilon = std::numeric_limits<double>::epsilon();

  for (unsigned int i = 0; i < size; ++i) {
    const double tmp = (*this)(i, i);
    if (fabs(tmp) > epsilon)
      (*fInverse)(i, i) = 1/tmp;
    else
      throw SingularMatrixException();
  }
  return fInverse;
}