Dear rooters,
I am trying to perform matrix calculations, in order to compute covariance matrices but this requires big matrices (maybe in the order of 50000). I started playing with TMatrix to do so, and so far I have the following code
#include <iostream>
#include <fstream>
#include <vector>
#include "TROOT.h"
#include "TGraph.h"
#include "TH1.h"
#include "TMath.h"
#include "TMatrixF.h"
#include "TMatrixFSym.h"
#include "TMatrixTSym.h"
#include "TString.h"
#include "TRandom.h"
#include "covariance_lib.cc"
using namespace std;
int covariance_linear_interpolation(TString filename){
// (1) Free-Up memory
gROOT->Reset();
// (2) Define the isolethargic binning - MeV
int E_MAX = 1, E_MIN = -11, N_BPDEC = 1000;
int ndec = E_MAX - E_MIN;
int nbins = (int) ndec*N_BPDEC;
float step = (float) ndec / nbins;
float En_USER[nbins+1];
for(int i=0; i <= nbins; i++) {
En_USER[i] = (float) pow(10., step * (float) i + E_MIN);
}
// (3) Read the evaluation and interpolate
int nbins_EVAL = count_lines(filename);
float xs_USER[nbins+1], En_EVAL[nbins_EVAL], xs_EVAL[nbins_EVAL];
int idx[nbins+1];
input_ascii(filename, En_EVAL, xs_EVAL);
linear_interpolation(nbins, nbins_EVAL, En_USER, xs_USER, En_EVAL, xs_EVAL, idx);
// (2) Define the size for the covariance matrices
const int size_EVAL = nbins_EVAL;
const int size_USER = nbins;
// (3) Define the types of T-objects that will be used
float covariance_EVAL[size_EVAL*(size_EVAL+1)/2];
for (int i=0; i<size_EVAL*(size_EVAL+1)/2; ++i){
covariance_EVAL[i] = gRandom->Gaus(0.1, 5);
}
TMatrixFSym COV_EVAL(size_EVAL, covariance_EVAL);
//
//COV_EVAL.Print();
/*
// (5) Get the evaluated cross section and store it in a graph
float En_EVAL[size_EVAL], xs_EVAL[size_EVAL];
for (int i=0; i<size_EVAL; ++i){
En_EVAL[i] = i+1;
xs_EVAL[i] = cos(i);
}
TGraph *g_XS_EVAL = new TGraph(size_EVAL, En_EVAL, xs_EVAL);
g_XS_EVAL->Draw("APL");
// (6) Calculate the interpolated cross section and it in a TGraph as well
float En_USER[size_USER], xs_USER[size_USER], idx[size_USER];
//for (int i=0; i<size_USER; ++i){
// En_USER[i] = 1.2*(i+1);
// xs_USER[i] = g_XS_EVAL->Eval(En_USER[i]);
//}
TGraph *g_XS_USER = new TGraph(size_USER, En_USER, xs_USER);
g_XS_USER->SetLineColor(kBlue); g_XS_USER->Draw("PLsame");
// (6) Define the sensitivity matrix and fill it
SMatrixStd_USERxEVAL G;
for (int i=0; i<size_USER; ++i){
for (int j=0; j<size_EVAL; ++j){
G(i,j) = pow(0.02, i)*j+i/100;
}
}
// (7) Calculate the covariance matrix - will be symmetric
SMatrixStd_USER V_USER;
V_USER = G*V*Transpose(G);
// (8) Print the matrices for checking
std::cout << V << std::endl;
std::cout << "**************************************************************************\n";
std::cout << G << std::endl;
std::cout << "**************************************************************************\n";
std::cout << V_USER << std::endl;
*/
return 0;
}
The covariance_lib.cc
follows
// Inputs : -name = name of the ascii file to be read - two column file with x and y values, rescpectively assumed
// -X_USER = array that contains the user defined energy binning
// -output = if output = TRUE exports the interpolation array
// if output = FALSE exports the index of the input values used in the interpolation
void linear_interpolation(int size_USER, int size_EVAL, float* X_USER, float* Y_USER, float* X_EVAL, float* Y_EVAL, int* index){
float x1, x2,y1, y2, a, b;
bool boo;
x1 = 0.; x2 = 0.; y1 = 0; y2 = 0.; a = 0.; b = 0.;
for(int i=0; i<=size_USER; ++i){
boo = true;
for (int j=0; j<=size_EVAL; j++){
if (boo==true){
if(X_USER[i] == X_EVAL[j]){
Y_USER[i] = Y_EVAL[j];
index[i] = j;
boo = false;
}
else if(X_USER[i]>X_EVAL[j]){
x1 = X_EVAL[j];
y1 = Y_EVAL[j];
index[i] = j;
}
else if( X_USER[i]<X_EVAL[j] && boo == true ){
x2 = X_EVAL[j];
y2 = Y_EVAL[j];
a = (y2-y1)/(x2-x1);
b = y2-a*x2;
Y_USER[i] = a*X_USER[i]+b;
//index[i] = j;
boo = false;
}
}
}//end of loop over EVAL
}//end of loop over USER
}//____linear_interpolation()
void input_ascii(TString filename, float* X_EVAL, float* Y_EVAL){
ifstream myfile;
myfile.open(filename);
float x = 0, y = 0;
int i = 0;
while(1){
myfile >> x >> y;
if (!myfile.good()) break;
X_EVAL[i] = x;
Y_EVAL[i] = y;
i++;
}
myfile.close();
}//___input_ascii()
int count_lines(TString filename){
float x, y;
int nlines;
// (1) Open the file to count its lines
ifstream myfile;
myfile.open(filename);
if(!myfile){
cout << "File not found" << endl;
}
nlines = 0;
while(1){
myfile >> x >> y;
if (!myfile.good()) break;
nlines++;
}
myfile.close();
return nlines;
}//___count_lines()
std::vector<float> isolethargic_binning(int min, int max, int bpd){
int ndec = max - min;
int nbins = (int) ndec*bpd;
float step = (float) ndec / nbins;
std::vector<float> xbins(nbins+1);
for(int i=0; i <= nbins; i++) {
xbins[i] = (float) pow(10., step * (float) i + min);
}
return xbins;
}//___isolethargic_binning()
A sample file can be found here
For this sample file, I have to create a 20000 x 20000 symmetric matrix, but it seems that I am out of memory. So when I compile and run the code (note that I am running 5.32 on lxplus) using
root [0] .L covariance_linear-interpolation.C++
root [1] covariance_linear_interpolation("ENDF_B-VII_1.dat")
root automatically quits. I suppose that this is happening because the matrix must be big.
I also tried defining the array using
float *covariance_EVAL = new float[size_EVAL*(size_EVAL+1)/2];
and in this case the code is running for more time but ends up with a crash
*** Break *** segmentation violation
===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
#0 0x00007f444a16482e in waitpid () from /lib64/libc.so.6
#1 0x00007f444a0f6479 in do_system () from /lib64/libc.so.6
#2 0x00007f444ae07b04 in TUnixSystem::StackTrace() () from /usr/lib64/root/libCore.so.5.34
#3 0x00007f444ae06f23 in TUnixSystem::DispatchSignals(ESignals) () from /usr/lib64/root/libCore.so.5.34
#4 <signal handler called>
#5 0x00007f444a141aaf in memcpy () from /lib64/libc.so.6
#6 0x00007f4446464dd0 in TMatrixTBase<float>::SetMatrixArray(float const*, char const*) () from /usr/lib64/root/libMatrix.so
#7 0x00007f44464903b9 in TMatrixTSym<float>::SetMatrixArray(float const*, char const*) () from /usr/lib64/root/libMatrix.so
#8 0x00007f4446491627 in TMatrixTSym<float>::TMatrixTSym(int, float const*, char const*) () from /usr/lib64/root/libMatrix.so
#9 0x00007f4442a4b1f6 in covariance_linear_interpolation (filename=<incomplete type>) at /eos/user/a/astamato/covariance/code/./covariance_linear-interpolation.C:47
#10 0x00007f4442a4b34d in G__covariance_linear_interpolation_C_ACLiC_dict__0_2358 (result7=0x7ffc958a0350, funcname=<value optimized out>, libp=<value optimized out>, hash=<value optimized out>) at /eos/user/a/astamato/covariance/code/covariance_linear_interpolation_C_ACLiC_dict.cxx:82
#11 0x00007f4449785ba9 in Cint::G__ExceptionWrapper(int (*)(G__value*, char const*, G__param*, int), G__value*, char*, G__param*, int) () from /usr/lib64/root/libCint.so.5.34
#12 0x00007f444982ded1 in G__execute_call () from /usr/lib64/root/libCint.so.5.34
#13 0x00007f444982ed22 in G__call_cppfunc () from /usr/lib64/root/libCint.so.5.34
#14 0x00007f444980da05 in G__interpret_func () from /usr/lib64/root/libCint.so.5.34
#15 0x00007f44497fc918 in G__getfunction () from /usr/lib64/root/libCint.so.5.34
#16 0x00007f44497db12e in G__getitem () from /usr/lib64/root/libCint.so.5.34
#17 0x00007f44497dfcb8 in G__getexpr () from /usr/lib64/root/libCint.so.5.34
#18 0x00007f444985abe7 in G__exec_statement () from /usr/lib64/root/libCint.so.5.34
#19 0x00007f44497c72d1 in ?? () from /usr/lib64/root/libCint.so.5.34
#20 0x00007f44497c75de in G__exec_tempfile_fp () from /usr/lib64/root/libCint.so.5.34
#21 0x00007f4449866645 in G__process_cmd () from /usr/lib64/root/libCint.so.5.34
#22 0x00007f444adc7c46 in TCint::ProcessLine(char const*, TInterpreter::EErrorCode*) () from /usr/lib64/root/libCore.so.5.34
#23 0x00007f444ad2b0e0 in TApplication::ProcessLine(char const*, bool, int*) () from /usr/lib64/root/libCore.so.5.34
#24 0x00007f444a97a05b in TRint::HandleTermInput() () from /usr/lib64/root/libRint.so.5.34
#25 0x00007f444ae04f8e in TUnixSystem::CheckDescriptors() () from /usr/lib64/root/libCore.so.5.34
#26 0x00007f444ae05273 in TUnixSystem::DispatchOneEvent(bool) () from /usr/lib64/root/libCore.so.5.34
#27 0x00007f444ad85bd6 in TSystem::InnerLoop() () from /usr/lib64/root/libCore.so.5.34
#28 0x00007f444ad8770b in TSystem::Run() () from /usr/lib64/root/libCore.so.5.34
#29 0x00007f444ad282df in TApplication::Run(bool) () from /usr/lib64/root/libCore.so.5.34
#30 0x00007f444a97b5b4 in TRint::Run(bool) () from /usr/lib64/root/libRint.so.5.34
#31 0x000000000040103c in main ()
===========================================================
The lines below might hint at the cause of the crash.
If they do not help you then please submit a bug report at
http://root.cern.ch/bugs. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#5 0x00007f444a141aaf in memcpy () from /lib64/libc.so.6
#6 0x00007f4446464dd0 in TMatrixTBase<float>::SetMatrixArray(float const*, char const*) () from /usr/lib64/root/libMatrix.so
#7 0x00007f44464903b9 in TMatrixTSym<float>::SetMatrixArray(float const*, char const*) () from /usr/lib64/root/libMatrix.so
#8 0x00007f4446491627 in TMatrixTSym<float>::TMatrixTSym(int, float const*, char const*) () from /usr/lib64/root/libMatrix.so
#9 0x00007f4442a4b1f6 in covariance_linear_interpolation (filename=<incomplete type>) at /eos/user/a/astamato/covariance/code/./covariance_linear-interpolation.C:47
===========================================================
Root > !!!Dictionary position not recovered because G__unloadfile() is used in a macro!!!
Any idea on how to see what’s happening and how to use TMatrix
for big matrices?
Thanks in advance!