    /*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;
    }