#include "math/square_mat.h"

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main() {
	// Kalman filter
	mat A0 = randu(3,3);
	mat A = A0*A0.T();
	
	//Test constructor
	fsqmat Fsq(A);
   cout << "Testing constructors:" << endl
       << "A = " << A << endl
       << "Fsq.to_mat() = " << Fsq.to_mat() << endl << endl;
       
   //Test inversion
	fsqmat iFsq(3);
	Fsq.inv(iFsq);
   cout << "inv(A) = " << inv(A) <<endl
   	<< "inv(Fsq).to_mat() = " << iFsq.to_mat() <<endl <<endl;
   	
   //Test logdet
   cout << "logdet(A) = " << log(det(A)) << endl
			<< "logdet(Fsq) = " << Fsq.logdet() << endl <<endl;
				
   //Test add
   fsqmat Fsq2(Fsq);
   Fsq2.add(Fsq);
   cout << "A+A = " << A+A << endl
			<<	"Fsq2.add(Fsq) = " << Fsq2.to_mat() << endl <<endl;
				
	vec v=randu(3);
   //Test qform
   cout << "vAv' = " << v*(A*v) << endl
			<<	"qform(Fsq,v) = " << Fsq.qform(v) << endl <<endl;
	
   //Test invqform
   cout << "v inv(A)v' = " << v*(inv(A)*v) << endl
			<<	"invqform(Fsq,v) = " << Fsq.invqform(v) << endl <<endl;
				
   //Test opupdate
   Fsq2=Fsq;
   Fsq2.opupdt(v,1.0);
   cout << "A+vv' = " << A+outer_product(v,v) << endl
			<<	"opupdt(Fsq,v) = " << Fsq2.to_mat() << endl <<endl;
				
}