#include <stdio.h>
#include <math.h>

/*
CGA:
Okay, I reformatted this and added p_min, and p_max args.
and added debug_messages
*/

int debug_messages = 0;

#define TINY 1.0e-25
#define ITMAX 100
static float sqrarg;
#define SQR(a) (sqrarg=(a),sqrarg*sqrarg)

void powell(p,p_min,p_max,xi,n,ftol,iter,fret,func,messages)
     float p[],p_min[],p_max[],**xi,ftol,*fret,(*func)();
     int n,*iter,messages;
{
  int i,ibig,j;
  float t,fptt,fp,del;
  float *pt,*ptt,*xit,*vector();
  void linmin(),nrerror(),free_vector();

  debug_messages = messages;

  pt=vector(1,n);
  ptt=vector(1,n);
  xit=vector(1,n);
  *fret=(*func)(p);
  for (j=1;j<=n;j++)
    pt[j]=p[j];
  for (*iter=1;;(*iter)++)
    {
      if ( debug_messages )
	{
	  printf( "POWELL ITERATION %d; XI:\n", *iter );  /* CGA */
	  for ( j = 1; j <= n; j++ ) /* CGA */
	    { /* CGA */
	      for ( i = 1; i <= n; i ++ ) /* CGA */
		printf( "%f ", xi[j][i] ); /* CGA */
	      printf( "\n" ); /* CGA */
	    } /* CGA */
	}
      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);
	  if ( debug_messages )
	    {
	      printf( "SEARCHING ALONG DIRECTION %d:\n", i ); /*CGA*/
	      for ( j = 1; j <= n; j++ ) /* CGA */
		printf( "%f ", xit[j] ); /* CGA */
	      printf( "\n\n" ); /* CGA */
	    }
	  linmin(p,p_min,p_max,xit,n,fret,func);
	  if (fabs(fptt-(*fret)) > del)
	    {
	      del=fabs(fptt-(*fret));
	      ibig=i;
	    }
	}
      if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))+TINY)
	{
	  free_vector(xit,1,n);
	  free_vector(ptt,1,n);
	  free_vector(pt,1,n);
	  return;
	}
      if (*iter == ITMAX || n==1)
	{
	  free_vector(xit,1,n);
	  free_vector(ptt,1,n);
	  free_vector(pt,1,n);
	  return;
	}
      for (j=1;j<=n;j++)
	{
	  ptt[j]=2.0*p[j]-pt[j];
	  /* CGA: HORRIBLE KLUDGE */
	  if ( ptt[j] < p_min[j] )
	    p[j] = p_min[j];
	  if ( ptt[j] > p_max[j] )
	    p[j] = p_max[j];
	  xit[j]=p[j]-pt[j];
	  pt[j]=p[j];
	}
      if ( debug_messages )
	{
	  printf( "TESTING WHETHER TO CHANGE DIRECTIONS\n" ); /* CGA */
	}
      fptt=(*func)(ptt);
      if (fptt < fp)
	{
	  t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
	  if (t < 0.0)
	    {
	      linmin(p,p_min,p_max,xit,n,fret,func);
	      for (j=1;j<=n;j++) 
		{
		  /* xi[j][ibig]=xit[j]; */
		  xi[j][ibig]=xi[j][n];
		  xi[j][n]=xit[j];
		}
	      if ( debug_messages )
		{
		  printf( "CHANGING DIRECTION %d\n", ibig ); /*CGA*/
		  for ( j = 1; j <= n; j++ ) /* CGA */
		    { /* CGA */
		      for ( i = 1; i <= n; i++ ) /* CGA */
			{ /* CGA */
			  if ( i == ibig ) /* CGA */
			    printf( " *" ); /* CGA */
			  printf( "%f ", xi[j][i] ); /* CGA */
			} /* CGA */
		      printf( "\n" ); /* CGA */
		    } /* CGA */
		}
	    }
	}
    }
}

#undef ITMAX
#undef SQR
