The code is direct simple.
void regressionPCA::Reset(const std::vector<std::vector<double> >&data,
const std::vector<double>& Ydata) {
unsigned int r(0),c(0),n(0);
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
r = data.size();
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
c = data[0].size();
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
setNrows(r);
setNcols(c);
n = r * c;
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
B.reset(new double[c]);
Y.reset(new double[r]);
U.reset(new double[n]);
V.reset(new double[c*c]);
S.reset(new double[c]);
std::cout << "Flag 1" << std::endl << std::flush;
double* Btmp = B.get();
double* Ytmp = Y.get();
double* Utmp = U.get();
double* Vtmp = V.get();
double* Stmp = S.get();
double *bptr = Utmp;
std::vector<std::vector<double> >::const_iterator it = data.begin(), end =
data.end();
while (it != end) {
bptr = std::copy(it->begin(), it->end(),bptr);
++it;
}
bptr = Ytmp;
std::cout << "Flag 2" << std::endl << std::flush;
std::copy(Ydata.begin(),Ydata.end(),bptr);
std::cout << "Flag 3" << std::endl << std::flush;
gsl_matrix_view Um = gsl_matrix_view_array(Utmp, getNrows(), getNcols());
gsl_vector_view Ym = gsl_vector_view_array(Ytmp, getNrows());
gsl_vector_view Bm = gsl_vector_view_array(Btmp, getNcols());
gsl_vector_view Sm = gsl_vector_view_array(Stmp, getNcols());
gsl_matrix_view Vm = gsl_matrix_view_array(Vtmp, getNcols(), getNcols());
gsl_linalg_SV_decomp_jacobi(&Um.matrix,&Vm.matrix,&Sm.vector);
gsl_linalg_SV_solve(&Um.matrix,&Vm.matrix,&Sm.vector,&Ym.vector,&Bm.vector);
std::cout << std::endl << std::endl << "Sv = " <<
gsl_vector_get(&Sm.vector,0) << "\t" << gsl_vector_get(&Sm.vector,1) <<
std::endl << std::flush;
std::cout << std::endl << std::endl << "V = " << std::endl;
std::cout << "\t" << gsl_matrix_get(&Vm.matrix,0,0) << "\t" <<
gsl_matrix_get(&Vm.matrix,0,1) << std::endl;
std::cout << "\t" << gsl_matrix_get(&Vm.matrix,1,0) << "\t" <<
gsl_matrix_get(&Vm.matrix,1,1) << std::endl;
std::cout << std::endl << std::endl << "Beta = " <<
gsl_vector_get(&Bm.vector,0) << "\t" << gsl_vector_get(&Bm.vector,1) <<
std::endl << std::flush;
};
This code executes fine the first two or three times the function it
invoked, but it never gets past the third execution without crashing. Now,
I thought that maybe it was my use of the boost::shared_array that might be
at fault, but the following program proves that that is not the case (as it
runs to completion).
void regressionPCA::Reset(const std::vector<std::vector<double> >&data,
const std::vector<double>& Ydata) {
unsigned int r(0),c(0),n(0);
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
r = data.size();
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
c = data[0].size();
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
setNrows(r);
setNcols(c);
n = r * c;
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl <<
std::flush;
B.reset(new double[c]);
Y.reset(new double[r]);
U.reset(new double[n]);
V.reset(new double[c*c]);
S.reset(new double[c]);
std::cout << "Flag 1" << std::endl << std::flush;
double* Btmp = B.get();
double* Ytmp = Y.get();
double* Utmp = U.get();
double* Vtmp = V.get();
double* Stmp = S.get();
double *bptr = Utmp;
std::vector<std::vector<double> >::const_iterator it = data.begin(), end =
data.end();
while (it != end) {
bptr = std::copy(it->begin(), it->end(),bptr);
++it;
}
bptr = Ytmp;
std::cout << "Flag 2" << std::endl << std::flush;
std::copy(Ydata.begin(),Ydata.end(),bptr);
std::cout << "Flag 3" << std::endl << std::flush;
gsl_matrix_view Um = gsl_matrix_view_array(Utmp, getNrows(), getNcols());
gsl_vector_view Ym = gsl_vector_view_array(Ytmp, getNrows());
gsl_vector_view Bm = gsl_vector_view_array(Btmp, getNcols());
gsl_vector_view Sm = gsl_vector_view_array(Stmp, getNcols());
gsl_matrix_view Vm = gsl_matrix_view_array(Vtmp, getNcols(), getNcols());
gsl_linalg_SV_decomp_jacobi(&Um.matrix,&Vm.matrix,&Sm.vector);
gsl_linalg_SV_solve(&Um.matrix,&Vm.matrix,&Sm.vector,&Ym.vector,&Bm.vector);
std::cout << std::endl << std::endl << "Sv = " <<
gsl_vector_get(&Sm.vector,0) << "\t" << gsl_vector_get(&Sm.vector,1) <<
std::endl << std::flush;
std::cout << std::endl << std::endl << "V = " << std::endl;
std::cout << "\t" << gsl_matrix_get(&Vm.matrix,0,0) << "\t" <<
gsl_matrix_get(&Vm.matrix,0,1) << std::endl;
std::cout << "\t" << gsl_matrix_get(&Vm.matrix,1,0) << "\t" <<
gsl_matrix_get(&Vm.matrix,1,1) << std::endl;
std::cout << std::endl << std::endl << "Beta = " <<
gsl_vector_get(&Bm.vector,0) << "\t" << gsl_vector_get(&Bm.vector,1) <<
std::endl << std::flush;
};
So, I went further, and if I replace the shared_array by naked pointers, and
apply operator delete[] to them at the very end of the function, the same
problem arises. So what can possibly be wrong in this code, or is it the
case GSL is doing something odd behind the scenes that causes the core dump?
Cheers
Ted
