A(3,1)=5; A(3,2)=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,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;

		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");

		int i,k,m,n;
		for(k=1;k<=3;k++){
			if(k==1)     { m=20; n=20; }
			else if(k==2){ m=30; n=10; }
			else if(k==3){ m=10; n=30; }
			cout << "  -----------------------------------------------------\n";
			cout << "  Testing with random matrices " << m << "x" << n << endl;
			cout << "  -----------------------------------------------------\n";
			for(i=1;i<=10;i++){
				A=Matrix<float>::Randn(m,n);
				isOK=Svd(A,U,S,V);
				if(isOK==0){
					cerr << "ERROR1 in SVD()\n";
				}
				B=U*S*Transpose(V);
				if((A-B).IsZero(1e-5)==0){
					cerr << "ERROR2 in SVD()\n";
					//(A-B).Print(cout,"A-U*S*V\'");
				}
				cout << "  Matrix " << i << " OK.\n";
			}
			
		}
		cout << "-----------------------------------------------------\n";
		cout << "Solving linear equations: A * x = b" << 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;
		Vector<float> b;
		b.Resize(3);
		b(1)=1; b(2)=2; b(3)=5;
		A.Print(cout,"  A");
		b.Print(cout,"  b");
		SvdSolve(A,b,x);
		x.Print(cout,"  x=A#b");

	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Pseudoinverse: Pinv(A)" << endl;
		cout << "-----------------------------------------------------\n";
		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; 
		A.Print(cout,"  A");
		Matrix<float> X;
		int retCode;
		X=Pinv(A,&retCode);
		(X).Print(cout,"  pinv(A)");
		(A*X).Print(cout,"  A*Pinv(A)");
		(X*A).Print(cout,"  Pinv(A)*A");
		(A*X*A).Print(cout,"  A*Pinv(A)*A");
		(X*A*X).Print(cout,"  Pinv(A)*A*Pinv(A)");
		Pinv(X,&retCode).Print(cout,"  Pinv(Pinv(A))");

		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;
		A.Print(cout,"  A");
		X=Pinv(A,&retCode);
		(X).Print(cout,"  pinv(A)");	
		(A*X).Print(cout,"  A*Pinv(A)");
		(X*A).Print(cout,"  Pinv(A)*A");
		(A*X*A).Print(cout,"  A*Pinv(A)*A");
		(X*A*X).Print(cout,"  Pinv(A)*A*Pinv(A)");
		Pinv(X,&retCode).Print(cout,"  Pinv(Pinv(A))");

		
		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");
		X=Pinv(A,&retCode);
		(X).Print(cout,"  pinv(A)");
		(A*X).Print(cout,"  A*Pinv(A)");
		(X*A).Print(cout,"  Pinv(A)*A");
		(A*X*A).Print(cout,"  A*Pinv(A)*A");
		(X*A*X).Print(cout,"  Pinv(A)*A*Pinv(A)");
		Pinv(X,&retCode).Print(cout,"  Pinv(Pinv(A))");

	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing solutions of linear equations: A*x=b" << endl;
		cout << "-----------------------------------------------------\n";
		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;
		A.Print(cout,"  A");
		Vector<float> b;
		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		Solve(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 Rank, Orthogonal/Null space" << endl;
		cout << "-----------------------------------------------------\n";
		A.Resize(3,3);
		A(1,1)=1; A(1,2)=2; A(1,3)=3;
		A(2,1)=1; A(2,2)=2; A(2,3)=3;
		A(3,1)=1; A(3,2)=2; A(3,3)=3;
		A.Print(cout,"  A");
		cout << "  Rank(A)=" << Rank(A) << endl;
		Orth(A).Print(cout,"  Orth(A)");
		Null(A).Print(cout,"  Null(A)");
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing condition numbers" << endl;
		cout << "-----------------------------------------------------\n";
#if !defined(USE_NRC_CODE)
		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 << "  Cond(A)=" << Cond(A) << endl;
		cout << "  CondLU(A)=" << CondLU(A) << endl;
#endif
		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");
		cout << "  Cond(A)=" << Cond(A) << endl;
	}
#endif

#if 1
	{
		/* A is the matrix whose eigenvalues and eigenvectors are sought */
		Matrix<float>   A(3,3,"0.846186767 0.468314366 0.236917463 "			"0.65696803 0.910127033 0.733477793 "			"0.593281737 0.283758285 0.218800796");

		cout << "-----------------------------------------------------\n";
		cout << "Testing the Hessenberg decomposition: A = P*H*P\'\n";
		cout << "-----------------------------------------------------\n";
		Matrix<float> H,P;
		Hess(A,H,P);
		A.Print(cout,"  A");
		H.Print(cout,"  H");
		P.Print(cout,"  P");
		(A-P*H*Transpose(P)).Print(cout,"  A-P*H*P\'=0");

		cout << "-----------------------------------------------------\n";
		cout << "Testing the Schur decomposition: A = Q*T*Q\'\n";
		cout << "-----------------------------------------------------\n";
		Matrix<float> T,Q;
		Schur(A,T,Q);
		A.Print(cout,"  A");
		T.Print(cout,"  T");
		Q.Print(cout,"  Q");
		(A-Q*T*Transpose(Q)).Print(cout,"  A-Q*T*Q\'=0");
		
		cout << "-----------------------------------------------------\n";
		cout << "Testing eigen analysis for general matrices\n";
		cout << "-----------------------------------------------------\n";
		Matrix<float> EVec,EVal;
		isOK=Eig(A,EVec,EVal);
		A.Print(cout,"  A");
		EVec.Print(cout,"  EVec");
		EVal.Print(cout,"  EVal");
		checkEigenVectors(A,EVec,EVal);

		A.Resize(2,2);
		A(1,1)=7; A(1,2)=10;
		A(2,1)=15; A(2,2)=22;
		isOK=Eig(A,EVec,EVal);
		A.Print(cout,"  A");
		EVec.Print(cout,"  EVec");
		EVal.Print(cout,"  EVal");
		checkEigenVectors(A,EVec,EVal);
	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		cout << "-----------------------------------------------------\n";
		cout << "Testing matrix multiplication, left and right division" << 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;

		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");
		mtl_ldivide(A,B).Print(cout,"  mtl_ldivide(A,B)=A.\\B");
		mtl_rdivide(A,B).Print(cout,"  mtl_rdivide(A,B)=A./B");
		mtl_mldivide(A,B).Print(cout,"  mtl_mldivide(A,B)=A\\B");
		mtl_mrdivide(A,B).Print(cout,"  mtl_mrdivide(A,B)=A/B");

	}
#endif

#if 1
	{
		Matrix<float> A;
		cout << "-----------------------------------------------------\n";
		cout << "Testing matrix square root, exponential, logarithm and power" << endl;
		cout << "-----------------------------------------------------\n";

		A.Resize(2,2);
		A(1,1)=1; A(1,2)=2;
		A(2,1)=3; A(2,2)=4;
		A=A*A;
		A.Print(cout,"  A");
		Sqrtm(A).Print(cout,"  Sqrtm(A)");

		A.Resize(2,2);
		A(1,1)=7; A(1,2)=10;
		A(2,1)=15; A(2,2)=22;
		A.Print(cout,"  A");
		Sqrtm(A).Print(cout,"  Sqrtm(A)");

		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");
		Expm(A).Print(cout,"  Expm(A)");

		A=Expm(A);
		A.Print(cout,"  A");
		Logm(A).Print(cout,"  Logm(A)");

		A.Resize(2,2);
		A(1,1)=1; A(1,2)=2; 
		A(2,1)=-3; A(2,2)=-4; 
		A.Print(cout,"  A");
		Pow(A,2.0).Print(cout,"  Pow(A,2.0)");
		Pow(A,-2.0).Print(cout,"  Pow(A,-2.0)");
		Pow(A,2.5).Print(cout,"  Pow(A,2.5)");
		Pow(A,-2.5).Print(cout,"  Pow(A,-2.5)");
		Powm(A,2.0).Print(cout,"  Powm(A,2.0)");
		Powm(A,-2.0).Print(cout,"  Powm(A,-2.0)");
		Powm(A,2.5).Print(cout,"  Powm(A,2.5)");
		Powm(A,-2.5).Print(cout,"  Powm(A,-2.5)");

	}
#endif

#if 1
	{
		Matrix<float> A;
		Matrix<float> B;
		Vector<float> x;
		int i;
		cout << "-----------------------------------------------------\n";
		cout << "Testing FFT" << endl;
		cout << "-----------------------------------------------------\n";
		int N=16;
		Vector<complex<float> > f(N);
		f(1)=complex<float>(1,0); f(2)=complex<float>(2,0); f(3)=complex<float>(3,0); f(4)=complex<float>(4,0);
		f(5)=complex<float>(5,0); f(6)=complex<float>(6,0); f(7)=complex<float>(7,0); f(8)=complex<float>(8,0);
		f.Print(cout,"  X");
		Vector<float> C(2*N);
		for(i=1;i<=N;i++){
			C[2*i-1]=f[i].real();
			C[2*i]=f[i].imag();
		}
		Vector<float> g=FFT(C,N);
		Vector<complex<float> > R(f.Size());
		for(i=1;i<=N;i++){
			R(i)=complex<float> (g[2*i-1],g[2*i]);
		}
		R.Print(cout,"  FFT(X,N)");
		
		for(i=1;i<=N;i++){
			C[2*i-1]=R[i].real();
			C[2*i]=R[i].imag();
		}
		Vector<float> h=IFFT(C,N);
		for(i=1;i<=N;i++){
			R(i)=complex<float> (h[2*i-1],h[2*i]);
		}
		R.Print(cout,"  IFFT(FFT(X))");
	}
#endif

#if 1
	{
		cout << "-----------------------------------------------------\n";
		cout << "Testing random number generators\n";
		cout << "   Outputs are saved in ./temp directory.\n";
		cout << "-----------------------------------------------------\n";

		int i;
		int n=10000;
		Vector<float> A(n);
		string resultDir="./temp/";
		string distName;
		if(CreateDir(resultDir.c_str())==0) {
			cerr << "ERROR: Failed in creating directory ./temp" << endl;
			exit(1);
		}

		distName="uniform";
		cout << "  Generating " << distName << " distributed random numbers\n";
		for(i=1;i<=n;i++){
			A[i]=(float)UniformRandom();
		}
		A.SaveMatlab(resultDir+distName+".dat");

		distName="gaussian";
		cout << "  Generating " << distName << " distributed random numbers\n";
		for(i=1;i<=n;i++){
			A[i]=(float)GaussianRandom();
		}
		A.SaveMatlab(resultDir+distName+".dat");

		distName="gamma";
		cout << "  Generating " << distName << " distributed random numbers\n";
		for(i=1;i<=n;i++){
			A[i]=(float)GammaRandom(5,1);
		}
		A.SaveMatlab(resultDir+distName+".dat");
	}
#endif

#if 1
	{
		cout << "-----------------------------------------------------\n";
		cout << "Testing Cholesky decomposition for complex matrices: A=L*L\'" << endl;
		cout << "-----------------------------------------------------\n";

		Matrix<complex<float> > A, l;
		Matrix<complex<float> > L;
		Vector<complex<float> > b,x;

		l.Resize(3,3);
		l(1,1)=complex<float> (1,1); l(1,2)=complex<float> (0,0); l(1,3)=complex<float> (0,0); 
		l(2,1)=complex<float> (2,-1); l(2,2)=complex<float> (3,-1); l(2,3)=complex<float> (0,0); 
		l(3,1)=complex<float> (4,1); l(3,2)=complex<float> (5,1); l(3,3)=complex<float> (6,1); 

		A=l*Transpose(l);

		// A=L*U;
		cout << "  A=L*U" << endl;
		A.Print(cout,"  A");
		Chol(A,L);
		L.Print(cout,"  L");
		(A-L*Transpose(L)).Print(cout,"  A-L*L\'=0");

		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		x.Resize(3);
		CholSolve(A,b,x);
		b.Print(cout,"  b");
		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		cout << "-----------------------------------------------------\n";
		cout << "Testing LU decomposition for complex matrices: A=L*U" << endl;
		cout << "-----------------------------------------------------\n";

		Matrix<complex<float> > A;
		Matrix<complex<float> > L,U,P;
		Vector<complex<float> > b,x;

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

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

		b.Resize(3);
		b(1)=1; b(2)=-2; b(3)=5;
		x.Resize(3);
		LuSolve(A,b,x);
		b.Print(cout,"  b");