/*File: lennard-jones_ga.c
  Purpose: to test Lennard-Jones problem using Tao/GA vectors using multiple processors.
  Date created: August 16, 2002
*/
     

/**************************************************************

Author: Limin Zhang, Ph.D.
        Mathematics Department
        Columbia Basin College
        Pasco, WA 99301
        Limin.Zhang@cbc2.org

Mentor: Jarek Naplocha, Ph.D.
        Environmental Molecular Science Laboratory
        Pacific Northwest National Laboratory
        Richland, WA 99352

Date: 4/22/2002

Purpose:
      to design and implement an application for  between TAO and
      global arrays.
**************************************************************/

/* Include MPI header files */
#include <mpi.h> 

/* Program usage: mpirun -machinefile m -np 1 lennard-jones_ga [-help] [all TAO options] */
/* Program usage: mpirun -machinefile m -np 4 lennard-jones_ga [-help] [all TAO options] */

/*  Include tao header files so we can use TAO  */
#include "tao.h" 
#include "taoapp_ga.h" 

/* Include ga header files so that we can use GA */
#include "ga.h" 
#include "taovec_ga.h" 
#include "macdecls.h" 


/*some standard header files */
#include <stdio.h> 
#include <math.h> 

#define NDIM 2 
#define NATOMS 16 

static char help[] = "This example demonstrates use of the TAO package to \n\
solve an unconstrained minimization problem on multiple processors.  We \n\
minimize the extended Lennard-Jones function.";

/*T 
   Concepts: TAO - Solving an unconstrained minimization problem
   Routines: TaoInitialize(); TaoFinalize(); TaoSetFromOptions();
   Routines: TaoGAApplicationCreate(); TaoGASetInitialSolutionVec();
   Routines: TaoCreate(); TaoGASetObjectiveAndGradientRoutine();
   Routines: TaoSetGAHessian(); TaoSetGAInitialVector(); 
   Routines: TaoSolveGAApplication(); TaoDestroy(); TaoGAAppDestroy();
   Routines: TaoGetTerminationReason(); TaoView();
   Processors: 1 
T*/


/* 
   User-defined application context - contains data needed by the 
   application-provided call-back routines that evaluate the function,
   gradient, and hessian.
*/
typedef struct
{
  int n;			/* dimension n = NDIM*NATOMS */
  int ndim;
  int natoms;
  int memHandle;
}
AppCtx;


/* -------------- User-defined routines ---------- */

int FormFunctionGradient (TAO_GA_APPLICATION gaapp, GAVec ga_X, double *f, GAVec ga_G, void *ptr);
int InitializeVariables(GAVec ga_X, AppCtx *user); 


int main (int argc, char **argv)
{
  double startTime;
  int info;			/* used to check for functions returning nonzeros */
  GAVec ga_x;        		/* solution vector */
  TAO_SOLVER tao;		/* TAO_SOLVER solver context */
  TAO_GA_APPLICATION taoapp;	/* TAO application context */
  TaoTerminateReason reason;
  AppCtx user;			/* user-defined application context */

  /*initialize GA and MPI */
  int heap = 4000, stack = 4000;
  MPI_Init (&argc, &argv);	/* initialize MPI */
  GA_Initialize ();		/* initialize GA */
  if (!MA_init(MT_F_DBL, stack, heap))
    GA_Error("MA_init failed", stack+heap);

  /* Initialize TAO */
  TaoInitialize (&argc, &argv, (char *) 0, help);

  startTime = MPI_Wtime();

  /* Initialize problem parameters */
  user.natoms = NATOMS;
  user.ndim = NDIM;
  user.n = user.natoms*user.ndim;

  /* Create working space */
  if (MA_push_stack(C_DBL, 2*user.n, "Vector buffers", &user.memHandle) == MA_FALSE)
    ga_error("MAIN::ma_alloc_get failed",2*user.n);

  /* Allocate Global Array vector for the solution */
  int dims[2];
  dims[0] = user.n;
  ga_x = NGA_Create (C_DBL, 1, dims, "GA_X", NULL);
  if (!ga_x) GA_Error ("lennard-jones.main::NGA_Create ga_x", ga_x);

  /* The TAO code begins here */
  /* Create TAO solver with desired solution method */
  info = TaoCreate (MPI_COMM_WORLD, "tao_cg_fr", &tao); CHKERRQ(info);
  info = TaoGAApplicationCreate (MPI_COMM_WORLD, &taoapp); CHKERRQ(info);

  /* Set initial vector */
  info = InitializeVariables(ga_x, &user); CHKERRQ(info);
  info = TaoGAAppSetInitialSolutionVec(taoapp, ga_x); CHKERRQ(info);

  /* Set routines for function, gradient */
  info = TaoGAAppSetObjectiveAndGradientRoutine (taoapp, FormFunctionGradient, (void *) &user); 
  CHKERRQ(info);

  /* Check for TAO command line options */
  info = TaoSetFromOptions (tao); CHKERRQ(info);

  /* SOLVE THE APPLICATION */
  info = TaoSolveGAApplication (taoapp, tao); CHKERRQ(info);

  /*  To View TAO solver information use */
  info = TaoView(tao); CHKERRQ(info);

  /* Get termination information */
  info = TaoGetTerminationReason (tao, &reason); CHKERRQ(info);

  if (reason <= 0)
    printf("Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");

  printf("TIME TAKEN = %lf\n", MPI_Wtime()-startTime);

  /*output the solutions */ 
  printf ("The solution is :\n");
  GA_Print (ga_x);

  /* Free TAO data structures */
  info = TaoDestroy (tao); CHKERRQ(info);
  info = TaoGAAppDestroy (taoapp); CHKERRQ(info);

  /* Free GA data structures */
  GA_Destroy (ga_x);
  if (!MA_pop_stack(user.memHandle))
    ga_error("Main::MA_pop_stack failed",0);

  /* Finalize TAO, GA, and MPI */
  TaoFinalize ();
  GA_Terminate ();
  MPI_Finalize ();

  return 0;
}

/* -------------------------------------------------------------------- */
/*  
    FormFunctionGradient - Evaluates the function, f(X), and gradient, G(X). 

    Input Parameters:
.   tao  - the TAO_SOLVER context
.   ga_X    - input vector
.   ptr  - optional user-defined context, as set by TaoSetFunctionGradient()
    
    Output Parameters:
.   ga_G - vector containing the newly evaluated gradient
.   f - function value

    Note:
    Some optimization methods ask for the function and the gradient evaluation
    at the same time.  Evaluating both at once may be more efficient that
    evaluating each separately. 
*/


int FormFunctionGradient (TAO_GA_APPLICATION gaapp, GAVec ga_X, double *f, GAVec ga_G, void *ptr)
{
  int lo, hi;			//the global coordinates
  AppCtx *user = (AppCtx *) ptr;
  int i,j;
  double *g, *x;
  double xx,yy,zz,temp,rij;
  

  MA_get_pointer(user->memHandle, &x);
  g = x + user->ndim*user->natoms;
    
  lo=0;
  hi=user->n-1; /* range of array indices */

  NGA_Get(ga_X, &lo, &hi, x, &hi);

  *f = 0;
  for (i=0; i < user->n; i++)
    g[i] = 0.0;

  if (user->ndim == 2) {
    for (j=1; j < user->natoms; j++) {
      for (i=0; i<j; i++) {
	xx = x[2*j] - x[2*i];
	yy = x[2*j+1] - x[2*i+1];
	rij = xx*xx + yy*yy;
	temp = 1.0/rij/rij/rij;
	*f += temp*(temp-2.0);
	temp *= 12.0*(temp-1.0)/rij;
	g[2*j] -= xx*temp;
	g[2*j+1] -= yy*temp;
	g[2*i] += xx*temp;
	g[2*i+1] += yy*temp;
      }
    }
  } else if (user->ndim == 3) {
    for (j=1; j < user->natoms; j++) {
      for (i=0; i < j; i++) {
	xx = x[3*j] - x[3*i];
	yy = x[3*j+1] - x[3*i+1];
	zz = x[3*j+2] - x[3*i+2];
	rij = xx*xx + yy*yy + zz*zz;
	temp = 1.0/rij/rij/rij;
	*f += temp*(temp-2.0);
	temp *= 12.0*(temp-1.0)/rij;
	g[3*j] -= xx*temp;
	g[3*j+1] -= yy*temp;
	g[3*j+2] -= zz*temp;
	g[3*i] += xx*temp;
	g[3*i+1] += yy*temp;
	g[3*i+2] += zz*temp;
      }
    }
  }
      
  NGA_Put(ga_G, &lo, &hi, g, &hi);
  return 0;
}

InitializeVariables">int InitializeVariables(GAVec ga_X, AppCtx *user) 
{
  double *x;
  double xx, yy, zz;
  int isqrtn, icrtn, left, i, j, k, il, jl, ctr;
  int lo, hi;


  lo = 0; hi = user->n - 1; /* range of array indices */
  MA_get_pointer(user->memHandle, &x);

  NGA_Get (ga_X, &lo, &hi, x, &hi);

  if (user->ndim == 2) {
    isqrtn = (int) sqrt( (double) user->natoms);
    left = user->natoms - isqrtn * isqrtn;
    xx = 0.0;
    yy = 0.0;
    for (j=0; j<=isqrtn + left/isqrtn; j++) {
      for (i=0; i < TaoMin(isqrtn, user->natoms - j*isqrtn); i++) {
	ctr = j*isqrtn + i;
	x[2*ctr] = xx;
	x[2*ctr+1] = yy;
	xx += 1.0;
      }
      yy += 1.0;
      xx = 0.0;
    }	     
  }
  else if (user->ndim == 3) {
    icrtn = (int) pow((user->natoms + 0.5),1.0/3.0);
    left = user->natoms - icrtn * icrtn * icrtn;
    xx = yy = zz = 0.0;
    for (k=0; k <= icrtn + left; k++) {
      jl = TaoMin(icrtn, (user->natoms - k*icrtn*icrtn)/icrtn+1);
      for (j=0; j<jl; j++) {
	il = TaoMin(icrtn, user->natoms - k*icrtn*icrtn - j*icrtn);
	for (i=0; i<il; i++) {
	  ctr = k*icrtn*icrtn + j*icrtn + i;
	  x[3*ctr] = xx;
	  x[3*ctr+1] = yy;
	  x[3*ctr+2] = zz;
	  xx += 1.0;
	}
	yy += 1.0;
	xx = 0.0;
      }
      zz += 1.0;
      yy = 0.0;
    }
  }

  NGA_Put(ga_X, &lo, &hi, x, &hi);
  return 0;
}