still segmentation fault in my SVD code

[Greg]

The code where the error occurs is below in svdfit. This is a
difficult one becuse changing the reporting (cout) can move the error.
dblvec and dblmat are std::vectors of doubles. I am using the
template below for debuging these vectors, but it is not showing
anything useful I can see. I output the dimensions of each vec to show
what it is, is the same as what it should be. Multifunc returns a
vector afunc (I can post more details) I checked the dimensions of the
other objects x, y etc. I feel there is somethin I doing wrong with
C++ data types rather than a algebraic syntax type mistake (like a bad
dimension). Any help would be appreciated.

void ffit::testfnfit(dblmat &x, dblvec &y, int ndata, dblvec &a, int
ma,
int max_iters, double &chisq) // , idxmultifunc *funcs)
{
DEBUG_LOG("test..");
int i;
dblvec sig, w;
sig.reserve(ndata);
sig.resize(ndata);
w.reserve(ma);
w.resize(ma);
dblmat u, v;
u.reserve(ndata);
u.resize(ndata);
for (i = 0; i < ma; i++)
u.resize(ma);
v.resize(ma);
for (i = 0; i < ma; i++)
v.resize(ma);
// DEBUG_LOG("test6");
for (i = 0; i < ndata; i++)
sig = 1.0;
// DEBUG_LOG();
std::cout << "calling svd " << ma << std::endl;
svdfit(x, y, sig, ndata, a, u, v, w, ma, max_iters, chisq); // ,
funcs);
}

void ffit::svdfit(dblmat &x, dblvec &y, dblvec &sig, const int &ndata,
dblvec &a, dblmat &u, dblmat &v, dblvec &w, const int &ma,
const int &max_iters, double &chisq)
{
int i, j;
dblvec afunc, b;
double tmp, wmax, thresh, s;
bool conv;
double const tol = 1.0E-5;

DEBUG_LOG("SVDFit");

afunc.reserve(ma);
afunc.resize(ma);
b.reserve(ndata);
b.resize(ndata);

std::cout << "sig " << sig.size() << "~" << ndata << std::endl;
std::cout << "u " << u.size() << "~" << ndata << " " << u[0].size() <<
"~" << ma << std::endl;
std::cout << "a " << a.size() << "~" << ma << std::endl;
std::cout << "y " << y.size() << "~" << ndata << std::endl;
for (i = 0; i < ndata; i++)
{
// funcs(i, ma, x, afunc);
multifn(i, ma, x, afunc);
DEBUG_LOG("SVDFits");
std::cout << i << " sig " << sig << " ";
tmp = 1/sig;
std::cout << "afsig " << tmp << " " << afunc[j] << std::endl;
for (j = 0; j < ma; j++)
u[j] = afunc[j]*tmp;
std::cout << "afsig i=" << i << std::endl;
// for (j = 0; j < ma; j++)
// std::cout << "a " << j << " " << afunc[j] << std::endl;
std::cout << "y " << y << std::endl;
std::cout << "afsig2" << std::endl;
b = y*tmp;
}
DEBUG_LOG("b4svdcmp");
svdcmp(u, ndata, ma, max_iters, w, v, conv);
std::cout << "w " << w.size() << " " << ma << std::endl;
wmax = 0.0;
for (j = 0; j < ma; j++)
{
if (w[j] > wmax)
wmax = w[j];
}
std::cout << "w " << w.size() << " " << ma << std::endl;
thresh = tol*wmax;
for (j = 0; j < ma; j++)
{
if (w[j] < thresh)
w[j] = 0.0;
}
std::cout << "a " << a.size() << std::endl;
svdksb(u, w, v, ndata, ma, b, a);
DEBUG_LOG("svdksb end");
chisq = 0.0;
for (i = 0; i < ndata; i++)
{
multifn(i, ma, x, afunc);
s = 0.0;
for (j = 0; j < ma; j++)
s = s + a[j]*afunc[j];
tmp = (y - s)/sig;
chisq = chisq + tmp*tmp;
}
DEBUG_LOG("svdfit end");
}

template < typename T >
class Vec : public std::vector< T > {
public:
Vec(): std::vector<T>() { }
Vec( int s ) : std::vector<T>(s) { }
T& operator[](int i)
{
return this -> at(i);
}
const T& operator[](int i) const
{
return this -> at(i);
}
};

typedef Vec<double> dblvec;
typedef Vec<dblvec> dblmat;


Output is as follows

[Session started at 2006-10-05 14:53:08 +0100.]
main
0
test fit
test..
calling svd 2
SVDFit
sig 10~10
u 10~10 2~2
a 2~2
y 10~10
fn 0 0 sz 1
fn 0 1 sz 1
multidone
SVDFits
terminate called after throwing an instance of 'std:ut_of_range'
what(): vector::_M_range_check
0 sig 1 afsig 1
nlps has exited due to signal 6 (SIGABRT).

 

[Greg]

tmp = 1/sig;

std::cout << "afsig " << tmp << " " << afunc[j] << std::endl;

Oops, sloppy debugging. There is one.
Now I'm getting an assertion error here...

for (j = 0; j < ma; j++)
u[j] = afunc[j]*tmp;

....it does not make sense it looks ok.

 

[jpalecek]

Greg napsal:

tmp = 1/sig;
std::cout << "afsig " << tmp << " " << afunc[j] << std::endl;

Oops, sloppy debugging. There is one.
Now I'm getting an assertion error here...

for (j = 0; j < ma; j++)
u[j] = afunc[j]*tmp;

...it does not make sense it looks ok.

The bug is in the test routine.

for (i = 0; i < ma; i++)
u.resize(ma);

should be

for (i = 0; i < ndata; i++)

Regards
Jiri Palecek

 

[Greg]

The bug is in the test routine.

for (i = 0; i < ma; i++)
u.resize(ma);

should be

for (i = 0; i < ndata; i++)

Regards
Jiri Palecek

Well done, I spotted it too believe it or not, although it has taken me
much longer than it has taken you. At least it was in the code I posted
and not lurking somewhere else. It tests ok now.

 

[BobR]

Greg wrote in message

[snip]
typedef Vec<double> dblvec;
typedef Vec<dblvec> dblmat;

void ffit::testfnfit(dblmat &x, dblvec &y, int ndata, dblvec &a, int
ma, int max_iters, double &chisq) // , idxmultifunc *funcs){
DEBUG_LOG("test..");
int i;
dblvec sig, w;
sig.reserve(ndata);
sig.resize(ndata);
w.reserve(ma);
w.resize(ma);

Hi Greg,
Just thought I'd rag on you a bit about using '.reserve()' and '.resize()',
one after the other.
In your case, you only need the '.resize()'[0]. You would use '.reserve()' if
you knew in advance that your vector would grow to a certain size, but
initially be smaller.

Check this out:

// includes here <iostream>, <ostream>, <vector>, etc.

void PrintVec( std::vector<int> const &vec, std:stream &sout){
sout<<" size="<<vec.size()<<" cap="<<vec.capacity()<<std::endl;
return;
} // PrintVec( vector<int> const &, std:stream &)

int main(){
using std::cout;
cout<<"\n--- VecInt(10) size test ---"<<std::endl;
std::vector<int> VecInt(10);
PrintVec( VecInt, cout);
VecInt.push_back( 1 );
PrintVec( VecInt, cout);
for(size_t i(0); i < 11; ++i){ VecInt.push_back( i );}
PrintVec( VecInt, cout);
VecInt.resize( 50 );
PrintVec( VecInt, cout);
VecInt.push_back( 1 );
PrintVec( VecInt, cout);
VecInt.resize( 40 );
PrintVec( VecInt, cout);
cout<<"--- VecInt(10) size test ---END"<<std::endl;
return 0;
} // main end

/* --- output ---
--- VecInt(10) size test ---
size=10 cap=10
size=11 cap=20
size=22 cap=40
size=50 cap=50
size=51 cap=100
size=40 cap=100
--- VecInt(10) size test ---END
*/

Note how '.resize(50)' set the capacity to 50 (instead of doubling the
capacity).
Then note that '.resize(40)' reduced the size, but not the capacity [1].

The way you did does not 'hurt' anything, but does add 'clutter', like doing:

int number(50); // initialize to 50
number = 50; // 'set' it to 50. (not needed, but also not an error)

[ if you knew all this, sorry. Maybe it will help some newbie. <G> ]

[0] - or initialize like Chris[val] showed you (that was you?) in other NG.
[1] - on my machine (win98, MinGW)