Optional return function

Dear co-rooters,

I am writing a function/macro that will be either used standalone or will be called from another function/macro.
This function makes some matrix computations using ublas and should either return a matrix or nothing (for instance a zero matrix).

The way I thought it would work, was to include a bool argument and if it’s true a matrix will be returned otherwise another zero matrix would be the output using the following check at the end of the function

if (return_TMatrix) return COV_USER;
else {zero_matrix<float> m_ZERO(1, 1); return m_ZERO;}

The thing is that when I run the code it breaks… I guess I am trying to do something strange?

Any idea or advice?

The code can be found below

//#include "/afs/cern.ch/user/a/astamato/public/lib_thanos/c_header.h"
//#include "/afs/cern.ch/user/a/astamato/public/lib_thanos/root_header.h"

// source /afs/cern.ch/sw/lcg/external/gcc/4.8/x86_64-slc6/setup.csh
// source /afs/cern.ch/sw/lcg/app/releases/ROOT/6.02.02/x86_64-slc6-gcc48-opt/root/bin/thisroot.csh

//ROOT6 on lxplus7
// source /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.12.06/x86_64-centos7-gcc48-opt/root/bin/thisroot.csh

#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>// std::max()

#include "TROOT.h"
#include "TGraph.h"
#include "TH1.h"
#include "TMath.h"
#include "TMatrixF.h"
#include "TMatrixFSym.h"
#include "TMatrixFSparse.h"
#include "TMatrixTSym.h"
#include "TString.h"
#include "TRandom.h"

#include </usr/include/boost/numeric/ublas/matrix.hpp>
#include </usr/include/boost/numeric/ublas/matrix_sparse.hpp>
#include </usr/include/boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include </usr/include/boost/numeric/ublas/io.hpp>

#include "covariance_lib.cc"

using namespace std;
using namespace boost::numeric::ublas;

// covariance_linear_interpolation("ENDF_B-VII_1_reduced.dat", -2, 0, 100, 0)

boost::numeric::ublas::symmetric_matrix<float> covariance_linear_interpolation(TString filename, int E_MIN, int E_MAX, int N_BPDEC, bool return_TMatrix){
// (1) Define the isolethargic binning - MeV
	int 	ndec  = E_MAX - E_MIN;
	int 	nbins = (int) ndec*N_BPDEC;
	float	step  = (float) ndec / nbins;
	float*	En_USER = new float[nbins];

	for(int i=0; i < nbins; i++) {
		En_USER[i] = (float) pow(10., step * (float) i + E_MIN);

// (2) Read the evaluation and interpolate
	int nbins_EVAL = count_lines(filename);
	float* xs_USER = new float[nbins];
	float* En_EVAL = new float[nbins_EVAL];
	float* xs_EVAL = new float[nbins_EVAL];
	int*   idx     = new int[nbins];
	input_ascii(filename, En_EVAL, xs_EVAL);	
	linear_interpolation(nbins, nbins_EVAL, En_USER, xs_USER, En_EVAL, xs_EVAL, idx);

// (3) Define the size for the covariance matrices
	const int size_EVAL = nbins_EVAL;
	const int size_USER = nbins;
// (4) Define the types of T-objects that will be used
	float *covariance_EVAL = new float[size_EVAL*(size_EVAL+1)/2];
	for (int i=0; i<size_EVAL*(size_EVAL+1)/2; ++i)
		covariance_EVAL[i] = 1.e-3*i;

	symmetric_matrix<float, upper> COV_EVAL (size_EVAL, size_EVAL);
	COV_EVAL = fill_symmetric( COV_EVAL, covariance_EVAL, size_EVAL);	

	delete[] covariance_EVAL;
	delete[] xs_USER;
	delete[] xs_EVAL;

// (5) Define the non-zero elements of the SParse TMatrix
	int*   icol  = new int[2*(nbins)];
	int*   irow  = new int[2*(nbins)];
	float* coeff = new float[2*(nbins)];

	for (int i=0; i<(nbins); ++i){
		icol[2*i]    = idx[i];
		icol[2*i+1]  = idx[i] + 1;	
		irow[2*i]    = i;
		irow[2*i+1]  = i;
		coeff[2*i]   = 1. + ( En_USER[i] - En_EVAL[ idx[i] ] )/( En_EVAL[idx[i]]-En_EVAL[idx[i+1]] );
		coeff[2*i+1] = coeff[2*i] - 1;
		if ( !isfinite(coeff[2*i]) ){// If En_USER == En_EVAL then sensitivity coefficient is 1
			coeff[2*i]   = 1;
			coeff[2*i+1] = 0;
	delete[] En_USER;
	delete[] En_EVAL;
	delete[] idx;
	mapped_matrix<float> G(size_USER, size_EVAL,2*nbins);
	G = fill_sparse( G, 2*nbins, irow, icol, coeff);
	delete[] icol;
	delete[] irow;
	delete[] coeff;

// (6) Calculate the covariance matrix - will be symmetric
	matrix<float> COV_USER;
	COV_USER = prod( G,COV_EVAL );
	COV_USER = prod( COV_USER, trans(G) );
	if (return_TMatrix) return COV_USER;
	else {zero_matrix<float> m_ZERO(1, 1); return m_ZERO;}


// 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

void input_ascii(TString filename, float* X_EVAL, float* Y_EVAL){

	std::ifstream myfile;
	float x = 0, y = 0;
	int   i = 0;
		myfile >> x >> y;
		if (!myfile.good()) break;
		X_EVAL[i]    = x;
		Y_EVAL[i]    = y;

int count_lines(TString filename){

	float	x, y;
	int nlines;	
	// (1) Open the file to count its lines
	std::ifstream myfile;
		cout << "File not found" << endl;
  	nlines = 0;
  		myfile >> x >> y;
  		if (!myfile.good()) break;
	return nlines;


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;


boost::numeric::ublas::symmetric_matrix<float> fill_symmetric (boost::numeric::ublas::symmetric_matrix<float> m_sym, float* filler, int size_FILLER){

	float* in = filler;
	for (size_t i=0; i<m_sym.size1(); ++ i)
    		for (size_t j = 0; j <= i && in != &filler[size_FILLER+1]; ++ j)
        		m_sym (i, j) = *in++;
	return m_sym;


boost::numeric::ublas::mapped_matrix<float> fill_sparse (boost::numeric::ublas::mapped_matrix<float> m_spar, int size_DATA, int* row, int* col, float* data){

	for (int i=0; i<size_DATA; ++ i){
		m_spar( row[i], col[i] ) = data[i];}
	return m_spar;


boost::numeric::ublas::diagonal_matrix<float> fill_diagonal (boost::numeric::ublas::diagonal_matrix<float> m_diag, float* data) {
    for (int i = 0; i < (int) m_diag.size1(); ++ i)
    	m_diag (i, i) = data[i];

    return m_diag;

A sample file can be found here. I am running ROOT6.12.04 on lxplus7 and compiling the code in ACLiC using

.L covariance_linear-interpolation.C++

The code can be run using the following

covariance_linear_interpolation("ENDF_B-VII_1_reduced.dat", -2, 0, 100, 0)

In your “covariance_lib.cc”, at least:

  for(int i = 0; i < size_USER; i++) {
    boo = true;
    for (int j = 0; j < size_EVAL; j++) {

BTW. “Invalid read of size 4/8”, when related to a “float” or “int” or “double” or “long int” array, usually means that you are trying to access array’s element with an improper index (e.g. beyond the end of the array).

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.