CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // testBug6181.cc 00002 // 00003 // The bug reported an erroneous result in inverting a 7x7 matrix. 00004 // 00005 // 7x7 represents the break point between the low-n methods and 00006 // the general-size method. This test program tests the case triggering 00007 // the bug report, and also inversion for matrices (both symmetric and 00008 // general) of sizes 1 - 9. 00009 00010 #include <iostream> 00011 #include <math.h> 00012 00013 #include "CLHEP/Matrix/Matrix.h" 00014 #include "CLHEP/Matrix/SymMatrix.h" 00015 00016 int test_inversion (int N) { 00017 00018 int i; 00019 int j; 00020 CLHEP::HepMatrix M(N,N,0); 00021 CLHEP::HepSymMatrix S(N,0); 00022 for(i=1;i<=N;++i) { 00023 for(j=1;j<=N;++j) { 00024 if(i<=j) { 00025 S (i,j) = 10*i+j; 00026 M (i,j) = 10*i+j; 00027 M (j,i) = M(i,j) + .1 * (i-j)*(i-j); 00028 } 00029 } 00030 } 00031 CLHEP::HepMatrix MM(N,N,0); 00032 CLHEP::HepSymMatrix SS(N,0); 00033 MM = M; 00034 SS = S; 00035 int ierr = 0; 00036 MM.invert(ierr); 00037 if (ierr) { 00038 std::cout<<"MM.invert failed!!!! N = " << N << 00039 " ierr = "<< ierr <<std::endl; 00040 return 1+100*N; 00041 } 00042 ierr = 0; 00043 SS.invert(ierr); 00044 if (ierr) { 00045 std::cout<<"SS.invert failed!!!! N = " << N << 00046 " ierr = "<< ierr <<std::endl; 00047 return 2+100*N; 00048 } 00049 CLHEP::HepMatrix MI(N,N,0); 00050 CLHEP::HepMatrix SI(N,N,0); 00051 CLHEP::HepMatrix MS(N,N,0); 00052 CLHEP::HepMatrix MSS(N,N,0); 00053 MI = MM*M; 00054 MS = S; 00055 MSS = SS; 00056 SI = MSS*MS; 00057 for(i=1;i<=N;++i) { 00058 for(j=1;j<=N;++j) { 00059 if(i!=j) { 00060 if (fabs(MI(i,j)) > 1.0e-12) { 00061 std::cout<<"MM.invert incorrect N = " << N << 00062 " error = "<< fabs(MI(i,j)) <<std::endl; 00063 return 3+100*N; 00064 } 00065 if (fabs(SI(i,j)) > 1.0e-12) { 00066 std::cout<<"SS.invert incorrect N = " << N << 00067 " error = "<< fabs(SI(i,j)) <<std::endl; 00068 return 4+100*N; 00069 } 00070 } 00071 if(i==j) { 00072 if (fabs(1-MI(i,j)) > 1.0e-12) { 00073 std::cout<<"MM.invert incorrect N = " << N << 00074 " error = "<< fabs(1-MI(i,j)) <<std::endl; 00075 return 3+100*N; 00076 } 00077 if (fabs(1-SI(i,j)) > 1.0e-12) { 00078 std::cout<<"SS.invert incorrect N = " << N << 00079 " error = "<< fabs(1-SI(i,j)) <<std::endl; 00080 return 4+100*N; 00081 } 00082 } 00083 } 00084 } 00085 // Finally, the case that actually (for N=7) triggered bug report 6181: 00086 M = S; 00087 MM = M; 00088 MM.invert(ierr); 00089 if (ierr) { 00090 std::cout<<"MM.invert for symmetric case failed!!!! N = " << N << 00091 " ierr = "<< ierr <<std::endl; 00092 return 5+100*N; 00093 } 00094 MI = MM*M; 00095 for(i=1;i<=N;++i) { 00096 for(j=1;j<=N;++j) { 00097 if(i!=j) { 00098 if (fabs(MI(i,j)) > 1.0e-12) { 00099 std::cout<<"MM.invert incorrect N = " << N << 00100 " error = "<< fabs(MI(i,j)) <<std::endl; 00101 return 6+100*N; 00102 } 00103 } 00104 if(i==j) { 00105 if (fabs(1-MI(i,j)) > 1.0e-12) { 00106 std::cout<<"MM.invert incorrect N = " << N << 00107 " error = "<< fabs(1-MI(i,j)) <<std::endl; 00108 return 7+100*N; 00109 } 00110 } 00111 } 00112 } 00113 return 0; 00114 } 00115 00116 00117 int main(int, char **) { 00118 int ret=0; 00119 for (int i=1; i<10; i++) { 00120 ret = test_inversion(i); 00121 if (ret) return ret; 00122 } 00123 return ret; 00124 } 00125