Transcript ppt - Cern
Test and Validation Studies of Mathematical Software Libraries A summary of my work as a technical student at CERN LCG AA Meeting, 22. September 2004 M. Hatlo / CERN Special Functions • • • • • Comparison of numerical results Performance GSL-NAG-C and GSL-TMath Bessel functions Gamma, Logarithm of Gamma, Error function and Complementary Error function M. Hatlo / CERN Results • GSL performed very well compared to NAG-C – Difference usually less than the estimated error • Bigger differences between GSL and TMath • Little difference in time for GSL • NAG-C and TMath are faster M. Hatlo / CERN M. Hatlo / CERN M. Hatlo / CERN Timing of Special Functions • 100 000 function calls for Bessel I0, I1, J0 and J1 • 50 000 function calls for the rest M. Hatlo / CERN Distributions • Comparison of numerical results • Performance in the evaluation • Generation of random numbers according to distribution – Comparison and Kolmogorov-Smirnov test • Normal distribution, Landau distribution, Gamma distribution, Poisson distribution and Chi-Square distribution M. Hatlo / CERN M. Hatlo / CERN Timing Results for some Distributions • 1 000 000 function calls M. Hatlo / CERN Random Numbers According to Distribution M. Hatlo / CERN Random Numbers • Two tests – Frequency test – Point test • Main generators from GSL – – – – – – – – – gsl_rng_mt19937 gsl_rng_cmrg gsl_rng_mrg gsl_rng_taus gsl_rng_taus2 gsl_rng_gfsr4 gsl_rng_ranlux389 gsl_rng_ranlux gsl_rng_ranlxd2 M. Hatlo / CERN Frequency Test Fill space in d dimensions with points formed from a sequence of random numbers. Look in a small volume and the frequency as the number of bins which maximize |Nodd-Neven|. gsl_rng_minstd M. Hatlo / CERN Frequency Test • With this frequency, look other places in the space and compute Nodd. • Nodd should be normal distributed M. Hatlo / CERN Results • • • • • 10 results for Nodd Kolmogorov-Smirnov test Taus, 8 dim and Ranlux389, 6 dim New test for poor results All passed M. Hatlo / CERN Point Test • Arrange a sequence of random numbers into (1) (d ) multidimensional points Pi ( pi ,..., pi ) • Define distance between two points as d Pi Pj dist k , k 1 dist k min pi( k ) p (jk ) ,1 pi( k ) p (jk ) 1 k d • Find all points Pi that are closer to P1 than the mean-n*standard deviation (n=3,4,5) • Calculate the distance between Pi+1 and P2 M. Hatlo / CERN Point test • For d 3 the distance should be normal distributed. • Use Kolmogorov-Smirnov test • All generators pass M. Hatlo / CERN Numerical Integration • Wrapper for existing gsl algorithms • Tested on a few number of integrals – Compare numerical results with analytical results • No difference larger than 10-7 (input tolerance), but need further testing M. Hatlo / CERN Integrals 1 1 1 0 log x 1 x dx C 0.5772157... log x (*) dx 4.0 x 0 1 (*) 2 x t e 2 /2 dt erf ( x 1 x ) 2 dx 0 1 x x 2 27 1 x p 1 (*) dx 1 x sin( p ) 0 dx 0 (1 x) x dx 0 a 2 sin 2 x b 2 cos 2 x ab 1 p x log 0 1 1 dx x ( p 1) 2 xe x dx 0 x e dx 2 0 2 2 2 0 e x1 dx 6 x x3 4 0 e x1 dx 15 sin x dx 0 x 2 (*) Lorenzianp eak Background M. Hatlo / CERN Performance in Numerical Integration • Quadrature routines • QAG – adaptive integration • QAGUI – adaptive integration from zero to infinity • QAGS – adaptive integration with singularities • QNG – non-adaptive Gauss-Kronrod integration • NB! Different integrals are used, marked with (*) on last slide M. Hatlo / CERN Linear Algebra • E. Myklebust summer student 2003 – A Comparative study of Numerical Linear Algebra Libraries • Particle Track Reconstruction, Kalman filter update equations – Multiplication, addition, inversion and transpose • Originally 2x2, 2x5, 5x2 and 5x5 • Extended to bigger matrices – 2x5, 4x10, 10x25 and 20x50 • CLHEP, uBLAS, LAPACK, GSL and ROOT • Used timer from SEAL base M. Hatlo / CERN Results (Linux, P4 1.8 MHz ) • High RMS for GSL and LAPACK in 2x5 • Error with ROOT for 20x50 M. Hatlo / CERN Conclusions • GSL performs reasonably good • Both tests of randomness were passed by all the main generators from GSL • More testing is needed for the numerical integration • All test programs are in the SEAL cvs repository • A test suite can be easily created and automatically run for new SEAL releases • A written report of my work will be put on the webpage when finished M. Hatlo / CERN Results UBLAS 1 UBLAS 2 UBLAS 3 LAPACK CLHEP GSL ROOT 2x5 0.0000291446 (1.59) 0.000036635 (2.00) 0.0000345755 (1.89) 0.0000262477 (1.44) 0.0000254471 1.39) 0.0000315174 (1.72) 0.000018277 (1.00) 4x10 0.000122262 (6.69) 0.000128379 (7.02) 0.00004496 (2.46) 0.0000558439 (3.06) 0.0000591721 (3.24) 0.0000387867 (2.12) 10x25 0.00162649 (89.0) 0.000313892 (17.2) 0.000559981 (30.6) 0.000401781 (22.0) 0.000277101 (15.2) 20x50 0.0148939 (815) 0.00125005 (68.4) 0.00302708 (166) 0.0014349 (78.5) M. Hatlo / CERN