?? 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;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -