upperTriangularMatrix.cc

Go to the documentation of this file.

#include "upperTriangularMatrix.h"
#include "matrixInversionExceptions.h"

using namespace oBLAS;


/*****************************************************
 *
 *  constructor of upperTriangularMatrix
 *
 *****************************************************/
upperTriangularMatrix::upperTriangularMatrix(int i) :
  ublas::triangular_matrix<double, ublas::upper>(i, i)
{
  Init();
}


/*****************************************************
 *
 *  copy constructor from ublas::matrix
 *
 *****************************************************/
upperTriangularMatrix::upperTriangularMatrix(const ublas::matrix<double>& r) :
  ublas::triangular_matrix<double, ublas::upper>(r)
{
  Init();
}


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


/*****************************************************
 *
 *  destructor of upperTriangularMatrix
 *
 *****************************************************/
upperTriangularMatrix::~upperTriangularMatrix()
{
  delete fInverse;
}


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


/************************************************************
 *
 * calculate and return inverse using recursive method from
 * V. Blobel and E. Lohrmann "Statistische und numerische
 * Methoden der Datenanalyse", Teubner Verlag, 1998
 *
 ************************************************************/
const upperTriangularMatrix*
upperTriangularMatrix::GetInverse()
{
  const unsigned int size = this->size1();
  if (fInverse == NULL)
    fInverse = new upperTriangularMatrix(size);

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

  // some abbreviatons ...

  upperTriangularMatrix& t = *this;
  upperTriangularMatrix& r = *fInverse;

  // first calculate diagonal elements ...
  for (unsigned int i = 0; i < size; ++i) {
    const double tmp = t(i, i);
    if (fabs(tmp) > epsilon)
      r(i, i) = 1/tmp;
    else
      throw SingularMatrixException();
  }

  // ... then off diagonal
  for (unsigned int i = 1; i < size; ++i) {
    const double diag = t(i, i);
    for (unsigned int k = 0; k < i; ++k) {
      double sum = 0;
      for (unsigned int j = k; j < i; ++j) {
        sum += t(j, i) * r(k, j);
      }
      r(k, i) = -sum / diag;
    }
  }
  return fInverse;
}