CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: testInversion.cc,v 1.2 2003/08/13 20:00:12 garren Exp $ 00003 // 00004 // This file is a part of CLHEP - a Class Library for High Energy Physics. 00005 // 00006 // This is a collection of parts of programs that are useful 00007 // for testing matrix inversion algorithms 00008 // 9/97, Mario Stanke 00009 00010 #include <time.h> 00011 #include <iostream> 00012 00013 #include "CLHEP/Matrix/defs.h" 00014 #include "CLHEP/Matrix/Matrix.h" 00015 #include "CLHEP/Matrix/SymMatrix.h" 00016 #include "CLHEP/Matrix/DiagMatrix.h" 00017 00018 using std::cout; 00019 using std::endl; 00020 00021 using namespace CLHEP; 00022 00023 int main() 00024 { 00025 //int n , i, j, k, ierr1, ierr2; 00026 int n, j, ierr1, ierr2; 00027 time_t zeit1, zeit2; 00028 00029 // ****compare the speed of inversion algorithms 00030 00031 HepSymMatrix A(5,1); 00032 00033 // for (j=1;j <= 100; j++) 00034 // for (k=1; k<=j; k++) 00035 // A(j,k)=rand()%9-5; 00036 00037 A(1,1)=6; 00038 A(1,2)=.8545; 00039 A(1,3)=-.684651; 00040 A(1,4)=.372547; 00041 A(1,5)=.68151; 00042 //A(1,6)=.068151; 00043 A(2,2)=4; 00044 A(2,3)=.5151; 00045 A(2,4)=.5151; 00046 A(2,5)=.651651; 00047 //A(2,6)=.9651651; 00048 A(3,3)=5; 00049 A(3,4)=.3; 00050 A(3,5)=.363; 00051 //A(3,6)=.7363; 00052 A(4,4)=-3; 00053 A(4,5)=-.23753; 00054 //A(4,6)=-.23753; 00055 A(5,5)=-5; 00056 //A(5,6)=-2; 00057 //A(6,6)=-3; 00058 00059 HepSymMatrix B(A); 00060 HepSymMatrix D(A); 00061 HepSymMatrix C(5,0); 00062 HepMatrix M(A); 00063 00064 cout << "M inverse" << M.inverse(ierr2) << endl; 00065 00066 C = B.inverse(ierr1); 00067 D.invert(ierr2); 00068 00069 cout << "B " << B << endl; 00070 cout << "B inverse" << C << endl; 00071 #ifndef INSTALLATION_CHECK 00072 cout << "B * inverse" << B * C << endl; 00073 #endif 00074 cout << "ierr1: " << ierr1 << endl; 00075 00076 cout << "D * inverse" << D * C << endl; 00077 cout << "ierr2: " << ierr2 << endl; 00078 cout << "M inverse" << M.inverse(ierr2) << endl; 00079 #ifndef INSTALLATION_CHECK 00080 cout << "M * inverse" << M * M.inverse(ierr2)<< endl; 00081 #endif 00082 cout << "ierr2: " << ierr2 << endl; 00083 00084 #ifndef INSTALLATION_CHECK 00085 n = 100000; 00086 #else 00087 n = 10; 00088 #endif 00089 zeit1 = time(NULL); 00090 cout << "Executing " << n << " inversions ..." << endl; 00091 for (j=0; j<n; j++) 00092 { 00093 B.invert(ierr1); 00094 if (ierr1) 00095 cout << "B: error in invert" << endl; 00096 } 00097 zeit2 = time(NULL); 00098 cout << "B: duration of inversion: " << zeit2-zeit1 << " secs" << endl; 00099 00100 zeit1 = time(NULL); 00101 cout << "Executing " << n << " inversions ..." << endl; 00102 for (j=0; j<n; j++) 00103 { 00104 D.invert(ierr2); 00105 if (ierr2) 00106 cout << "D: error in invert" << endl; 00107 } 00108 zeit2 = time(NULL); 00109 cout << D << endl; 00110 cout << "D: duration of inversion: " << zeit2-zeit1 << " secs" << endl; 00111 00112 00113 /***** check correctness and compare results of inversion algorithms 00114 00115 double dist1, dist2; 00116 HepSymMatrix A(5,1), B(5), C(5), I(5,1); 00117 HepSymMatrix M(5); 00118 n = 200000; 00119 00120 for (j=1; j <= n ; j++) 00121 { 00122 A(1,1)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00123 A(1,2)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00124 A(1,3)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00125 A(1,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00126 A(1,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00127 //A(1,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00128 A(2,2)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00129 A(2,3)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00130 A(2,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00131 A(2,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00132 //A(2,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00133 A(3,3)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00134 A(3,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00135 A(3,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00136 //A(3,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00137 A(4,4)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00138 A(4,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00139 //A(4,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00140 A(5,5)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00141 //A(5,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00142 //A(6,6)=(((rand()%7-3)))+ (rand()%3-1)*1e-10; 00143 00144 M=B=C=A; 00145 00146 B.invert(ierr2); 00147 M.old_invert(ierr1); 00148 00149 dist2 = norm_infinity(B*C-I); 00150 dist1 = norm_infinity(M*C-I); 00151 00152 00153 if (ierr1 != ierr2) 00154 { 00155 cout << C << endl; 00156 cout << "B " << B << endl; 00157 cout << "B*C" << B*C << endl; 00158 cout << "M*C" << M*C << endl; 00159 cout << "M " << M << endl; 00160 cout << "determinant " << C.determinant() << endl; 00161 cout << "dist2 " << dist2 << endl; 00162 cout << "dist1 " << dist1 << endl; 00163 00164 cout << "ierr2 " << ierr2 <<endl; 00165 cout << "ierr1 " << ierr1 <<endl; 00166 00167 cout << "j " << j << endl; 00168 } 00169 00170 if (ierr2==0 && dist2 > 1e-4) 00171 { 00172 cout << "bunch failed to invert but did not indicate" << endl; 00173 cout << C << endl; 00174 cout << "B " << B << endl; 00175 cout << "B*C" << B*C << endl; 00176 cout << "M*C" << M*C << endl; 00177 cout << "determinant " << C.determinant() << endl; 00178 cout << "dist2 " << dist2 << endl; 00179 cout << "dist1 " << dist1 << endl; 00180 00181 cout << "ierr2 " << ierr2 <<endl; 00182 cout << "ierr1 " << ierr1 <<endl; 00183 00184 cout << "j " << j << endl; 00185 } 00186 if (ierr2==1) 00187 { 00188 // bunch thinks it is singular 00189 if (norm_infinity(M*C-I) < 1e-6) 00190 { 00191 cout << "bunch said it is singular but old could invert" 00192 << endl; 00193 cout << C << endl; 00194 cout << "B*C" << B*C << endl; 00195 cout << "M*C" << M*C << endl; 00196 cout << "determinant " << C.determinant() << endl; 00197 00198 cout << "dist2" << dist2 << endl; 00199 cout << "ierr2 " << ierr2 <<endl; 00200 cout << "ierr1 " << ierr1 <<endl; 00201 00202 cout << "j " << j << endl; 00203 } 00204 } 00205 00206 if (ierr1==0 && dist1 > 1e-4 && ierr2==0) 00207 { 00208 cout << "old failed to invert but did not indicate" << endl; 00209 00210 cout << C << endl; 00211 cout << "B*C" << B*C << endl; 00212 cout << "M*C" << M*C << endl; 00213 cout << "determinant " << C.determinant() << endl; 00214 00215 cout << "ierr2 " << ierr2 <<endl; 00216 cout << "ierr1 " << ierr1 <<endl; 00217 00218 cout << "dist1 " << dist1 << endl; 00219 cout << "j " << j << endl; 00220 return 0; 00221 00222 } 00223 if (ierr1==1) 00224 { 00225 // old thinks it is singular 00226 if (norm_infinity(B*C-I) < 1e-6) 00227 { 00228 cout << "old said it is singular but bunch could invert" 00229 << endl; 00230 00231 cout << C << endl; 00232 cout << "B*C" << B*C << endl; 00233 cout << "M*C" << M*C << endl; 00234 cout << "determinant " << C.determinant() << endl; 00235 00236 cout << "dist1" << dist1 << endl; 00237 cout << "dist2" << dist2 << endl; 00238 cout << "ierr2 " << ierr2 <<endl; 00239 cout << "ierr1 " << ierr1 <<endl; 00240 00241 cout << "j " << j << endl; 00242 return 0; 00243 00244 } 00245 } 00246 00247 } 00248 */ 00249 00250 /*** one tough symmetric matrix from real physical data 00251 00252 sm(1,1)=5347.51; 00253 sm(1,2)=-142756; 00254 sm(1,3)= -1.86624e+06; 00255 sm(1,4)=0.0164743; 00256 sm(1,5)=0.0915348; 00257 sm(1,6)=0.421851; 00258 sm(2,2)=3.81277; 00259 sm(2,3)=4.98668e+07; 00260 sm(2,4)=-0.0697533; 00261 sm(2,5)=12.8555; 00262 sm(2,6)=-2.01124; 00263 sm(3,3)=6.52498e+08; 00264 sm(3,4)=3.87491; 00265 sm(3,5)=365.965; 00266 sm(3,6)=93.3686; 00267 sm(4,4)=7.77672e-05; 00268 sm(4,5)=0.0032134; 00269 sm(4,6)=0.00194407; 00270 sm(5,5)=0.132845; 00271 sm(5,6)=0.0803294; 00272 sm(6,6)=0.0485992; 00273 00274 */ 00275 00276 /**** a group of near singular (nonsingular) matrices 00277 int n=5; 00278 HepSymMatrix sm(n); // nxn Hilbert Matrix 00279 for(i=1;i<=n;i++) 00280 for(k=i;k<=n;k++) 00281 sm(i,k)=1./(i+k-1); 00282 */ 00283 return 0; 00284 } // end of main 00285 00286 00287