JMatrix2S.hh

Go to the documentation of this file.

#ifndef __JMATRIX2S__
#define __JMATRIX2S__
 
#include <cmath>
#include <algorithm>
 
#include "JMath/JMatrix2D.hh"
#include "JLang/JException.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JMATH {}
namespace JPP { using namespace JMATH; }
 
namespace JMATH {
 
  using JLANG::JDivisionByZero;
 
 
  /**
   * 2 x 2 symmetric matrix
   */
  class JMatrix2S :
    public JMatrix2D
  {
  public:
    /**
     * Default constructor.
     */
    JMatrix2S() :
      JMatrix2D()
    {}
 
 
    /**
     * Contructor.
     *
     * \param   A      matrix
     */
    JMatrix2S(const JMatrix2D& A) :
      JMatrix2D(A)
    {}
 
 
    /**
     * Contructor.
     * The lower triangle is internally set.
     *
     * \param   __a00   (0,0)
     * \param   __a01   (0,1)
     * \param   __a11   (1,1)
     */
    JMatrix2S(const double __a00, const double __a01,
              /*               */ const double __a11) :
      JMatrix2D(__a00, __a01,
                __a01, __a11)
    {}
 
 
    /**
     * Invert matrix
     */
    void invert()
    {
      using std::swap;
 
      // LDU decomposition
 
      int    p0 = 0; // permute row 0
 
      double val;
 
      val = a00;
 
      if (fabs(a10) > fabs(val)) {
      
        p0  = 1;
        val = a10;
 
        swap(a00,a10);
        swap(a01,a11);
      }
 
      if (val == 0) {
        throw JDivisionByZero("LDU decomposition zero pivot");
      }
 
      a10 /= val;
      a11 -= a10 * a01;
 
      // invert D
 
      if (a11 == 0) {
        throw JDivisionByZero("D matrix not invertable");
      }
 
      a00  = 1.0 / a00;
      a11  = 1.0 / a11;
 
      // invert U
 
      a01 *= -a00;
      a01 *=  a11;
 
      // invert L
 
      a10  = -a10;
 
      // U^-1 x L^-1
 
      a00 += a01 * a10;
      a10 *= a11;
 
      switch (p0) {
 
      case 1:
        swap(a00,a01);
        swap(a10,a11);
        break;
      }
    }
  };
}
 
#endif