#include "../ezmtl/Matrix.h"
#ifdef _WIN32
#include <io.h> // for access()
#include <direct.h> // for mkdir()
#else
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#define _access access
#define _mkdir mkdir
#endif

template <typename T> int checkEigenVectors(Matrix<T>& A,Matrix<T>& EVec,Matrix<T>& EVal);
static int CreateDir(const char *dirName);
int main(){

	int isOK;
	{
		cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n";
		cout << "CAUTION: Microsoft Visual C++ may produce different answers \n";
		cout << "         between the Release and Debug modes\n";
		cout << "         when the single precision is used.\n";
		cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n";
		cout << "\n";
	}

#if 1

#if 1
	{
		Matrix<float> A;
		Vector<float> x;		x.Resize(2);		Vector<float> y(3);		x(1)=32; x(2)=15;		y(1)=3; y(2)=5; y(3)=8;
		Matrix<float> B(y,2,2);
		cout << "-----------------------------------------------------\n";
		cout << "General matrix operations\n";
		cout << "-----------------------------------------------------\n";
		B.Print(cout,"B");
		float tt=(float)Ipow(10,-5);	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		A.Resize(2,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=3; A(2,2)=5; A(2,3)=1;
		cout << "-----------------------------------------------------\n";
		cout << "Testing reshape" << endl;
		cout << "-----------------------------------------------------\n";
		A.Print(cout,"A");
		Reshape(A,3,2).Print(cout,"Reshape(A,3,2)");
	}
#endif
		
#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		A.Resize(2,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=3; A(2,2)=5; A(2,3)=1;
		x.Resize(2);
		Vector<float> y(3);
		x(1)=32; x(2)=15;
		y(1)=3; y(2)=5; y(3)=8;
		cout << "-----------------------------------------------------\n";		cout << "Testing x\'*A*y" << endl;
		cout << "-----------------------------------------------------\n";		float z=(float)MultiplyXtAY(x,A,y);		x.Print(cout,"  x");		A.Print(cout,"  A");		y.Print(cout,"  y");		cout << "  x\'*A*y=" << z << endl << endl;
	}
#endif
#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		A.Resize(3,3);
		A(1,1)=3; A(1,2)=2; A(1,3)=1; 
		A(2,1)=2; A(2,2)=4; A(2,3)=2; 
		A(3,1)=1; A(3,2)=2; A(3,3)=5; 
		Matrix<float> L;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Cholesky decomposition for symmetric positive definite matrices: A=L*L\'" << endl;
		cout << "-----------------------------------------------------\n";
		A.Print(cout,"  A");
		Chol(A,L);
		L.Print(cout,"  L");
		(A-L*Transpose(L)).Print(cout,"  A-L*L\'=0");
		Vector<float> b(3);
		b(1)=1; b(2)=-2; b(3)=5;
		CholSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		A.Resize(3,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=4; A(2,2)=5; A(2,3)=6; 
		A(3,1)=7; A(3,2)=8; A(3,3)=0; 
		
		cout << "-----------------------------------------------------\n";
		cout << "Testing LU decomposition for square matrices: A=L*U" << endl;
		cout << "-----------------------------------------------------\n";
		A.Print(cout,"A");

		Matrix<float> L,U,P;
		// A=L*U;
		cout << "  A=L*U" << endl;
		Lu(A,L,U);
		L.Print(cout,"  L");
		U.Print(cout,"  U");
		(A-L*U).Print(cout,"  A-L*U=0");
		
		// L*U=P*A
		cout << "  P*A=L*U" << endl;
		Lu(A,L,U,P);
		L.Print(cout,"  L");
		U.Print(cout,"  U");
		P.Print(cout,"  P");
		(L*U-P*A).Print(cout,"  L*U-P*A=0");

		Vector<float> b(3);
		b(1)=1; b(2)=-2; b(3)=5;
		x.Resize(3);
		LuSolve(A,b,x);
		x.Print(cout,"  x");
	
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Determinant" << endl;
		cout << "-----------------------------------------------------\n";
		A.Resize(3,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=4; A(2,2)=5; A(2,3)=6; 
		A(3,1)=7; A(3,2)=8; A(3,3)=0;
		A.Print(cout,"A");
		cout << "  det(A)=" << Det(A) << endl << endl;
		Matrix<float> C(2,2);
		C(1,1)=1; C(1,2)=2; 
		C(2,1)=3; C(2,2)=4;
		C.Print(cout,"C");
		cout << "  det(C)=" << Det(C) << endl << endl;
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Inverse" << endl;
		cout << "-----------------------------------------------------\n";
		A.Resize(3,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=4; A(2,2)=5; A(2,3)=6; 
		A(3,1)=7; A(3,2)=8; A(3,3)=0; 
		int retCode;
		A.Print(cout,"A");
		Inv(A,&retCode).Print(cout,"Inv(A)");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		Vector<float> b;
		A.Resize(3,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=3; A(2,2)=5; A(2,3)=1; 
		A(3,1)=1; A(3,2)=-1; A(3,3)=0; 
		cout << "-----------------------------------------------------\n";
		cout << "Testing QR decomposition for any matrices: A=Q*R" << endl;
		cout << "-----------------------------------------------------\n";
		Matrix<float> Q,R;
		isOK=Qr(A,Q,R);
		A.Print(cout,"  A");
		Q.Print(cout,"  Q");
		R.Print(cout,"  R");
		(A-Q*R).Print(cout,"  A-QR=0");
		b.Resize(3);
		b(1)=1; b(2)=2; b(3)=3;
		QrSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");

		A.Resize(2,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3;
		A(2,1)=4; A(2,2)=5; A(2,3)=6;
		isOK=Qr(A,Q,R);
		A.Print(cout,"  A");
		Q.Print(cout,"  Q");
		R.Print(cout,"  R");
		(A-Q*R).Print(cout,"  A-QR=0");
		b.Resize(2);
		b(1)=1; b(2)=2;
		QrSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");

		A.Resize(3,2);
		A(1,1)=1; A(1,2)=2;
		A(2,1)=3; A(2,2)=4;
		A(3,1)=5; A(3,2)=6;
		isOK=Qr(A,Q,R);
		A.Print(cout,"  A");
		Q.Print(cout,"  Q");
		R.Print(cout,"  R");
		(A-Q*R).Print(cout,"  A-QR=0");
		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		QrSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing LDL factorization: A=L*D*L\'" << endl;
		cout << "-----------------------------------------------------\n";
		Matrix<float> L(3,3),D(3,3);
		A.Resize(3,3);
		A(1,1)=3; A(1,2)=2; A(1,3)=1; 
		A(2,1)=2; A(2,2)=4; A(2,3)=2; 
		A(3,1)=1; A(3,2)=2; A(3,3)=5;

		A.Print(cout,"  A");
		Ldl(A,L,D);
		L.Print(cout,"  L");
		D.Print(cout,"  D");
		(A-L*D*Transpose(L)).Print(cout,"  A-L*D*L\'=0");

		Vector<float> b;
		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		LdlSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing QRCP factorization: " << endl;
		cout << "-----------------------------------------------------\n";

		Matrix<float> L(3,3),D(3,3);
		A.Resize(3,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3;
		A(2,1)=4; A(2,2)=5; A(2,3)=6;
		A(3,1)=7; A(3,2)=8; A(3,3)=0;

		Matrix<float> QRCP;
		Vector<float> diag;
		Vector<int> px;

		A.Print(cout,"  A");
		Qrcp(A,QRCP,diag,px);
		QRCP.Print(cout,"  QRCP");
		diag.Print(cout,"  diag");
		px.Print(cout,"  px");

		Vector<float> b;
		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		QrcpSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing BKP factorization" << endl;
		cout << "-----------------------------------------------------\n";
		Matrix<float> BKP(3,3);
		Vector<int> pivot, blocks;
		A.Resize(3,3);
		A(1,1)=3; A(1,2)=2; A(1,3)=1; 
		A(2,1)=2; A(2,2)=4; A(2,3)=2; 
		A(3,1)=1; A(3,2)=2; A(3,3)=5;

		A.Print(cout,"  A");
		Bkp(A,BKP,pivot,blocks);
		BKP.Print(cout,"  BKP");
		pivot.Print(cout,"  pivot");
		blocks.Print(cout,"  blocks");

		Vector<float> b;
		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		BkpSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		Matrix<float> EVal, EVec;
		float z;
		Matrix<float> ans;
		int i;

		cout << "-----------------------------------------------------\n";
		cout << "Testing Eigen analysis by the Jacobi transformation" << endl;
		cout << "-----------------------------------------------------\n";		A.Resize(2,2);
		A(1,1)=5; A(1,2)=1;
		A(2,1)=1; A(2,2)=5;
		A.Print(cout,"  A");
		Eig(A, EVec, EVal);
		EVal.Print(cout,"  EVal");
		EVec.Print(cout,"  EVec");
		for(i=1;i<=2;i++){
			Vector<float> x=EVec.Col(i);
			z=(float)Multiply3(Transpose(x),A,x);
			cout << "  EVec(" << i << ")\'*A*Evec(" << i << ")=" << z << endl;
		}
		ans=A*EVec-EVec*EVal;
		ans.Print(cout,"  A*EVec-EVec*EVal=0");
		cout << endl;
		
		A.Resize(3,3);		A(1,1)=3; A(1,2)=2; A(1,3)=1; 		A(2,1)=2; A(2,2)=4; A(2,3)=2; 		A(3,1)=1; A(3,2)=2; A(3,3)=5; 		
		cout << "-----------------------------------------------------\n";		cout << "Testing Eigen analysis by the Householder transformation and QL algorithm" << endl;		cout << "-----------------------------------------------------\n";
		A.Print(cout,"  A");		EigHHolder(A, EVec, EVal);		EVal.Print(cout,"  EVal");		EVec.Print(cout,"  EVec");		for(i=1;i<=3;i++){			Vector<float> x=EVec.Col(i);			z=(float)Multiply3(Transpose(x),A,x);			cout << "  EVec(" << i << ")\'*A*Evec(" << i << ")=" << z << endl;		}
		ans=A*EVec-EVec*EVal;		ans.Print(cout,"  A*EVec-EVec*EVal=0");
		cout << endl;
	}
#endif	
#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Simultaneous diagonalization: E\'*A*E=I, E\'*B*E=D" << endl;
		cout << "-----------------------------------------------------\n";		Matrix<float> E(3,3),D(3,3);
		A.Resize(3,3);
		A(1,1)=5; A(1,2)=2; A(1,3)=1; 
		A(2,1)=2; A(2,2)=5; A(2,3)=1; 
		A(3,1)=1; A(3,2)=1; A(3,3)=2; 
		B.Resize(3,3);
		B(1,1)=3; B(1,2)=2; B(1,3)=1; 
		B(2,1)=2; B(2,2)=2; B(2,3)=1; 
		B(3,1)=1; B(3,2)=1; B(3,3)=3; 
		
		A.Print(cout,"  A");
		B.Print(cout,"  B");		SimDiag(A, B, E, D);
		E.Print(cout,"  E");
		D.Print(cout,"  D");
		(Transpose(E)*A*E).Print(cout,"  E\'*A*E=I");
		(Transpose(E)*B*E).Print(cout,"  E\'*B*E=D");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Singular value decomposition: A=U*S*V\'" << endl;
		cout << "-----------------------------------------------------\n";
		Matrix<float> U, S, V;

		A.Resize(2,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3; 
		A(2,1)=4; A(2,2)=5; A(2,3)=6;

		Svd(A,U,S,V);
		A.Print(cout,"  A");
		U.Print(cout,"  U");
		S.Print(cout,"  S");
		V.Print(cout,"  V");
		(A-U*S*Transpose(V)).Print(cout, "  A-U*S*V\'=0");

		A.Resize(3,2);
		A(1,1)=1; A(1,2)=2;
		A(2,1)=3; A(2,2)=4;