#include <math.h>
#include <stdio.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 200

void dpowell(double p[], double **xi, int n, double ftol, int *iter, double *fret,
	double (*func)(double []))
{
	void dlinmin(double p[], double xi[], int n, double *fret,
		double (*func)(double []));
	int i,ibig,j;
	double del,fp,fptt,t,*pt,*ptt,*xit;

	pt=dvector(1,n);
	ptt=dvector(1,n);
	xit=dvector(1,n);
	*fret=(*func)(p);
//	printf("in powell, the function value: %20.12f\n", *fret);
	for (j=1;j<=n;j++) pt[j]=p[j];
	for (*iter=1;;++(*iter)) {
		fp=(*fret);
		ibig=0;
		del=0.0;
		for (i=1;i<=n;i++) {
			for (j=1;j<=n;j++) xit[j]=xi[j][i];
			fptt=(*fret);
			dlinmin(p,xit,n,fret,func);
//			double last1=1.0;
//			for(int k=1;k<=n;k++) {last1-=p[k]; printf("%6.0f ", p[k]*100000);}
//			printf("%6.0f val:%12.8f\n",last1*100000,*fret);
//			printf("new fret: %20.12f  check value: %20.12f\n",*fret,(*func)(p));
			if (fabs(fptt-(*fret)) > del) {
				del=fabs(fptt-(*fret));
//			if (fptt-(*fret) > del) {
//				del=fptt-(*fret);
				ibig=i;
			}
		}
//		printf("converge criteria: %20.15e\n", (2.0*fabs(fp-(*fret))/(fabs(fp)+fabs(*fret))));
		if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
//		if (2.0*(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
			free_dvector(xit,1,n);
			free_dvector(ptt,1,n);
			free_dvector(pt,1,n);
			return;
		}
		if (*iter == ITMAX) {//nrerror("powell exceeding maximum iterations.");
			printf("powell exceeding maximum iterations.\n");
			free_dvector(xit,1,n);
			free_dvector(ptt,1,n);
			free_dvector(pt,1,n);
			return;
		}

		for (j=1;j<=n;j++) {
			ptt[j]=2.0*p[j]-pt[j];
			xit[j]=p[j]-pt[j];
			pt[j]=p[j];
		}
		fptt=(*func)(ptt);
//		printf("fptt: %20.12f  fp: %20.12f\n", fptt,fp);
		if (fptt < fp) {
			t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
			if (t < 0.0) {
				dlinmin(p,xit,n,fret,func);
//				printf("fret in the last part: %20.12f check: %20.12f\n",*fret, (*func)(p));
				for (j=1;j<=n;j++) {
					xi[j][ibig]=xi[j][n];
					xi[j][n]=xit[j];
				}
			}
		}
		for(i=1;i<=n;i++) printf("%6.0f ",p[i]*100000);
		double last=1.0; for(i=1;i<=n;i++) last-=p[i];
		printf("%6.0f ",last*100000);
		printf(" Val: %20.12f and fret: %20.12f\n\n",(*func)(p), *fret);
		fflush(stdout);
	}
}
#undef ITMAX
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */
