lowerTriangularMatrix.cc

Go to the documentation of this file.

#include "lowerTriangularMatrix.h"
#include "matrixInversionExceptions.h"

using namespace oBLAS;


lowerTriangularMatrix::lowerTriangularMatrix(const int i) :
  ublas::triangular_matrix<double, ublas::lower>(i, i)
{ }


lowerTriangularMatrix::lowerTriangularMatrix(const ublas::matrix<double>& r)
  : ublas::triangular_matrix<double, ublas::lower>(r)
{ }


lowerTriangularMatrix::~lowerTriangularMatrix()
{
  delete fInverse;
}


void
lowerTriangularMatrix::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)(i, j) = r(i, j);
}


/* calculate and return inverse using recursive method from
   V. Blobel and E. Lohrmann "Statistische und numerische
   Methoden der Datenanalyse", Teubner Verlag, 1998 */

const lowerTriangularMatrix*
lowerTriangularMatrix::GetInverse()
  const
{
  const unsigned int size = this->size1();
  if (!fInverse)
    fInverse = new lowerTriangularMatrix(size);

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

  // some abbreviatons ...

  const lowerTriangularMatrix& t = *this;
  lowerTriangularMatrix& 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(i, j) * r(j, k);
      }
      r(i, k) = -sum / diag;
    }
  }

  return fInverse;
}