Skip to content
Snippets Groups Projects
Commit 8dc3fbcf authored by Lorenzo Moneta's avatar Lorenzo Moneta
Browse files

update tests for using also CLHEP 

git-svn-id: http://root.cern.ch/svn/root/trunk@14838 27541ba8-7e3a-0410-8455-c3a389f83636
parent 691c43ec
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
smatrix/test/Makefile
+ 9
− 2
View file @ 8dc3fbcf
...@@ -11,8 +11,15 @@ include ../../test/Makefile.arch ...@@ -11,8 +11,15 @@ include ../../test/Makefile.arch
#------------------------------------------------------------------------------ #------------------------------------------------------------------------------
ifeq ($(PLATFORM),macosx) ifeq ($(PLATFORM),macosx)
CXXFLAGS+= -g -funroll-loops #unroll loop better on gcc > 4
#CXXFLAGS+= -g CXXFLAGS+= -g -funroll-loops
endif
# if have clhep
#CLHEPBASE=/Users/moneta/mathlibs/CLHEP-1.9.2.2
ifneq ($(CLHEPBASE),)
CXXFLAGS+= -I$(CLHEPBASE)/include -DHAVE_CLHEP
EXTRALIBS+= $(CLHEPBASE)/lib/libCLHEP.a
endif endif
ifneq ($(PLATFORM),win32) ifneq ($(PLATFORM),win32)
......
smatrix/test/matrixOperations.C
+ 60
− 21
View file @ 8dc3fbcf
...@@ -4,11 +4,12 @@ bool drawSingleGraph = false; ...@@ -4,11 +4,12 @@ bool drawSingleGraph = false;
int topX=10; int topX=10;
int topY=50; int topY=50;
void matrixOperations_do(std::string type = "");
void matrixOperations(std::string type = "") { void matrixOperations_do(std::string type = "", bool clhep=false);
matrixOperations_do(type); void matrixOperations(std::string type = "",bool clhep=false) {
matrixOperations_do(type,clhep);
// matrixOperations_do("slc3_ia32_gcc323"); // matrixOperations_do("slc3_ia32_gcc323");
// matrixOperations_do("win32_vc71"); // matrixOperations_do("win32_vc71");
...@@ -96,10 +97,10 @@ TGraphErrors * hd = h1; ...@@ -96,10 +97,10 @@ TGraphErrors * hd = h1;
TLegend * tleg = new TLegend(0.78, 0.25, 0.97 ,0.45); TLegend * tleg = new TLegend(0.78, 0.25, 0.97 ,0.45);
tleg->AddEntry(h1, "SMatrix", "p"); tleg->AddEntry(h1, "SMatrix", "p");
tleg->AddEntry(h2, "TMatrix", "p"); tleg->AddEntry(h2, "TMatrix", "p");
if (h3) tleg->AddEntry(h3, "SMatrix_sym", "p"); if (h3 != 0) tleg->AddEntry(h3, "SMatrix_sym", "p");
if (h4) tleg->AddEntry(h4, "TMatrix_sym", "p"); if (h4 != 0) tleg->AddEntry(h4, "TMatrix_sym", "p");
if (h5) tleg->AddEntry(h5, "HepMatrix", "l"); if (h5 != 0) tleg->AddEntry(h5, "HepMatrix", "l");
if (h6) tleg->AddEntry(h6, "HepMatrix_sym", "l"); if (h6 != 0) tleg->AddEntry(h6, "HepMatrix_sym", "l");
tleg->Draw(); tleg->Draw();
...@@ -139,7 +140,7 @@ void GetData(std::string s,double * x, double * y, double * ey) { ...@@ -139,7 +140,7 @@ void GetData(std::string s,double * x, double * y, double * ey) {
void matrixOperations_do(std::string type ) { void matrixOperations_do(std::string type, bool clhep) {
systemName = type; systemName = type;
...@@ -163,46 +164,63 @@ void matrixOperations_do(std::string type ) { ...@@ -163,46 +164,63 @@ void matrixOperations_do(std::string type ) {
double cmat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 }; double cmat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 };
double ymat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 }; double ymat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 };
double wmat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 }; double wmat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 };
double zmat[N]; // = { 2., 3., 4., 5,8,10,300,800,5000 };
double es[N]; double es[N];
double et[N]; double et[N];
double ec[N]; double ec[N];
double ey[N]; double ey[N];
double ew[N]; double ew[N];
double ez[N];
c1->cd(1); c1->cd(1);
GetData("SMatrix_dot",x,smat,es); GetData("SMatrix_dot",x,smat,es);
GetData("TMatrix_dot",x,tmat,et); GetData("TMatrix_dot",x,tmat,et);
//GetData("HepMatrix_dot",x,cmat,ec); if (clhep) GetData("HepMatrix_dot",x,cmat,ec);
TGraphErrors * g10 = new TGraphErrors(nb,x, smat,0, es); TGraphErrors * g10 = new TGraphErrors(nb,x, smat,0, es);
TGraphErrors * g20 = new TGraphErrors(nb,x, tmat,0, et); TGraphErrors * g20 = new TGraphErrors(nb,x, tmat,0, et);
//TGraphErrors * g30 = new TGraphErrors(nb,x, cmat,0, ec); TGraphErrors * g30 = 0;
DrawData("#vec{v} #upoint #vec{w}",g10,g20); TGraphErrors * g40 = 0;
TGraphErrors * g50 = 0;
if (clhep) g50 = new TGraphErrors(nb,x, cmat,0, ec);
DrawData("#vec{v} #upoint #vec{w}",g10,g20,g30,g40,g50);
c1->cd(2); c1->cd(2);
GetData("SMatrix_M*V+",x,smat,es); GetData("SMatrix_M*V+",x,smat,es);
GetData("TMatrix_M*V+",x,tmat,et); GetData("TMatrix_M*V+",x,tmat,et);
GetData("SMatrix_sym_M*V+",x,ymat,ey); GetData("SMatrix_sym_M*V+",x,ymat,ey);
GetData("TMatrix_sym_M*V+",x,wmat,ew); GetData("TMatrix_sym_M*V+",x,wmat,ew);
//GetData("HepMatrix_M*V+",x,cmat,ec); if (clhep) GetData("HepMatrix_M*V+",x,cmat,ec);
if (clhep) GetData("HepMatrix_sym_M*V+",x,zmat,ez);
TGraphErrors * g11 = new TGraphErrors(nb,x, smat,0,es); TGraphErrors * g11 = new TGraphErrors(nb,x, smat,0,es);
TGraphErrors * g21 = new TGraphErrors(nb,x, tmat,0,et); TGraphErrors * g21 = new TGraphErrors(nb,x, tmat,0,et);
TGraphErrors * g31 = new TGraphErrors(nb,x, ymat,0,ey); TGraphErrors * g31 = new TGraphErrors(nb,x, ymat,0,ey);
TGraphErrors * g41 = new TGraphErrors(nb,x, wmat,0,ew); TGraphErrors * g41 = new TGraphErrors(nb,x, wmat,0,ew);
DrawData("M #upoint #vec{v} + #vec{w}",g11,g21,g31,g41); TGraphErrors * g51 = 0;
TGraphErrors * g61 = 0;
if (clhep) g51 = new TGraphErrors(nb,x, cmat,0, ec);
if (clhep) g61 = new TGraphErrors(nb,x, zmat,0, ez);
DrawData("M #upoint #vec{v} + #vec{w}",g11,g21,g31,g41,g51,g61);
c1->cd(3); c1->cd(3);
GetData("SMatrix_prod",x,smat,es); GetData("SMatrix_prod",x,smat,es);
GetData("TMatrix_prod",x,tmat,et); GetData("TMatrix_prod",x,tmat,et);
GetData("SMatrix_sym_prod",x,ymat,ey); GetData("SMatrix_sym_prod",x,ymat,ey);
GetData("TMatrix_sym_prod",x,wmat,ew); GetData("TMatrix_sym_prod",x,wmat,ew);
//GetData("HepMatrix_M*V+",x,cmat,ec); if (clhep) {
GetData("HepMatrix_M*V+",x,cmat,ec);
GetData("HepMatrix_sym_M*V+",x,zmat,ez);
}
TGraphErrors * g12 = new TGraphErrors(nb,x, smat,0,es); TGraphErrors * g12 = new TGraphErrors(nb,x, smat,0,es);
TGraphErrors * g22 = new TGraphErrors(nb,x, tmat,0,et); TGraphErrors * g22 = new TGraphErrors(nb,x, tmat,0,et);
TGraphErrors * g32 = new TGraphErrors(nb,x, ymat,0,ey); TGraphErrors * g32 = new TGraphErrors(nb,x, ymat,0,ey);
TGraphErrors * g42 = new TGraphErrors(nb,x, wmat,0,ew); TGraphErrors * g42 = new TGraphErrors(nb,x, wmat,0,ew);
DrawData("v^{T} * M * v",g11,g21,g31,g41); TGraphErrors * g52 = 0;
TGraphErrors * g62 = 0;
if (clhep) g52 = new TGraphErrors(nb,x, cmat,0, ec);
if (clhep) g62 = new TGraphErrors(nb,x, zmat,0, ez);
DrawData("v^{T} * M * v",g11,g21,g31,g41,g52,g62);
c1->cd(4); c1->cd(4);
...@@ -210,36 +228,57 @@ void matrixOperations_do(std::string type ) { ...@@ -210,36 +228,57 @@ void matrixOperations_do(std::string type ) {
GetData("TMatrix_M*M",x,tmat,et); GetData("TMatrix_M*M",x,tmat,et);
GetData("SMatrix_sym_M*M",x,ymat,ey); GetData("SMatrix_sym_M*M",x,ymat,ey);
GetData("TMatrix_sym_M*M",x,wmat,ew); GetData("TMatrix_sym_M*M",x,wmat,ew);
//GetData("HepMatrix_M*M",x,cmat,ec); if (clhep) {
GetData("HepMatrix_M*M",x,cmat,ec);
GetData("HepMatrix_sym_M*M",x,zmat,ez);
}
TGraphErrors * g14 = new TGraphErrors(nb,x, smat,0,es); TGraphErrors * g14 = new TGraphErrors(nb,x, smat,0,es);
TGraphErrors * g24 = new TGraphErrors(nb,x, tmat,0,et); TGraphErrors * g24 = new TGraphErrors(nb,x, tmat,0,et);
TGraphErrors * g34 = new TGraphErrors(nb,x, ymat,0,ey); TGraphErrors * g34 = new TGraphErrors(nb,x, ymat,0,ey);
TGraphErrors * g44 = new TGraphErrors(nb,x, wmat,0,ew); TGraphErrors * g44 = new TGraphErrors(nb,x, wmat,0,ew);
DrawData("A * B + C",g14,g24,g34,g44); TGraphErrors * g54 = 0;
TGraphErrors * g64 = 0;
if (clhep) g54 = new TGraphErrors(nb,x, cmat,0, ec);
if (clhep) g64 = new TGraphErrors(nb,x, zmat,0, ez);
DrawData("A * B + C",g14,g24,g34,g44,g54,g64);
c1->cd(5); c1->cd(5);
GetData("SMatrix_At*M*A",x,smat,es); GetData("SMatrix_At*M*A",x,smat,es);
GetData("TMatrix_At*M*A",x,tmat,et); GetData("TMatrix_At*M*A",x,tmat,et);
GetData("SMatrix_sym_At*M*A",x,ymat,ey); GetData("SMatrix_sym_At*M*A",x,ymat,ey);
GetData("TMatrix_sym_At*M*A",x,wmat,ew); GetData("TMatrix_sym_At*M*A",x,wmat,ew);
//GetData("HepMatrix_At*M*A",x,cmat,ec); if (clhep) {
GetData("HepMatrix_At*M*A",x,cmat,ec);
GetData("HepMatrix_sym_At*M*A",x,zmat,ez);
}
TGraphErrors * g15 = new TGraphErrors(nb,x, smat,0,es); TGraphErrors * g15 = new TGraphErrors(nb,x, smat,0,es);
TGraphErrors * g25 = new TGraphErrors(nb,x, tmat,0,et); TGraphErrors * g25 = new TGraphErrors(nb,x, tmat,0,et);
TGraphErrors * g35 = new TGraphErrors(nb,x, ymat,0,ey); TGraphErrors * g35 = new TGraphErrors(nb,x, ymat,0,ey);
TGraphErrors * g45 = new TGraphErrors(nb,x, wmat,0,ew); TGraphErrors * g45 = new TGraphErrors(nb,x, wmat,0,ew);
DrawData("A * M * A^{T}",g15,g25,g35,g45); TGraphErrors * g55 = 0;
TGraphErrors * g65 = 0;
if (clhep) g55 = new TGraphErrors(nb,x, cmat,0, ec);
if (clhep) g65 = new TGraphErrors(nb,x, zmat,0, ez);
DrawData("A * M * A^{T}",g15,g25,g35,g45,g55,g65);
c1->cd(6); c1->cd(6);
GetData("SMatrix_inv",x,smat,es); GetData("SMatrix_inv",x,smat,es);
GetData("TMatrix_inv",x,tmat,et); GetData("TMatrix_inv",x,tmat,et);
GetData("SMatrix_sym_inv",x,ymat,ey); GetData("SMatrix_sym_inv",x,ymat,ey);
GetData("TMatrix_sym_inv",x,wmat,ew); GetData("TMatrix_sym_inv",x,wmat,ew);
//GetData("HepMatrix_inv",x,cmat,ec); if (clhep) {
GetData("HepMatrix_inv",x,cmat,ec);
GetData("HepMatrix_sym_inv",x,zmat,ez);
}
TGraphErrors * g16 = new TGraphErrors(nb,x, smat,0,es); TGraphErrors * g16 = new TGraphErrors(nb,x, smat,0,es);
TGraphErrors * g26 = new TGraphErrors(nb,x, tmat,0,et); TGraphErrors * g26 = new TGraphErrors(nb,x, tmat,0,et);
TGraphErrors * g36 = new TGraphErrors(nb,x, ymat,0,ey); TGraphErrors * g36 = new TGraphErrors(nb,x, ymat,0,ey);
TGraphErrors * g46 = new TGraphErrors(nb,x, wmat,0,ew); TGraphErrors * g46 = new TGraphErrors(nb,x, wmat,0,ew);
DrawData("A^{-1}",g16,g26,g36,g46); TGraphErrors * g56 = 0;
TGraphErrors * g66 = 0;
if (clhep) g56 = new TGraphErrors(nb,x, cmat,0, ec);
if (clhep) g66 = new TGraphErrors(nb,x, zmat,0, ez);
DrawData("A^{-1}",g16,g26,g36,g46,g56,g66);
// TCanvas::Update() draws the frame, after which one can change it // TCanvas::Update() draws the frame, after which one can change it
......
smatrix/test/matrix_op.h
+ 86
− 8
View file @ 8dc3fbcf
...@@ -126,6 +126,7 @@ void testGMV(const M & mat, const V & v1, const V & v2, double & time, V & resul ...@@ -126,6 +126,7 @@ void testGMV(const M & mat, const V & v1, const V & v2, double & time, V & resul
} }
} }
// general matrix matrix op // general matrix matrix op
template<class A, class B, class C> template<class A, class B, class C>
void testMM(const A & a, const B & b, const C & c, double & time, C & result) { void testMM(const A & a, const B & b, const C & c, double & time, C & result) {
...@@ -208,7 +209,8 @@ void testATBA_S(const A & a, const B & b, double & time, C & result) { ...@@ -208,7 +209,8 @@ void testATBA_S(const A & a, const B & b, double & time, C & result) {
//result = Transpose(a) * b * a; //result = Transpose(a) * b * a;
//result = a * b * Transpose(a); //result = a * b * Transpose(a);
//result = a * b * a; //result = a * b * a;
btmp(0,0) = gV[l]; btmp(0,0) = gV[l];
// A at = Transpose(a);
C tmp = btmp * Transpose(a); C tmp = btmp * Transpose(a);
result = a * tmp; result = a * tmp;
} }
...@@ -218,6 +220,7 @@ void testATBA_S(const A & a, const B & b, double & time, C & result) { ...@@ -218,6 +220,7 @@ void testATBA_S(const A & a, const B & b, double & time, C & result) {
template<class A, class B, class C> template<class A, class B, class C>
void testATBA_S2(const A & a, const B & b, double & time, C & result) { void testATBA_S2(const A & a, const B & b, double & time, C & result) {
B btmp = b; B btmp = b;
test::Timer t(time,"At*M*A"); test::Timer t(time,"At*M*A");
for (int l = 0; l < NLOOP; l++) for (int l = 0; l < NLOOP; l++)
{ {
...@@ -432,56 +435,131 @@ void testMT_T(const A & a, double & time, C & result) { ...@@ -432,56 +435,131 @@ void testMT_T(const A & a, double & time, C & result) {
//smatrix //smatrix
template<class V> template<class V>
double testDot_C(const V & v1, const V & v2, double & time) { double testDot_C(const V & v1, const V & v2, double & time) {
V vtmp = v2;
test::Timer t(time,"dot "); test::Timer t(time,"dot ");
double result=0; double result=0;
for (int l = 0; l < 10*NLOOP; l++) for (int l = 0; l < 10*NLOOP; l++)
{ {
result += dot(v1,v2); vtmp[0] = gV[l];
result = dot(v1,v2);
} }
return result; return result;
} }
template<class M, class V> template<class M, class V>
double testInnerProd_C(const M & a, const V & v, double & time) { double testInnerProd_C(const M & a, const V & v, double & time) {
V vtmp = v;
test::Timer t(time,"prod"); test::Timer t(time,"prod");
double result=0; double result=0;
for (int l = 0; l < NLOOP; l++) for (int l = 0; l < NLOOP; l++)
{ {
V tmp = a*v; vtmp[0] = gV[l];
result += dot(v,tmp); V tmp = a*vtmp;
result = dot(vtmp,tmp);
} }
return result; return result;
} }
// matrix assignmnent(index starts from 1)
template<class M>
void testMeq_C(const M & m, double & time, M & result) {
M mtmp = m;
test::Timer t(time,"M=M ");
for (int l = 0; l < NLOOP; l++)
{
mtmp(1,1) = gV[l];
result = mtmp;
}
}
// matrix sum
template<class M>
void testMad_C(const M & m1, const M & m2, double & time, M & result) {
M mtmp = m2;
test::Timer t(time,"M+M ");;
for (int l = 0; l < NLOOP; l++)
{
mtmp(1,1) = gV[l];
result = m1; result += mtmp;
}
}
// matrix * constant
template<class M>
void testMscale_C(const M & m1, double a, double & time, M & result) {
M mtmp = m1;
test::Timer t(time,"a*M ");;
for (int l = 0; l < NLOOP; l++)
{
mtmp(1,1) = gV[l];
result = mtmp * a;
}
}
// clhep matrix matrix op (index starts from 1)
template<class A, class B, class C>
void testMM_C(const A & a, const B & b, const C & c, double & time, C & result) {
B btmp = b;
test::Timer t(time,"M*M ");
for (int l = 0; l < NLOOP; l++)
{
btmp(1,1) = gV[l];
result = a * btmp + c;
}
}
//inversion //inversion
template<class M> template<class M>
void testInv_C( const M & a, double & time, M& result){ void testInv_C( const M & a, double & time, M& result){
M mtmp = a;
test::Timer t(time,"inv "); test::Timer t(time,"inv ");
int ifail = 0; int ifail = 0;
for (int l = 0; l < NLOOP; l++) for (int l = 0; l < NLOOP; l++)
{ {
result = a.inverse(ifail); mtmp(1,1) = gV[l];
result = mtmp.inverse(ifail);
if (ifail) {std::cout <<"error inverting" << mtmp << std::endl; return; }
} }
} }
// general matrix matrix op // general matrix matrix op
template<class A, class B, class C> template<class A, class B, class C>
void testATBA_C(const A & a, const B & b, double & time, C & result) { void testATBA_C(const A & a, const B & b, double & time, C & result) {
B btmp = b;
test::Timer t(time,"At*M*A"); test::Timer t(time,"At*M*A");
for (int l = 0; l < NLOOP; l++) for (int l = 0; l < NLOOP; l++)
{ {
btmp(1,1) = gV[l];
//result = a.T() * b * a; //result = a.T() * b * a;
result = a * b * a.T(); result = a * btmp * a.T();
} }
} }
template<class A, class B, class C> template<class A, class B, class C>
void testATBA_C2(const A & a, const B & b, double & time, C & result) { void testATBA_C2(const A & a, const B & b, double & time, C & result) {
B btmp = b;
test::Timer t(time,"At*M*A"); test::Timer t(time,"At*M*A");
for (int l = 0; l < NLOOP; l++) for (int l = 0; l < NLOOP; l++)
{ {
result = b.similarity(a); btmp(1,1) = gV[l];
result = btmp.similarity(a);
}
}
template<class A, class C>
void testMT_C(const A & a, double & time, C & result) {
A atmp = a;
test::Timer t(time,"Transp");
for (int l = 0; l < NLOOP; l++)
{
atmp(1,1) = gV[l];
result = atmp.T();
} }
} }
......
smatrix/test/testOperations.cxx
+ 54
− 12
View file @ 8dc3fbcf
#define ENABLE_TEMPORARIES
#include "Math/SVector.h" #include "Math/SVector.h"
#include "Math/SMatrix.h" #include "Math/SMatrix.h"
...@@ -356,12 +357,15 @@ int test_tmatrix_op() { ...@@ -356,12 +357,15 @@ int test_tmatrix_op() {
MnMatrix C1(NDIM1,NDIM1); testATBA_T(B,C0,t_ama,C1); MnMatrix C1(NDIM1,NDIM1); testATBA_T(B,C0,t_ama,C1);
//if (k == 0) C1.Print(); //if (k == 0) C1.Print();
MnMatrix C2(NDIM1,NDIM1); testInv_T(C1,t_inv,C2); MnMatrix C2(NDIM1,NDIM1); testInv_T(C1,t_inv,C2);
//if (k == 0) C2.Print();
MnVector v3(NDIM1); testVeq(v,t_veq,v3); MnVector v3(NDIM1); testVeq(v,t_veq,v3);
MnVector v4(NDIM1); testVad_T(v2,v3,t_vad,v4); MnVector v4(NDIM1); testVad_T(v2,v3,t_vad,v4);
MnVector v5(NDIM1); testVscale_T(v4,2.0,t_vsc,v5); MnVector v5(NDIM1); testVscale_T(v4,2.0,t_vsc,v5);
MnMatrix C3(NDIM1,NDIM1); testMeq(C,t_meq,C3); MnMatrix C3(NDIM1,NDIM1); testMeq(C,t_meq,C3);
MnMatrix C4(NDIM1,NDIM1); testMad_T(C2,C3,t_mad,C4); MnMatrix C4(NDIM1,NDIM1); testMad_T(C2,C3,t_mad,C4);
//if (k == 0) C4.Print();
MnMatrix C5(NDIM1,NDIM1); testMscale_T(C4,0.5,t_msc,C5); MnMatrix C5(NDIM1,NDIM1); testMscale_T(C4,0.5,t_msc,C5);
//if (k == 0) C5.Print();
MnMatrix C6(NDIM1,NDIM1); testMT_T(C5,t_tra,C6); MnMatrix C6(NDIM1,NDIM1); testMT_T(C5,t_tra,C6);
#ifdef DEBUG #ifdef DEBUG
...@@ -508,8 +512,10 @@ int test_hepmatrix_op() { ...@@ -508,8 +512,10 @@ int test_hepmatrix_op() {
std::cout << " HepMatrix operations test " << first << " x " << second << std::endl; std::cout << " HepMatrix operations test " << first << " x " << second << std::endl;
std::cout << "************************************************\n"; std::cout << "************************************************\n";
double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_vsc, t_msc, t_ama = 0; double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_vsc, t_msc, t_ama, t_tra = 0;
double totTime1, totTime2;
double r1,r2; double r1,r2;
int npass = NITER; int npass = NITER;
...@@ -525,6 +531,8 @@ int test_hepmatrix_op() { ...@@ -525,6 +531,8 @@ int test_hepmatrix_op() {
MnVector v(NDIM1); MnVector v(NDIM1);
MnVector p(NDIM2); MnVector p(NDIM2);
TStopwatch w;
{ {
// fill matrices with random data // fill matrices with random data
fillRandomMat(r,A,first,second,1); fillRandomMat(r,A,first,second,1);
...@@ -548,30 +556,46 @@ int test_hepmatrix_op() { ...@@ -548,30 +556,46 @@ int test_hepmatrix_op() {
} }
#endif #endif
w.Start();
MnVector v1(NDIM1); testMV(A,v,t_mv,v1); MnVector v1(NDIM1); testMV(A,v,t_mv,v1);
MnVector v2(NDIM1); testGMV(A,v,v1,t_gmv,v2); MnVector v2(NDIM1); testGMV(A,v,v1,t_gmv,v2);
MnMatrix C0(NDIM1,NDIM1); testMM(A,B,C,t_mm,C0); MnMatrix C0(NDIM1,NDIM1); testMM_C(A,B,C,t_mm,C0);
MnMatrix C1(NDIM1,NDIM1); testATBA_C(B,C0,t_ama,C1); MnMatrix C1(NDIM1,NDIM1); testATBA_C(B,C0,t_ama,C1);
//std::cout << " C1 = " << C1 << std::endl; //std::cout << " C1 = " << C1 << std::endl;
MnMatrix C2(NDIM1,NDIM1); testInv_C(C1,t_inv,C2); MnMatrix C2(NDIM1,NDIM1); testInv_C(C1,t_inv,C2);
//std::cout << " C2 = " << C2 << std::endl;
MnVector v3(NDIM1); testVeq(v,t_veq,v3); MnVector v3(NDIM1); testVeq(v,t_veq,v3);
MnVector v4(NDIM1); testVad(v2,v3,t_vad,v4); MnVector v4(NDIM1); testVad(v2,v3,t_vad,v4);
MnVector v5(NDIM1); testVscale(v4,2.0,t_vsc,v5); MnVector v5(NDIM1); testVscale(v4,2.0,t_vsc,v5);
MnMatrix C3(NDIM1,NDIM1); testMeq(C,t_meq,C3); MnMatrix C3(NDIM1,NDIM1); testMeq_C(C,t_meq,C3);
MnMatrix C4(NDIM1,NDIM1); testMad(C2,C3,t_mad,C4); MnMatrix C4(NDIM1,NDIM1); testMad_C(C2,C3,t_mad,C4);
MnMatrix C5(NDIM1,NDIM1); testMscale(C4,0.5,t_msc,C5); //std::cout << " C4 = " << C4 << std::endl;
MnMatrix C5(NDIM1,NDIM1); testMscale_C(C4,0.5,t_msc,C5);
//std::cout << " C5 = " << C5 << std::endl;
MnMatrix C6(NDIM1,NDIM1); testMT_C(C5,t_tra,C6);
r1 = testDot_C(v3,v5,t_dot); r1 = testDot_C(v3,v5,t_dot);
r2 = testInnerProd_C(C5,v5,t_prd); r2 = testInnerProd_C(C6,v5,t_prd);
#ifdef DEBUG
if (k == 0) {
std::cout << " C6 = " << C6 << std::endl;
std::cout << " v5 = " << v5 << std::endl;
}
#endif
// MnMatrix C2b(NDIM1,NDIM1); testInv_T2(C1,t_inv2,C2b); // MnMatrix C2b(NDIM1,NDIM1); testInv_T2(C1,t_inv2,C2b);
w.Stop();
totTime1 = w.RealTime();
totTime2 = w.CpuTime();
} }
// tr.dump(); // tr.dump();
std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
std::cout << " r1 = " << r1 << " r2 = " << r2 << std::endl; std::cout << " r1 = " << r1 << " r2 = " << r2 << std::endl;
return 0; return 0;
...@@ -601,6 +625,7 @@ int test_hepmatrix_sym_op() { ...@@ -601,6 +625,7 @@ int test_hepmatrix_sym_op() {
double t_meq, t_mad, t_mv, t_gmv, t_mm, t_prd, t_inv, t_msc, t_ama = 0; double t_meq, t_mad, t_mv, t_gmv, t_mm, t_prd, t_inv, t_msc, t_ama = 0;
double totTime1, totTime2;
double r1; double r1;
...@@ -615,6 +640,8 @@ int test_hepmatrix_sym_op() { ...@@ -615,6 +640,8 @@ int test_hepmatrix_sym_op() {
MnVector v(NDIM1); MnVector v(NDIM1);
#define N NDIM1 #define N NDIM1
TStopwatch w;
{ {
// fill matrices with random data // fill matrices with random data
fillRandomSym(r,A,first,1); fillRandomSym(r,A,first,1);
...@@ -631,15 +658,16 @@ int test_hepmatrix_sym_op() { ...@@ -631,15 +658,16 @@ int test_hepmatrix_sym_op() {
} }
#endif #endif
w.Start();
MnVector v1(N); testMV(A,v,t_mv,v1); MnVector v1(N); testMV(A,v,t_mv,v1);
MnVector v2(N); testGMV(A,v,v1,t_gmv,v2); MnVector v2(N); testGMV(A,v,v1,t_gmv,v2);
MnMatrix C0(N,N); testMM(A,B,C,t_mm,C0); MnMatrix C0(N,N); testMM_C(A,B,C,t_mm,C0);
MnSymMatrix C1(N); testATBA_C2(C0,B,t_ama,C1); MnSymMatrix C1(N); testATBA_C2(C0,B,t_ama,C1);
MnSymMatrix C2(N); testInv_C(A,t_inv,C2); MnSymMatrix C2(N); testInv_C(A,t_inv,C2);
MnSymMatrix C3(N); testMeq(C2,t_meq,C3); MnSymMatrix C3(N); testMeq_C(C2,t_meq,C3);
MnSymMatrix C4(N); testMad(A,C3,t_mad,C4); MnSymMatrix C4(N); testMad_C(A,C3,t_mad,C4);
MnSymMatrix C5(N); testMscale(C4,0.5,t_msc,C5); MnSymMatrix C5(N); testMscale_C(C4,0.5,t_msc,C5);
r1 = testInnerProd_C(C5,v2,t_prd); r1 = testInnerProd_C(C5,v2,t_prd);
...@@ -648,10 +676,15 @@ int test_hepmatrix_sym_op() { ...@@ -648,10 +676,15 @@ int test_hepmatrix_sym_op() {
if (k == 0) { if (k == 0) {
} }
#endif #endif
w.Stop();
totTime1 = w.RealTime();
totTime2 = w.CpuTime();
} }
//tr.dump(); //tr.dump();
std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
std::cout << " r1 = " << r1 << std::endl; std::cout << " r1 = " << r1 << std::endl;
return 0; return 0;
...@@ -667,7 +700,7 @@ int test_hepmatrix_sym_op() { ...@@ -667,7 +700,7 @@ int test_hepmatrix_sym_op() {
MATRIX_SIZE=N; \ MATRIX_SIZE=N; \
TEST_TYPE=0; test_smatrix_op<N,N>(); \ TEST_TYPE=0; test_smatrix_op<N,N>(); \
TEST_TYPE=1; test_tmatrix_op<N,N>(); \ TEST_TYPE=1; test_tmatrix_op<N,N>(); \
TEST_TYPE=2; test_hepmatrix_op<N,N>(); TEST_TYPE=2; test_hepmatrix_op<N,N>(); \
TEST_TYPE=3; test_smatrix_sym_op<N,N>(); \ TEST_TYPE=3; test_smatrix_sym_op<N,N>(); \
TEST_TYPE=4; test_tmatrix_sym_op<N,N>(); \ TEST_TYPE=4; test_tmatrix_sym_op<N,N>(); \
TEST_TYPE=5; test_hepmatrix_sym_op<N,N>(); TEST_TYPE=5; test_hepmatrix_sym_op<N,N>();
...@@ -745,9 +778,18 @@ int main(int argc , char *argv[] ) { ...@@ -745,9 +778,18 @@ int main(int argc , char *argv[] ) {
typeNames[0] = "SMatrix"; typeNames[0] = "SMatrix";
typeNames[1] = "TMatrix"; typeNames[1] = "TMatrix";
//typeNames[2] = "HepMatrix"; #if !defined(HAVE_CLHEP) && defined (TEST_SYM)
typeNames[2] = "SMatrix_sym"; typeNames[2] = "SMatrix_sym";
typeNames[3] = "TMatrix_sym"; typeNames[3] = "TMatrix_sym";
#elif defined(HAVE_CLHEP) && defined (TEST_SYM)
typeNames[2] = "HepMatrix";
typeNames[3] = "SMatrix_sym";
typeNames[4] = "TMatrix_sym";
typeNames[5] = "HepMatrix_sym";
#elif defined(HAVE_CLHEP) && !defined (TEST_SYM)
typeNames[2] = "HepMatrix";
#endif
#endif #endif
#ifndef TEST_ALL_MATRIX_SIZES #ifndef TEST_ALL_MATRIX_SIZES
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment