ezmtltest.cpp

矩阵运算的模板类

		x.Print(cout,"  x");
	}
#endif

#if 1
	{
		Matrix<complex<float> > A;
		Vector<complex<float> > b,x;
		Matrix<complex<float> > Q,R;

		cout << "-----------------------------------------------------\n";
		cout << "Testing QR decomposition for complex matrices: A=Q*R" << endl;
		cout << "-----------------------------------------------------\n";

		A.Resize(3,3);
		A(1,1)=complex<float> (1,1); A(1,2)=complex<float> (2,0); A(1,3)=complex<float> (3,1); 
		A(2,1)=complex<float> (3,0); A(2,2)=complex<float> (5,1); A(2,3)=complex<float> (1,0); 
		A(3,1)=complex<float> (1,0); A(3,2)=complex<float> (-1,0); A(3,3)=complex<float> (0,0); 
		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<complex<float> > A;
		cout << "-----------------------------------------------------\n";
		cout << "Testing condition numbers for complex matrices" << endl;
		cout << "-----------------------------------------------------\n";
#if !defined(USE_NRC_CODE)
		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.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
	{
		cout << "-----------------------------------------------------\n";
		cout << "Testing eigen analysis for Hermitian matrices\n";
		cout << "-----------------------------------------------------\n";

		Matrix<complex<float> > A, l;
		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);

		Matrix<complex<float> > EVec,EVal;
		isOK=Eig(A,EVec,EVal);
		if(isOK==0){
			cerr << "ERROR: EigComplex()\n";
		}
		A.Print(cout,"  A");
		EVec.Print(cout,"  EVec");
		EVal.Print(cout,"  EVal");
		int i;
		for(i=1;i<=3;i++){
			Vector<complex<float> > x=EVec.Col(i);
			Vector<complex<float> > y=Transpose(x);
			complex<float> z = Multiply3(Transpose(x),A,x);
			cout << "  EVec(" << i << ")\'*A*Evec(" << i << ")=" << z << endl;
		}
		Matrix<complex<float> > ans=A*EVec-EVec*EVal;
		ans.Print(cout,"  A*EVec-EVec*EVal=0");
		cout << endl;
	}
#endif

#if 1
	{
		int i;
		Matrix<complex<float> > A;

		A.Resize(3,3);
#if 1
		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);
#else
		A(1,1)=complex<float> (1,0); A(1,2)=complex<float> (2,0); A(1,3)=complex<float> (3,0); 
		A(2,1)=complex<float> (4,0); A(2,2)=complex<float> (5,0); A(2,3)=complex<float> (6,0); 
		A(3,1)=complex<float> (7,0); A(3,2)=complex<float> (8,0); A(3,3)=complex<float> (0,0);
#endif

		cout << "-----------------------------------------------------\n";
		cout << "Testing the Hessenberg decomposition for complex matrices: A = P*H*P\'\n";
		cout << "-----------------------------------------------------\n";
		Matrix<complex<float> > H,P;
		HessComplex(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 for complex matrices: A = Q*T*Q\'\n";
		cout << "-----------------------------------------------------\n";
		Matrix<complex<float> > T,Q;
		isOK=Schur(A,T,Q);
		if(isOK==0){
			cerr << "ERROR: SchurComplex()\n";
		}
		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 complex matrices\n";
		cout << "-----------------------------------------------------\n";
		Matrix<complex<float> > EVec,EVal;
		isOK=Eig(A,EVec,EVal);
		if(isOK==0){
			cerr << "ERROR: EigComplex()\n";
		}
		A.Print(cout,"  A");
		EVec.Print(cout,"  EVec");
		EVal.Print(cout,"  EVal");	
		for(i=1;i<=3;i++){
			Vector<complex<float> > x=EVec.Col(i);
			Vector<complex<float> > y=Transpose(x);
			complex<float> z = Multiply3(Transpose(x),A,x);
			cout << "  EVec(" << i << ")\'*A*Evec(" << i << ")=" << z << endl;
		}
		Matrix<complex<float> > ans=A*EVec-EVec*EVal;
		ans.Print(cout,"  A*EVec-EVec*EVal=0");
		cout << endl;

		cout << "  -----------------------------------------------------\n";
		cout << "  Testing with random matrices\n";
		cout << "  -----------------------------------------------------\n";
		for(int k=1;k<=10;k++){
			A=Matrix<complex<float> >::CRandn(10,10);
			isOK=Eig(A,EVec,EVal);
			if(isOK==0){
				cerr << "ERROR: EigComplex()\n";
			}
			for(i=1;i<=3;i++){
				Vector<complex<float> > x=EVec.Col(i);
				Vector<complex<float> > y=Transpose(x);
				complex<float> z = Multiply3(Transpose(x),A,x);
				if(Abs(z-EVal(i,i))>1e-5){
					cerr<<"ERROR2: Eig()\n";
				}
			}
			ans=A*EVec-EVec*EVal;
			if(ans.IsZero(1e-5)==0){
				cerr<<"ERROR1: Eig()\n";
			}
			cout << "  Matrix " << k << " OK.\n";
		}
	}
#endif

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

		A.Resize(2,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); 

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

		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<complex<float> >::CRandn(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";
			}	
		}
	}
#endif

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

#endif

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

		A.Resize(2,2);
		A(1,1)=complex<float> (1,1); A(1,2)=complex<float> (2,2); 
		A(2,1)=complex<float> (3,-3); A(2,2)=complex<float> (4,-4); 

		A = A*A;
		A.Print(cout,"  A");
		Sqrtm(A).Print(cout,"  Sqrtm(A)");

		A.Resize(2,2);
		A(1,1)=complex<float> (1,1); A(1,2)=complex<float> (2,2); 
		A(2,1)=complex<float> (3,-3); A(2,2)=complex<float> (4,-4); 

		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(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)=complex<float> (1,1); A(1,2)=complex<float> (2,2); 
		A(2,1)=complex<float> (3,-3); A(2,2)=complex<float> (4,-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<complex<float> > A;
		cout << "-----------------------------------------------------\n";
		cout << "Testing Determinant" << endl;
		cout << "-----------------------------------------------------\n";
		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.Print(cout,"A");
		cout << "  det(A)=" << complex<float>(Det(A)) << endl << endl;

		A.Resize(2,2);
		A(1,1)=complex<float> (1,1); A(1,2)=complex<float> (2,2); 
		A(2,1)=complex<float> (3,-3); A(2,2)=complex<float> (4,-4); 
		A.Print(cout,"A");
		cout << "  det(A)=" << complex<float>(Det(A)) << endl << endl;
	}
#endif


	return 0;

}


template <typename T> int checkEigenVectors(Matrix<T>& A,Matrix<T>& EVec,Matrix<T>& EVal)
{
	if(EVal.IsDiagonal()){
		cout << "  >>> The A matrix has real eigenvalues.\n";
		Matrix<T> ans=A*EVec-EVec*EVal;
		ans.Print(cout,"  A*EVec-EVec*EVal=0");
	}
	else{
		cout << "  >>> The A matrix has complex eigenvalues.\n";
		Matrix<complex<T> > CVec, CVal;
		Matrix<T>::Compact2Complex(EVec,EVal,CVec,CVal);
		Matrix<complex<T> > C(NumRows(A),NumCols(A));
		for(int i=1;i<=NumRows(A);i++){
			for(int j=1;j<=NumCols(A);j++){
				C(i,j)=complex<T>(A(i,j),0);
			}
		}
		Matrix<complex<T> > ans=C*CVec-CVec*CVal;
		ans.Print(cout,"  C*CVec-CVec*CVal=0");

		Matrix<T> AVec,AVal;
		Matrix<T>::Complex2Compact(CVec,CVal,AVec,AVal);
		if((EVec-AVec).IsZero()==false){
			cout << "ERROR: Error in Complex2Compact().\n";
			CVec.Print(cout,"CVec");
			AVec.Print(cout,"AVec");
			assert(0);
		}
		if((EVal-AVal).IsZero()==false){
			cout << "ERROR: Error in Complex2Compact().\n";
			CVal.Print(cout,"CVal");
			AVal.Print(cout,"AVal");
			assert(0);
		}
	}
	return 1;
}


static int CreateDir(const char *dirName)
{
	/* Check for existence */
	if( (_access( dirName, 0 )) != -1 ) {
		/* Check for write permission */
		if( (_access( dirName, 2 )) != -1 ) return 1;
		else {
			fprintf(stderr, "Directory %s does not have write permission\n",dirName);
			return 0;
		}
	}
	else{
		fprintf(stderr, "Creating the directory %s\n",dirName);
#ifdef _WIN32
		if(_mkdir(dirName) == -1){
#else
		if(_mkdir(dirName,0755) == -1){
#endif
			fprintf(stderr, "Failed in creating Directory %s\n",dirName);
			return 0;
		}
	}
	return 1;
}