mrqmin.cpp (fixed in v. 2.11)
Saul Teukolsky
10-08-2002, 01:06 PM
Version 2.10 of mrqmin.cpp uses a fixed value of ma, the dimension of the vector of fitting coefficients. If you try different fits in the same program, you can't increase ma because of the static variables declared in mrqmin. The fix is to use new and delete to define the static variables. This will be done in the next bug fix release, v. 2.11. A temporary workaround is just to rerun the program when you change ma.
Note that it is not safe simply to delete the keyword "static". This may actually work on some machines, but will fail (correctly) on most.
This problem does not occur in the Fortran or C routines.
David Wilkinson
10-16-2002, 08:46 AM
A solution for these static variables which does not require multiple new's (not exception safe) is
int ma=a.size();
static Mat_DP oneda(1, 1);
static Vec_DP atry(1), beta(1), da(1);
if(alambda < 0.0)
{
oneda=Mat_DP(ma, 1);
atry=Vec_DP(ma);
beta=Vec_DP(ma);
da=Vec_DP(ma);
//etc
}
if(atry.size() != ma)
nrerror("Initial call must have alambda < 0");
It might also be a good idea to free this memory on the final call
if(alambda == 0.0)
{
oneda=Mat_DP(1, 1);
atry=Vec_DP(1);
beta=Vec_DP(1);
da=Vec_DP(1);
// etc
}
I have used 1 rather than 0 as the default dimension because (due to the design of NR_Vec and NRMat) assignment from a NRVec or NRMat of zero size may not be well-defined on all compilers (it causes new[] to be called with a zero argument).
BTW, why is oneda a matrix?
David Wilkinson
Sherlock
10-20-2002, 09:57 PM
Good ideas.
As for the last question. I think oneda[][] is a matrix rather than a vector because it is passed as the second argument to the Recipe gaussj(), and that argument is of matrix type. In the more general case, the second argument of gaussj() is a collection of right-hand sde vectors composed into a matrix, but here there is only a single right-hand side vector.
kibitzer
10-21-2002, 08:31 AM
I believe there is a problem with the proposed alternative code. A statement like
atry=Vec_DP(ma);
will give a syntax error.