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 <gsl/gsl_linalg.h>

 

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<std::vector<double> >&, const std::vector<double>&);

  void Reset(const std::vector<std::vector<double> >&, 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 <gsl/gsl_linalg.h>

#include <boost/smart_ptr/shared_array.hpp>

 

regressionPCA::regressionPCA(const std::vector<std::vector<double> >&data,

                             const std::vector<double>& Ydata) {

  Reset(data,Ydata);

}

 

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;

  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<std::vector<double> >::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;

};