Hi Ted,
On 03/25/2012 03:38 AM, Ted Byers wrote:
boost::shared_array<double> Btmp(new double[c]),
Ytmp(new double[r]),
Utmp(new double[n]),
Vtmp(new double[c*c]),
Stmp(new double[c]);
double* B = Btmp.get();
double* Y = Ytmp.get();
double* U = Utmp.get();
double* V = Vtmp.get();
double* S = Stmp.get();
I am a little bit astonished about these lines. The pointer B is only
valid as long as Btmp is in scope and this scope exits after your Reset
method has finished (and therefore the object Btmp will be deleted since
it is a object created on the stack). So after reset has finished,
B,Y,... are all dangling pointer, which should never be accessed. Doing
so may arise a core dump earlier or later. You have to store at least
one instance of the shared_array to have access to B, Y, ... .
As it seems, your first post to the list have the same problem. I found
this section in your code (maybe you should post an example without the
misleading comments, since the program won't compile if one will remove
the comments in the first example. This was very missleading in your
first post, so I was unable to see your problem there!
Maybe you should rewrite your class:
#ifndef REG_CLASS_H
#define REG_CLASS_H
#include <iosfwd>
#include <vector>
#include
class regressionPCA {
private:
unsigned int nrows, ncols;
regressionPCA(void) {}; // makes default construction impossible
// Add this to your class, holding the pointers as long as they were
needed!
boost::shared_array<double> B, Y, U, V, S;
public:
regressionPCA(const std::vector&, const
std::vector<double>&);
void Reset(const std::vector&, const
std::vector<double>&);
inline void setNrows(unsigned int v) { nrows = v;};
inline void setNcols(unsigned int v) { ncols = v;};
inline unsigned int getNrows(void) const { return nrows; };
inline unsigned int getNcols(void) const { return ncols; };
};
#endif
---------
#include "reg.class.h"
#include <iostream>
#include <algorithm>
#include
#include
regressionPCA::regressionPCA(const std::vector&data,
const std::vector<double>& Ydata) {
Reset(data,Ydata);
}
void regressionPCA::Reset(const std::vector&data,
const std::vector<double>& Ydata) {
unsigned int r(0),c(0),n(0);
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
r = data.size();
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
c = data[0].size();
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
setNrows(r);
setNcols(c);
n = r * c;
std::cout << "r: " << r << "\tc: " << c << "\tn: " << n << std::endl;
/* THIS CODE IN THIS COMMENT IS WRONG, THIS WILL CREATE A STACK OBJECT
ONLY!
Refering to your first post on the list, B is not the B of your class.
The B of your class is hidden by the B you have created on the stack here.
boost::shared_array<double> B(new double[c]); // wrong!
*/
// This is what you want to do:
B.reset( new double[c] );
Y.reset( new double[n] );
// ... do this with the rest too!
// Use the pointer with B.get() for external libs, this was
quite correct.
// double *bptr = U.get();
double *bptr = U;
std::vector::const_iterator it = data.begin(),
end = data.end();
while (it != end) {
bptr = std::copy(it->begin(), it->end(),bptr);
++it;
}
bptr = Y;
std::copy(Ydata.begin(),Ydata.end(),bptr);
gsl_matrix_view Um = gsl_matrix_view_array(U, getNrows(), getNcols());
gsl_vector_view Ym = gsl_vector_view_array(Y, getNrows());
gsl_vector_view Bm = gsl_vector_view_array(B, getNcols());
gsl_vector_view Sm = gsl_vector_view_array(S, getNcols());
gsl_matrix_view Vm = gsl_matrix_view_array(V, 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::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;
};
