Context Navigation

square_mat.h @ 384

Revision 384, 8.4 kB (checked in by mido, 15 years ago)
possibly broken?
Property svn:eol-style set to `native`

Rev	Line
[7]	1	/*!
	2	* \file
	3	* \brief Matrices in decomposed forms (LDL', LU, UDU', etc).
	4	* \author Vaclav Smidl.
	5	*
	6	* -----------------------------------
	7	* BDM++ - C++ library for Bayesian Decision Making under Uncertainty
	8	*
	9	* Using IT++ for numerical operations
	10	* -----------------------------------
	11	*/
	12
	13	#ifndef DC_H
	14	#define DC_H
	15
[262]	16
[180]	17	#include "../itpp_ext.h"
[7]	18
[219]	19	/*!
	20	\defgroup math Auxiliary math functions
	21	@{
	22	*/
	23
[7]	24	using namespace itpp;
	25
[37]	26	//! Auxiliary function dydr; dyadic reduction
	27	void dydr( double * r, double f, double Dr, double Df, double R, int jl, int jh, double *kr, int m, int mx );
	28
	29	//! Auxiliary function ltuinv; inversion of a triangular matrix;
	30	//TODO can be done via: dtrtri.f from lapack
	31	mat ltuinv( const mat &L );
	32
[7]	33	/*! \brief Virtual class for representation of double symmetric matrices in square-root form.
	34
[37]	35	All operations defined on this class should be optimized for the chosen decomposition.
[7]	36	*/
[22]	37	class sqmat
	38	{
	39	public:
	40	/*!
[75]	41	* Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$.
[22]	42	* @param v Vector forming the outer product to be added
	43	* @param w weight of updating; can be negative
[7]	44
[22]	45	BLAS-2b operation.
	46	*/
	47	virtual void opupdt ( const vec &v, double w ) =0;
[7]	48
[22]	49	/*! \brief Conversion to full matrix.
	50	*/
[7]	51
[168]	52	virtual mat to_mat() const =0;
[7]	53
[75]	54	/! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = CV*C'\f$
[22]	55	@param C multiplying matrix,
	56	*/
[26]	57	virtual void mult_sym ( const mat &C ) =0;
	58
[75]	59	/! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'V*C\f$
[26]	60	@param C multiplying matrix,
	61	*/
	62	virtual void mult_sym_t ( const mat &C ) =0;
[7]	63
	64
[22]	65	/*!
	66	\brief Logarithm of a determinant.
[7]	67
[22]	68	*/
[26]	69	virtual double logdet() const =0;
[22]	70
	71	/*!
[75]	72	\brief Multiplies square root of \f$V\f$ by vector \f$x\f$.
[22]	73
	74	Used e.g. in generating normal samples.
	75	*/
[32]	76	virtual vec sqrt_mult (const vec &v ) const =0;
[22]	77
	78	/*!
[75]	79	\brief Evaluates quadratic form \f$x= v'Vv\f$;
[22]	80
	81	*/
[32]	82	virtual double qform (const vec &v ) const =0;
[22]	83
[75]	84	/*!
	85	\brief Evaluates quadratic form \f$x= v'inv(V)v\f$;
	86
	87	*/
	88	virtual double invqform (const vec &v ) const =0;
	89
[22]	90	// //! easy version of the
[7]	91	// sqmat inv();
	92
[22]	93	//! Clearing matrix so that it corresponds to zeros.
	94	virtual void clear() =0;
[7]	95
[22]	96	//! Reimplementing common functions of mat: cols().
	97	int cols() const {return dim;};
[7]	98
[22]	99	//! Reimplementing common functions of mat: cols().
	100	int rows() const {return dim;};
	101
[28]	102	//! Destructor for future use;
	103	virtual ~sqmat(){};
[32]	104	//! Default constructor
	105	sqmat(const int dim0): dim(dim0){};
[270]	106	//! Default constructor
	107	sqmat(): dim(0){};
[22]	108	protected:
[33]	109	//! dimension of the square matrix
[22]	110	int dim;
[7]	111	};
	112
	113
	114	/*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix.
	115
	116	This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
	117	*/
[22]	118	class fsqmat: public sqmat
	119	{
	120	protected:
[33]	121	//! Full matrix on which the operations are performed
[22]	122	mat M;
	123	public:
	124	void opupdt ( const vec &v, double w );
[168]	125	mat to_mat() const;
[26]	126	void mult_sym ( const mat &C);
	127	void mult_sym_t ( const mat &C);
[75]	128	//! store result of \c mult_sym in external matrix \f$U\f$
[32]	129	void mult_sym ( const mat &C, fsqmat &U) const;
[75]	130	//! store result of \c mult_sym_t in external matrix \f$U\f$
[32]	131	void mult_sym_t ( const mat &C, fsqmat &U) const;
[22]	132	void clear();
[7]	133
[26]	134	//! Default initialization
[270]	135	fsqmat(){}; // mat will be initialized OK
[26]	136	//! Default initialization with proper size
	137	fsqmat(const int dim0); // mat will be initialized OK
[22]	138	//! Constructor
	139	fsqmat ( const mat &M );
[37]	140	//! Constructor
[183]	141	fsqmat ( const fsqmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");};
	142	//! Constructor
[37]	143	fsqmat ( const vec &d ):sqmat(d.length()){M=diag(d);};
[22]	144
[28]	145	//! Destructor for future use;
	146	virtual ~fsqmat(){};
	147
	148
[22]	149	/*! \brief Matrix inversion preserving the chosen form.
	150
	151	@param Inv a space where the inverse is stored.
	152
	153	*/
	154	virtual void inv ( fsqmat &Inv );
	155
[26]	156	double logdet() const {return log ( det ( M ) );};
[32]	157	double qform (const vec &v ) const {return ( v* ( M*v ) );};
[75]	158	double invqform (const vec &v ) const {return ( v* ( itpp::inv(M)*v ) );};
[37]	159	vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;};
[22]	160
[85]	161	//! Add another matrix in fsq form with weight w
	162	void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;};
	163
[37]	164	//! Access functions
	165	void setD (const vec &nD){M=diag(nD);}
	166	//! Access functions
[51]	167	vec getD (){return diag(M);}
	168	//! Access functions
[37]	169	void setD (const vec &nD, int i){for(int j=i;j<nD.length();j++){M(j,j)=nD(j-i);}} //Fixme can be more general
	170
[85]	171
[33]	172	//! add another fsqmat matrix
[22]	173	fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;};
[33]	174	//! subtrack another fsqmat matrix
[22]	175	fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;};
[33]	176	//! multiply by a scalar
[22]	177	fsqmat& operator = ( double x ) {M=x;return *this;};
	178	// fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
[33]	179	//! print full matrix
[32]	180	friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );
[26]	181
[7]	182	};
	183
[180]	184	/*! \brief Matrix stored in LD form, (commonly known as UD)
[33]	185
[75]	186	Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored.
[33]	187	All inplace operations modifies only these and the need to compose and decompose the matrix is avoided.
	188	*/
[22]	189	class ldmat: sqmat
	190	{
	191	public:
	192	//! Construct by copy of L and D.
	193	ldmat ( const mat &L, const vec &D );
	194	//! Construct by decomposition of full matrix V.
[26]	195	ldmat (const mat &V );
[180]	196	//! Construct by restructuring of V0 accordint to permutation vector perm.
	197	ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){ ldform(V0.L.get_cols(perm), V0.D);};
[22]	198	//! Construct diagonal matrix with diagonal D0
	199	ldmat ( vec D0 );
[26]	200	//!Default constructor
[22]	201	ldmat ();
[26]	202	//! Default initialization with proper size
	203	ldmat(const int dim0);
[180]	204
[28]	205	//! Destructor for future use;
	206	virtual ~ldmat(){};
	207
[22]	208	// Reimplementation of compulsory operatios
[7]	209
[22]	210	void opupdt ( const vec &v, double w );
[168]	211	mat to_mat() const;
[26]	212	void mult_sym ( const mat &C);
	213	void mult_sym_t ( const mat &C);
[33]	214	//! Add another matrix in LD form with weight w
[22]	215	void add ( const ldmat &ld2, double w=1.0 );
[26]	216	double logdet() const;
[32]	217	double qform (const vec &v ) const;
[75]	218	double invqform (const vec &v ) const;
[22]	219	void clear();
[26]	220	int cols() const;
	221	int rows() const;
[32]	222	vec sqrt_mult ( const vec &v ) const;
[7]	223
[180]	224
[22]	225	/*! \brief Matrix inversion preserving the chosen form.
	226	@param Inv a space where the inverse is stored.
	227	*/
[32]	228	virtual void inv ( ldmat &Inv ) const;
[8]	229
[75]	230	/*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class.
[33]	231	@param C matrix to multiply with
[26]	232	@param U a space where the inverse is stored.
[22]	233	*/
[32]	234	void mult_sym ( const mat &C, ldmat &U) const;
[12]	235
[75]	236	/*! \brief Symmetric multiplication of \f$U\f$ by a transpose of a general matrix \f$C\f$, result of which is stored in the current class.
[33]	237	@param C matrix to multiply with
[26]	238	@param U a space where the inverse is stored.
	239	*/
[32]	240	void mult_sym_t ( const mat &C, ldmat &U) const;
[26]	241
	242
[75]	243	/*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$
[26]	244
[75]	245	The new decomposition fullfills: \f$A'diag(D)A = self.L'diag(self.D)self.L\f$
[22]	246	@param A general matrix
	247	@param D0 general vector
	248	*/
[26]	249	void ldform (const mat &A,const vec &D0 );
[22]	250
[32]	251	//! Access functions
	252	void setD (const vec &nD){D=nD;}
[33]	253	//! Access functions
	254	void setD (const vec &nD, int i){D.replace_mid(i,nD);} //Fixme can be more general
	255	//! Access functions
[32]	256	void setL (const vec &nL){L=nL;}
	257
[98]	258	//! Access functions
	259	const vec& _D() const {return D;}
	260	//! Access functions
	261	const mat& _L() const {return L;}
	262
[33]	263	//! add another ldmat matrix
[22]	264	ldmat& operator += ( const ldmat &ldA );
[33]	265	//! subtract another ldmat matrix
[22]	266	ldmat& operator -= ( const ldmat &ldA );
[33]	267	//! multiply by a scalar
[22]	268	ldmat& operator *= ( double x );
	269
[33]	270	//! print both \c L and \c D
[32]	271	friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );
[22]	272
	273	protected:
[75]	274	//! Positive vector \f$D\f$
[22]	275	vec D;
[75]	276	//! Lower-triangular matrix \f$L\f$
[22]	277	mat L;
	278
[7]	279	};
	280
	281	//////// Operations:
[33]	282	//!mapping of add operation to operators
[22]	283	inline ldmat& ldmat::operator += ( const ldmat &ldA ) {this->add ( ldA );return *this;}
[33]	284	//!mapping of negative add operation to operators
[22]	285	inline ldmat& ldmat::operator -= ( const ldmat &ldA ) {this->add ( ldA,-1.0 );return *this;}
[33]	286	//!access function
[26]	287	inline int ldmat::cols() const {return dim;}
[33]	288	//!access function
[26]	289	inline int ldmat::rows() const {return dim;}
[7]	290
[219]	291	/! @} /
	292
[7]	293	#endif // DC_H

Note: See TracBrowser for help on using the browser.

Context Navigation

root/library/bdm/math/square_mat.h @ 384

Download in other formats: