Hi i am new to NR.
when i try to use amoeba function written in NR C book

i discover the function do not change initial guess matrix p through pointer **p

all trying do not change matrix p which is stored initial parameters

1. where is wrong? as i have a number of years not writing c language, i can not find where is wrong.

2. some initial guess happen NMAX exceeded. How to avoid?


#include "stdafx.h"
/*
#include "../include/tbb/task_scheduler_init.h"
#include "../include/tbb/parallel_invoke.h"
#include "../include/tbb/parallel_reduce.h"
#include "../include/tbb/blocked_range2d.h"

#include "../include/tbb/parallel_for.h"
#include "../include/tbb/tick_count.h"

#include "../boost/scoped_ptr.hpp"
*/
#include <string>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#include <cmath>
#include <complex>
//#define TINY 1.0e-10 //A small number.
#define TINY 1.0e-5 //A small number.
#define NMAX 10000 //Maximum allowed number of function evaluations.
#define GET_PSUM \
for (j=1;j<=ndim;j++) {\
for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];\
psum[j]=sum;}
#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
#define NR_END 1
//#define FTOL 1.0e-6
#define FTOL 1.0e-6
#define FREE_ARG char*

float *vector(long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
float amotry(float **p, float y[], float psum[], int ndim, float (*funk)(float []), int ihi, float fac);
void amoeba(float **p, float y[], int ndim, float ftol,float (*funk)(float []), int *nfunk);
float funk(float x[]);

float *vector(long nl, long nh)/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;

v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) printf("allocation failure in vector()");
return v-nl+NR_END;
}
void free_vector(float *v, long nl, long nh)/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
float **matrix(long nrl, long nrh, long ncl, long nch)
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;

/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) printf("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;

/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))) ;
if (!m[nrl]) printf("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

/* return pointer to array of pointers to rows */
return m;
}
void amoeba(float **p, float y[], int ndim, float ftol,float (*funk)(float []), int *nfunk)
{
//float amotry(float **p, float y[], float psum[], int ndim,
//float (*funk)(float []), int ihi, float fac);
int i,ihi,ilo,inhi,j,mpts=ndim+1;
float rtol,sum,swap,ysave,ytry,*psum;
psum=vector(1,ndim);
*nfunk=0;
GET_PSUM
/*
for (j=1;j<=ndim;j++) {
for (sum=0.0,i=1;i<=mpts;i++)
sum += p[i][j];
psum[j]=sum;
}*/
for (;;) {
ilo=1;
//First we must determine which point is the highest (worst), next-highest, and lowest
//(best), by looping over the points in the simplex.
ihi = y[1]>y[2] ? (inhi=2,1) : (inhi=1,2);
for (i=1;i<=mpts;i++) {
if (y[i] <= y[ilo]) ilo=i;
if (y[i] > y[ihi]) {
inhi=ihi;
ihi=i;
} else if (y[i] > y[inhi] && i != ihi) inhi=i;
}
rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo])+TINY);
//Compute the fractional range from highest to lowest and return if satisfactory.
if (rtol < ftol) { //If returning, put best point and value in slot 1.
SWAP(y[1],y[ilo])
for (i=1;i<=ndim;i++)
SWAP(p[1][i],p[ilo][i])
break;
}
if (*nfunk >= NMAX)
{
printf("NMAX exceeded");
while(true){}
}
*nfunk += 2;
//Begin a new iteration. First extrapolate by a factor −1 through the face of the simplex
//across from the high point, i.e., reflect the simplex from the high point.
ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0);
if (ytry <= y[ilo])
{
//Gives a result better than the best point, so try an additional extrapolation by a
//factor 2.
ytry=amotry(p,y,psum,ndim,funk,ihi,2.0);
}
else if (ytry >= y[inhi]) {
//The reflected point is worse than the second-highest, so look for an intermediate
//lower point, i.e., do a one-dimensional contraction.
ysave=y[ihi];
ytry=amotry(p,y,psum,ndim,funk,ihi,0.5);
if (ytry >= ysave) { //Can't seem to get rid of that high point. Better
for (i=1;i<=mpts;i++) { //contract around the lowest (best) point

if (i != ilo) {
for (j=1;j<=ndim;j++)
p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
y[i]=(*funk)(psum);
}
}
*nfunk += ndim; //Keep track of function evaluations.
//GET_PSUM Recompute psum.
}
} else
--(*nfunk); //Correct the evaluation count.
} //Go back for the test of doneness and the next
free_vector(psum,1,ndim); //iteration.
}
float amotry(float **p, float y[], float psum[], int ndim, float (*funk)(float []), int ihi, float fac)
{
int j;
float fac1,fac2,ytry,*ptry;
ptry=vector(1,ndim);
fac1=(1.0-fac)/ndim;
fac2=fac1-fac;
for (j=1;j<=ndim;j++)
{
printf("%d ptry[j] = %f " ,j, ptry[j]);
ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
}
ytry=(*funk)(ptry); //Evaluate the function at the trial point.
if (ytry < y[ihi]) { //If it's better than the highest, then replace the highest.
y[ihi]=ytry;
for (j=1;j<=ndim;j++) {
psum[j] += ptry[j]-p[ihi][j];
p[ihi][j]=ptry[j];
}
}
free_vector(ptry,1,ndim);
return ytry;
}
float L;
float funk(float *x)
{
//printf("%f", L);
return (x[1]*-L*x[1]/3)/x[1];
}
int _tmain(int argc, _TCHAR* argv[])
{
L = 0.278327923;
int nfunc;
int MP = 3;
int NP = 2;
int ndim = NP;
float **p;
float *x;
float *y;
x = vector(1, NP);
y = vector(1, MP);
p = matrix(1, MP, 1, NP);

p[1][1] = 0.5;
p[1][2] = 0.4;
p[2][1] = 0.5;
p[2][2] = 0.4;
p[3][1] = 0.5;
p[3][2] = 0.4;
/*
p[1][1] = 0.5;
p[1][2] = 0.3;
p[2][1] = 0.3;
p[2][2] = 0.5;
p[3][1] = 0.5;
p[3][2] = 0.5;
*/
float pt[3];
for (int i=1;i<=MP;i++) {
//for (int j=1;j<=NP;j++)
//{
pt[1] = p[i][1];
pt[2] = p[i][2];
//x[j]=p[i][j]=(i == (j+1) ? 0.5 : 0.0);
//}
y[i]=funk(pt);
}

printf("y[1] %f", y[1]);
printf("\n");
printf("y[2] %f", y[2]);
printf("\n");
printf("y[3] %f", y[3]);
printf("\n");
//amoeba(float **p, float y[], int ndim, float ftol,float (*funk)(float []), int *nfunk)
amoeba(p, y, ndim, FTOL,funk, &nfunc);

for (int i=1;i<=MP;i++) {
printf("%3d ",i);
for (int j=1;j<=NP;j++)
{
printf("%12.6f ",p[i][j]);
}
printf("y: ");
printf("%12.6f\n",y[i]);
}

while(true){}
return 0;
}