Hi everyone,

I am trying to call the amoeba routine to minimize a function of 54 variables. The program halts with a segmentation fault right after amoeba is called. By inserting print statements in the amoeba code, I see that the fault occurs during the GET_PSUM statement near the beginning. Has anyone else experienced this?

Thanks,
Daryl

...the fault occurs during the GET_PSUM...

Are you using C? (The reason that I ask, is that I don't think it would happen with the C++ version. Of course, I could be wrong.)

Anyhow...

I think that segfaults like this are most commonly caused by addressing memory outside the bounds of an array.

In the C version, GET_PSUM assumes (but has no way of checking) that your array size is [ndim+1 x ndim], where ndim is the third argument in the call to amoeba. How did you declare the matrix whose name is the first argument in your call to amoeba?

Regards,

Dave

Thanks for your response Dave. I'm using the C version of amoeba, and this is how I'm declaring my matrix:

int i,j;
float x,obj,var[nVar],p[nVar+1][nVar],init[nVar],y[nVar+1];
float ftol=0.01;
int nfunk;

for(i=0;i<nVar+1;i++)
{
for(j=0;j<nVar;j++)
{
x=rand()/((double)RAND_MAX+1);
x=x*2.0-1.0; /* uniformly between -1 and 1 */
p[i][j]=x;
init[j]=x;
}
y[i]=objective(init); /* initial simplex */
}

amoeba(p,y,nVar,ftol,objective,&nfunk);


Maybe the problem is in the way I'm passing the parameters to amoeba.

Daryl

Here's the prototype for amoeba():

void amoeba(float **p, float y[], int ndim, float ftol,
float (*funk)(float []), int *iter);

The first argument to amoeba must be "Pointer to pointer to float." This is not the same as a 2-D array. It's a pointer to pointer. See Footnote.

The second argument must be "Pointer to float." (That's the same as float [], and, in fact, as far as a C compiler is concerned, it could be the name of an array, but there is more to the story, as I describe below.)

The third is an int (your nVar)

The fourth is a float. (the tolerance).

The fifth is a pointer to a function that returns a float, The function itself has a single argument: a pointer to a float.

The last is a simple pointer to int.

Now, here's the biggie if you are going to use amoeba() or other Numerical Recipes functions that operate on vectors and matrices:

The pointers to float and pointer-to-pointer to float point to "Numerical Recipes" style vectors and matrices. They point to dynamically-allocated blocks of memory, and the pointer values have been adjusted so that accessing their elements must be done using index values that start at one, not zero.


If you are an experienced C programmer, you must suppress your antipathy toward having arrays starting with index value 1 and just live with it. (Or you could use the C++ versions of these functions, which use standard C++ vector objects that are like arrays that start with index zero. NR version 3 is all C++, and is even easier to use for C++ programmers, but takes a little getting used to for us mere mortals.)


Numerical recipes functions vector() and matrix() are designed to make it relatively easier to allocate storage for the stuff, but you still have to remember to make your loops start at 1, not zero.

Something like the following. I have commented out your array declarations and I show how to set things up with NR vector and matrix function calls:


#include <stdio.h>
#include <stdlib.h>
.
. other system includes
.
#define NRANSI
#include "nrutil.h"
#include "nr.h"
.
.
.
float objective(float x[]); /* designed to operate on NR-style vector */

int main()
{
.
.
.
int nVar = 54; /* doesn't have to be constant (but it could be) */
float x;
int i, j;
/*
* replace arrays with pointers and use NR functions to
* set them up
*/
/*float var[nVar], p[nVar + 1][nVar], init[nVar], y[nVar + 1]; */
float *var;
float **p;
float *init;
float *y;
var = vector(1, nVar);
p = matrix(1, nVar+1, 1, nVar);
init = vector(1, nVar);
y = vector(1, nVar+1);
.
.
.
for (i = 1; i <= nVar + 1; i++) {
for (j = 1; j <= nVar; j++) {
x = rand() ... whatever
p[i][j] = x;
init[j] = x;
}
y[i] = ... whatever
.
.
.
.
.
.
amoeba(p, y, nVar, ftol, objective, &nfunk);
.
.
.
free_vector(y, 1, nVar);
free_matrix(p, 1, nVar+1, 1, nVar);
free_vector(init, 1, nVar);
free_vector(y, 1, nVar+1);

return 0;
} /* end of main() */



Regards,

Dave

Footnote: Didn't you get compiler warnings when you used a 2-D array name as the first argument in your call to function amoeba()? If it said something like, "incompatible pointer type," that's what it was talking about. Depending on the compiler, the message may have been about "different levels of indirection," or "suspicious pointer conversion," or some such thing, but there should have been some indication of a problem there.

Regardless of that, however, there is no way the compiler could have known that the function required pointer values that would lead to one-based arrays rather than the native C zero-based arrays.

Ah, the infamous 0-1 allocation issue. Yes, the compiler did complain about incompatible pointer types. I made the changes you suggested and the program works fine now.

Thanks for your help Dave. I appreciate you taking the time to help me fix my code.

Daryl