LinearAlgebra.cc

Go to the documentation of this file.

#include "LinearAlgebra.h"

#include <iostream>
#include <string>
#include <iomanip>
#include <sstream>

#include <cmath>
#include <limits>


namespace FdEnergyDepositFinderKG {

  void
  LowerTriangularMatrix::Zero()
  {
    fill(fElements.begin(), fElements.end(), 0);
  }


  unsigned int
  LowerTriangularMatrix::GetSize()
    const
  {
    return 0.5*(sqrt(1 + 8*fElements.size()) - 1);
  }


  ColumnVector
  LowerTriangularMatrix::operator*(const ColumnVector& rhs)
    const
  {
    const unsigned int n = GetSize();
    //std::cout << "Matrix size: " << n << "Column vector size: " << rhs.GetSize() << std::endl;

    if (!n) { // hack for "only direct light" case to avoit too many ifs in CalculateProfile
      ColumnVector zero(rhs.GetSize());
      zero.Zero();
      return zero;
    }
    if (n != rhs.GetSize())
      throw utl::OutOfBoundException("operator* size mismatch");
    ColumnVector retval(n);

    for (unsigned int row = 0; row < n; ++row) {
      const double* const rowArray = &fElements[(row*(row+1))/2];
      double sum = 0;
      for (unsigned int col = 0; col <= row; ++col)
        sum += rowArray[col] * rhs.fElements[col];
      retval.fElements[row] = sum;
    }

    return retval;
  }


  LowerTriangularMatrix
  LowerTriangularMatrix::operator+(const LowerTriangularMatrix& rhs)
    const
  {
    const unsigned int n = GetSize();
    if (n != rhs.GetSize())
      throw utl::OutOfBoundException("operator+ size mismatch");
    LowerTriangularMatrix retval(n);

    for (unsigned int row = 0; row < n; ++row) {
      const unsigned int idx_row = (row*(row+1))/2;
      for (unsigned int col = 0; col <= row; ++col) {
        const unsigned int idx = idx_row + col;
        retval.fElements[idx] = fElements[idx] + rhs.fElements[idx];
      }
    }
    return retval;
  }


  LowerTriangularMatrix
  LowerTriangularMatrix::GetInverse()
    const
  {
    const unsigned int n = GetSize();

    // some abbreviatons ...
    const LowerTriangularMatrix& t = *this;
    LowerTriangularMatrix r(n);

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

    // first calculate diagonal elements ...
    for (unsigned int i = 0; i < n; ++i) {
      const double tmp = t(i, i);
      if (fabs(tmp) > epsilon)
        r(i, i) = 1/tmp;
      else {
        std::ostringstream errMsg;
        errMsg << "singular lower triangular matrix "
                  "diagonal element " << i << " too small: "
               << tmp;
        throw SingularMatrixException(errMsg.str());
      }
    }

    // ... then off diagonal
    for (unsigned int i = 1; i < n; ++i) {
      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 r;
  }


  void
  LowerTriangularMatrix::Print()
    const
  {
    std::ostringstream dump;
    dump << std::resetiosflags(std::ios::scientific)
         << std::setiosflags(std::ios::fixed) << std::setprecision(2)
         << std::setfill(' ');
    const unsigned int n = GetSize();
    dump << "LowerTriangularMatrix(" << n << "):\n";
    for (unsigned int row = 0; row < n; ++row) {
      dump << "[ ";
      for (unsigned int col = 0; col < n; ++col) {
        dump << std::setw(8);
        if (row >= col)
          dump << (*this)(row, col) << " ";
        else
          dump << 0 << " ";
      }
      dump << " ]\n";
    }
    std::cout << dump.str() << std::flush;
  }


  LowerTriangularMatrix
  LowerTriangularMatrix::operator+(const DiagonalMatrix& rhs)
    const
  {
    const unsigned int n = GetSize();
    if (n != rhs.GetSize())
      throw utl::OutOfBoundException("operator+ size mismatch");
    LowerTriangularMatrix retval(*this);

    for (unsigned int row = 0; row < n; ++row) {
      const unsigned int idx = (row*(row+1))/2 + row;

      retval.fElements[idx] += rhs.fElements[row];
    }
    return retval;
  }


  /**************************************************************/

  const double&
  DiagonalMatrix::operator()(const int row, const int col)
    const
  {
    if (row != col)
      throw utl::OutOfBoundException("operator(i, j): i != j");
    return fElements[row];
  }


  double&
  DiagonalMatrix::operator()(const int row, const int col)
  {
    if (row != col)
      throw utl::OutOfBoundException("operator(i, j): i != j");
    return fElements[row];
  }


  void
  DiagonalMatrix::Zero()
  {
    fill(fElements.begin(), fElements.end(), 0);
  }


  void
  DiagonalMatrix::Print()
    const
  {
    std::ostringstream dump;
    dump << std::resetiosflags(std::ios::scientific)
         << std::setiosflags(std::ios::fixed) << std::setprecision(2)
         << std::setfill(' ');
    const unsigned int n = GetSize();
    dump << "DiagonalMatrix(" << n << "):\n";
    for (unsigned int row = 0; row < n; ++row) {
      dump << "[ ";
      for (unsigned int col = 0; col < n; ++col) {
        dump << std::setw(8);
        if (row == col)
          dump << (*this)(row) << ' ';
        else
          dump << "0 ";
      }
      dump << " ]\n";
    }
    std::cout << dump.str() << std::flush;
  }


  DiagonalMatrix
  DiagonalMatrix::operator+(const DiagonalMatrix& rhs)
    const
  {
    const unsigned int n = GetSize();
    if (n != rhs.GetSize())
      throw utl::OutOfBoundException("operator+ size mismatch");
    DiagonalMatrix retval(n);

    for (unsigned int row = 0; row < n; ++row)
      retval.fElements[row] = fElements[row] + rhs.fElements[row];

    return retval;
  }


  ColumnVector
  DiagonalMatrix::operator*(const ColumnVector& rhs)
    const
  {
    const unsigned int n = GetSize();
    if (n != rhs.GetSize())
      throw utl::OutOfBoundException("operator* size mismatch");
    ColumnVector retval(n);

    for (unsigned int row = 0; row < n; ++row)
      retval.fElements[row] = fElements[row] * rhs.fElements[row];

    return retval;
  }


  DiagonalMatrix
  DiagonalMatrix::GetInverse()
    const
  {
    const unsigned int n = GetSize();

    // some abbreviatons ...
    const DiagonalMatrix& t = *this;
    DiagonalMatrix r(n);

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

    // calculate diagonal elements ...
    for (unsigned int i = 0; i < n; ++i) {
      const double tmp = t(i, i);
      if (fabs(tmp) > epsilon)
        r(i, i) = 1/tmp;
      else {
        std::ostringstream errMsg;
        errMsg << "singular diagonal matrix "
                  "element " << i << " too small: "
               << tmp;
        throw SingularMatrixException(errMsg.str());
      }
    }

    return r;
  }


  /**************************************************************/

  ColumnVector
  ColumnVector::operator+(const ColumnVector& rhs)
    const
  {
    const unsigned int n = GetSize();
    if (n != rhs.GetSize())
      throw utl::OutOfBoundException("operator+ size mismatch");

    ColumnVector retval(n);
    for (unsigned int row = 0; row < n; ++row)
      retval.fElements[row] = fElements[row] + rhs.fElements[row];
    return retval;
  }


  void
  ColumnVector::Zero()
  {
    fill(fElements.begin(), fElements.end(), 0);
  }

}