/***************************************************************
 * PrMatrix6.h
 *
 * This class provides a 6x6 matrix.
 *
 ****************************************************************/

/*
 * modification history
 *----------------------
 *
 * 11/20/97: K.C. Chang: added inline methods.
 * 11/05/97: K.C. Chang: created.
 */

#ifndef _PrMatrix6_h
#define _PrMatrix6_h

#include "PrGlobalDefn.h"
#include <assert.h>

#include "PrMatrix3.h"
#include "PrVector3.h"
#include "PrVector6.h"

class PrMatrix6
{
  public:
    PrMatrix6() {} // initialized to [0]
    PrMatrix6(const PrMatrix6 &m);
    PrMatrix6(const PrMatrix3 &m11,const PrMatrix3 &m12,
	      const PrMatrix3 &m21,const PrMatrix3 &m22);
    ~PrMatrix6() {}

    PrMatrix6 &operator=(const PrMatrix6 &m);
    PrMatrix3 &operator[](int i);
    const PrMatrix3 &operator[](int i) const;

    void zero();
    Float &elementAt(int i,int j);
    const Float &elementAt(int i,int j) const;

    PrMatrix6 operator-() const;
    void negate(PrMatrix6 &dest) const;

    PrMatrix6 operator~() const;
    void inverse(PrMatrix6 &dest) const;
    //(symmetric and positive definite)
    void inverseSPD(PrMatrix6 &dest) const;

    // y = (*this) x
    // solve for x given y
    void solve(const PrVector6 &y,PrVector6 &x) const;
    //(symmetric and positive definite)
    void solveSPD(const PrVector6 &y,PrVector6 &x) const;

    PrMatrix6 operator+(const PrMatrix6 &m) const;
    void add(const PrMatrix6 &m,PrMatrix6 &dest) const;

    PrMatrix6 operator-(const PrMatrix6 &m) const;
    void subtract(const PrMatrix6 &m,PrMatrix6 &dest) const;

    PrMatrix6 operator*(const PrMatrix6 &m) const;
    void multiply(const PrMatrix6 &m,PrMatrix6 &dest) const;

    PrVector6 operator*(const PrVector6 &v) const;
    void multiply(const PrVector6 &v,PrVector6 &dest) const;

    PrMatrix6 operator*(Float s) const;
    void multiply(Float s,PrMatrix6 &dest) const;

    PrMatrix6 &operator+=(const PrMatrix6 &m);
    PrMatrix6 &operator-=(const PrMatrix6 &m);
    PrMatrix6 &operator*=(const PrMatrix6 &m);
    PrMatrix6 &operator*=(Float s);

    PrMatrix6 transpose() const;
    void transpose(PrMatrix6 &dest) const;

    void diagonal(const PrVector3 &v1,const PrVector3 &v2); // set the diagonal
    PrVector6 diagonal() const;
    void diagonal(PrVector6 &dest) const;

    void display(char *name=NULL);

  private:
    PrMatrix3 mat3_[4];

    void LUdecomp(PrMatrix6 &lu) const;
    //(symmetric and positive definite)
    void LUdecompSPD(PrMatrix6 &lu) const;

    // y = LU x
    // find x given LU and y 
    void BackSub(const PrMatrix6 &lu,const PrVector6 &y,PrVector6 &x) const;
    //(symmetric and positive definite)
    void BackSubSPD(const PrMatrix6 &lu,const PrVector6 &y,PrVector6 &x) const;
};

// inline methods

inline PrMatrix6::PrMatrix6(const PrMatrix6 &m)
{
    mat3_[0] = m.mat3_[0];
    mat3_[1] = m.mat3_[1];
    mat3_[2] = m.mat3_[2];
    mat3_[3] = m.mat3_[3];
}

inline PrMatrix6::PrMatrix6(const PrMatrix3 &m11,const PrMatrix3 &m12,
			    const PrMatrix3 &m21,const PrMatrix3 &m22)
{
    mat3_[0] = m11;
    mat3_[1] = m12;
    mat3_[2] = m21;
    mat3_[3] = m22;
}

inline PrMatrix6 &PrMatrix6::operator=(const PrMatrix6 &m)
{
    if (this != &m)
      {
	  mat3_[0] = m.mat3_[0];  
	  mat3_[1] = m.mat3_[1];
	  mat3_[2] = m.mat3_[2];
	  mat3_[3] = m.mat3_[3];
      }
    return (*this);
}

inline void PrMatrix6::zero()
{
    mat3_[0].zero();
    mat3_[1].zero();
    mat3_[2].zero();
    mat3_[3].zero();
}

inline Float &PrMatrix6::elementAt(int i,int j)
{
    assert(i>=0 && i<6 && j>=0 && j<6);
    if (i<3)
      {
	if (j<3)
	  return mat3_[0][i][j];
	else
	  return mat3_[1][i][j-3];
      }
    else
      {
	if (j<3)
	  return mat3_[2][i-3][j];
	else
	  return mat3_[3][i-3][j-3];
      }
}

inline const Float &PrMatrix6::elementAt(int i,int j) const
{
    assert(i>=0 && i<6 && j>=0 && j<6);
    if (i<3)
      {
	if (j<3)
	  return mat3_[0][i][j];
	else
	  return mat3_[1][i][j-3];
      }
    else
      {
	if (j<3)
	  return mat3_[2][i-3][j];
	else
	  return mat3_[3][i-3][j-3];
      }
}

inline PrMatrix3 &PrMatrix6::operator[](int i)
{
    assert(i>=0 && i<4);
    return (mat3_[i]);
}

inline const PrMatrix3 &PrMatrix6::operator[](int i) const
{
    assert(i>=0 && i<4);
    return (mat3_[i]);
}


#endif // _PrMatrix6_h
