/*!
 * \file
 * \brief Matrices in decomposed forms (LDL', LU, UDU', etc).
 * \author Vaclav Smidl.
 *
 * -----------------------------------
 * BDM++ - C++ library for Bayesian Decision Making under Uncertainty
 *
 * Using IT++ for numerical operations
 * -----------------------------------
 */

#ifndef DC_H
#define DC_H

#include <itpp/itbase.h>

using namespace itpp;

/*! \brief Virtual class for representation of double symmetric matrices in square-root form.

All operations defined on this class should be optimized for the chosed decomposition.
*/
class sqmat {
public:
	/*!
	 * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$.
	 * @param  v Vector forming the outer product to be added
	 * @param  w weight of updating; can be negative

	 BLAS-2b operation.
	 */
	virtual void opupdt( const vec &v, double w ) =0;

	/*! \brief Conversion to full matrix.
	*/

	virtual mat to_mat() =0;

	/*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$
	@param C multiplying matrix,
	@param trans if true, product $V = C'*V*C$ will be computed instead;
	*/
	virtual void mult_sym( const mat &C, bool trans=true ) =0;

	
	/*!
	\brief Logarithm of a determinant.
	
	*/
	virtual double logdet() =0;

	/*!
	\brief Evaluates quadratic form $x= v'*V*v$;
	
	*/
	virtual double qform(vec &v) =0;

//	//! easy version of the 
//	sqmat inv();

	friend std::ostream &operator<< ( std::ostream &os, sqmat &sq );

	//! Clearing matrix so that it corresponds to zeros.
	virtual void clear() =0;
	
	//! Reimplementing common functions of mat: cols().
	virtual int cols() =0;

	//! Reimplementing common functions of mat: cols().
	virtual int rows() =0;

};


/*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix.

This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
*/
class fsqmat: sqmat {
	void opupdt( const vec &v, double w );
	mat to_mat();
	void mult_sym( const mat &C, bool trans=false );
   void mult_sym( const mat &C, fsqmat &U, bool trans=false );
	void inv(fsqmat &Inv);
	void clear();
	
	//! Constructor
	fsqmat(const mat &M);
	
	/*! \brief Matrix inversion preserving the chosen form.

	@param Inv a space where the inverse is stored.
	
	*/
	virtual void inv(fsqmat* Inv);
};

class ldmat: sqmat {
public:

	//! Construct by copy of L and D.
	ldmat( const mat &L, const vec &D );
	//! Construct by decomposition of full matrix V.
	ldmat( mat V );
	ldmat ();

	// Reimplementation of compulsory operatios

	void opupdt( const vec &v, double w );
	mat to_mat();
   void mult_sym( const mat &C, bool trans=false );
	void add ( const ldmat &ld2, double w=1.0 );
	double logdet();
	double qform(vec &v);
//	sqmat& operator -= ( const sqmat & ld2 );
	void clear();
	int cols();
	int rows();

	/*! \brief Matrix inversion preserving the chosen form.

	@param Inv a space where the inverse is stored.
	
	*/
	virtual void inv(ldmat &Inv);
	
	/*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class.

	@param Inv a space where the inverse is stored.
	
	*/
   void mult_sym( const mat &C, ldmat &U, bool trans=false );

	ldmat& operator += (const ldmat &ldA);
	ldmat& operator -= (const ldmat &ldA);
	ldmat& operator *= (double x);
	
protected:
	vec D;
	mat L;

};

//////// Operations:

inline ldmat& ldmat::operator += (const ldmat &ldA)  {this->add(ldA);return *this;}
inline ldmat& ldmat::operator -= (const ldmat &ldA)  {this->add(ldA,-1.0);return *this;}
inline int ldmat::cols(){return L.cols();}
inline int ldmat::rows(){return L.rows();}

#endif // DC_H