/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "stdlib.h"
#include "string.h"
#include "compute_expand_chunk_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "group.h"
#include "comm.h"
#include "force.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

enum{COMPUTE,FIX};

#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4

/* ---------------------------------------------------------------------- */

ComputeExpandChunkAtom::ComputeExpandChunkAtom(LAMMPS *lmp, int narg, char **arg) :
  Compute(lmp, narg, arg)
{
  if (narg < 5) error->all(FLERR,"Illegal compute expandchunk/atom command");

  MPI_Comm_rank(world,&me);
  
  // ID of compute chunk/atom

  int n = strlen(arg[3]) + 1;
  idchunk = new char[n];
  strcpy(idchunk,arg[3]);
  
  // parse remaining values until one isn't recognized

  which = new int[narg-4];
  argindex = new int[narg-4];
  ids = new char*[narg-4];
  value2index = new int[narg-4];
  nvalues = 0;

  peratom_flag = 1;
  size_peratom_cols = 0;
  
  for (int iarg = 4; iarg < narg; iarg++) {
    if (strncmp(arg[iarg],"c_",2) == 0 ||
        strncmp(arg[iarg],"f_",2) == 0) {
      if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
      else if (arg[iarg][0] == 'f') which[nvalues] = FIX;

      int n = strlen(arg[iarg]);
      char *suffix = new char[n];
      strcpy(suffix,&arg[iarg][2]);

      char *ptr = strchr(suffix,'[');
      if (ptr) {
        if (suffix[strlen(suffix)-1] != ']')
          error->all(FLERR,"Illegal compute expandchunk/atom command");
        argindex[nvalues] = atoi(ptr+1);
        *ptr = '\0';
      } else argindex[nvalues] = 0;

      n = strlen(suffix) + 1;
      ids[nvalues] = new char[n];
      strcpy(ids[nvalues],suffix);
      nvalues++;
      delete [] suffix;

    } else error->all(FLERR,"Illegal compute expandchunk/atom command");
  }
	
  if (nvalues == 1) size_peratom_cols = 0;
  else size_peratom_cols = nvalues;	

  init();
  
  nchunk = cchunk->setup_chunks();
  cchunk->compute_ichunk();
  int *ichunk = cchunk->ichunk;
  
  for (int i = 0; i < nvalues; i++) {
    if (which[i] == COMPUTE) {
      int icompute = modify->find_compute(ids[i]);
      if (icompute < 0)
        error->all(FLERR,"Compute ID for compute expandchunk/atom does not exist");
      if (modify->compute[icompute]->vector_flag) {
        if (argindex[i])
          error->all(FLERR,"Compute expandchunk/atom compute does not calculate a global array");
        if (nchunk != modify->compute[icompute]->size_vector)
          error->all(FLERR,"Compute expandchunk/atom compute vector does not have size nchunk");
      } else if (modify->compute[icompute]->array_flag) {
        if (argindex[i] == 0)
          error->all(FLERR,"Compute expandchunk/atom compute does not calculate a global vector");
        if (argindex[i] > modify->compute[icompute]->size_array_cols)
          error->all(FLERR,"Compute expandchunk/atom compute array is accessed out-of-range");
        if (nchunk != modify->compute[icompute]->size_array_rows)
          {
			char buffer [500];
			sprintf(buffer,"Compute expandchunk/atom compute array (%d) does not have size nchunk = %d",modify->compute[icompute]->size_array_rows,nchunk);
			error->all(FLERR,buffer);
			//error->all(FLERR,"Compute expandchunk/atom compute array does not have size nchunk");
		  }
      } else error->all(FLERR,"Compute expandchunk/atom compute does not calculate "
                        "global vector or array");
    } else if (which[i] == FIX) {
      int ifix = modify->find_fix(ids[i]);
      if (ifix < 0)
        error->all(FLERR,"Fix ID for compute expandchunk/atom does not exist");
      if (modify->fix[ifix]->vector_flag) {
        if (argindex[i])
          error->all(FLERR,"Compute expandchunk/atom fix does not calculate a global array");
        if (nchunk != modify->fix[ifix]->size_vector)
          error->all(FLERR,"Compute expandchunk/atom fix vector does not have size nchunk");
      } else if (modify->fix[ifix]->array_flag) {
        if (argindex[i] == 0)
          error->all(FLERR,"Compute expandchunk/atom fix does not calculate a global vector");
        if (argindex[i] > modify->fix[ifix]->size_array_cols)
          error->all(FLERR,"Compute expandchunk/atom fix array is accessed out-of-range");
        if (nchunk != modify->fix[ifix]->size_array_rows)
          error->all(FLERR,"Compute expandchunk/atom fix array does not have size nchunk");
      } else error->all(FLERR,"Compute expandchunk/atom fix does not calculate "
                        "global vector or array");
    }
  }
  
  // chunk-based data
  maxchunk = 0;
  allocate();

  nmax = 0;
  vector = NULL;
  array = NULL;
}

/* ---------------------------------------------------------------------- */

ComputeExpandChunkAtom::~ComputeExpandChunkAtom()
{
  delete [] idchunk;
  delete [] index;
  memory->destroy(vector);
  memory->destroy(array);
  
  delete [] which;
  delete [] argindex;
  for (int m = 0; m < nvalues; m++) delete [] ids[m];
  delete [] ids;
  delete [] value2index;
}

/* ---------------------------------------------------------------------- */

void ComputeExpandChunkAtom::init()
{
  int icompute = modify->find_compute(idchunk);
  if (icompute < 0)
    error->all(FLERR,"Chunk/atom compute does not exist for compute expandchunk/atom");
  cchunk = (ComputeChunkAtom *) modify->compute[icompute];
  if (strcmp(cchunk->style,"chunk/atom") != 0)
    error->all(FLERR,"Compute expandchunk/atom does not use chunk/atom compute");

  // set indices and check validity of all computes,fixes

  for (int m = 0; m < nvalues; m++) {
    if (which[m] == COMPUTE) {
      int icompute = modify->find_compute(ids[m]);
      if (icompute < 0)
        error->all(FLERR,"Compute ID for compute expandchunk/atom does not exist");
      value2index[m] = icompute;
    } else if (which[m] == FIX) {
      int ifix = modify->find_fix(ids[m]);
      if (ifix < 0)
        error->all(FLERR,"Fix ID for compute expandchunk/atom does not exist");
      value2index[m] = ifix;
    }
  }
}

/* ---------------------------------------------------------------------- */

void ComputeExpandChunkAtom::compute_peratom()
{
  invoked_peratom = update->ntimestep;
  
  // compute chunk/atom assigns atoms to chunk IDs
  // extract ichunk index vector from compute
  // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms

  nchunk = cchunk->setup_chunks();
  cchunk->compute_ichunk();
  int *ichunk = cchunk->ichunk;

  if (nchunk > maxchunk) allocate();
  
  if (atom->nlocal > nmax) {
    nmax = atom->nmax;
    if (nvalues == 1) {
      memory->destroy(vector);
      memory->create(vector,nmax,"expandchunk/atom:vector");
      vector_atom = vector;
    } else {
      memory->destroy(array);
      memory->create(array,nmax,nvalues,"expandchunk/atom:array");
      array_atom = array;
    }
  }

  // fill vector or array with per-atom values

  int *mask = atom->mask;
  int nlocal = atom->nlocal;

  if (nvalues == 1) {
	/*.............*/
	buf = vector;
	Compute *compute = modify->compute[value2index[0]];

    if (argindex[0] == 0) {
      if (!(compute->invoked_flag & INVOKED_VECTOR)) {
        compute->compute_vector();
        compute->invoked_flag |= INVOKED_VECTOR;
      }
      double *cvector = compute->vector;
      
      for (int i = 0; i < nlocal; i++) {
        if (mask[i] & groupbit) {
          buf[i] = cvector[(ichunk[i])-1];
        } else buf[i] = 0.00;
	  }

    } else {
      if (!(compute->invoked_flag & INVOKED_ARRAY)) {
        compute->compute_array();
        compute->invoked_flag |= INVOKED_ARRAY;
      }
      double **carray = compute->array;
      int icol = argindex[0]-1;
      
      for (int i = 0; i < nlocal; i++) {
        if (mask[i] & groupbit) {
          buf[i] = carray[(ichunk[i])-1][icol];
        } else buf[i] = 0.00;
	  }
    }
	
	/*.............*/
	
	/*buf = vector;
    for (int i = 0; i < nlocal; i++) {
      if (mask[i] & groupbit) {
        buf[i] = static_cast<int> (ichunk[i]);
      } else buf[i] = static_cast<int> (0);
	}*/
  } else {
    
    if (nmax) buf = &array[0][0];
    else buf = NULL;
    for (int n = 0; n < nvalues; n++) {
	  /*..........*/
	  Compute *compute = modify->compute[value2index[n]];
	  	
	  if (argindex[n] == 0) {
      if (!(compute->invoked_flag & INVOKED_VECTOR)) {
        compute->compute_vector();
        compute->invoked_flag |= INVOKED_VECTOR;
      }
      double *cvector = compute->vector;
      
      for (int i = 0; i < nlocal; i++) {
        if (mask[i] & groupbit) {
          buf[n+i*nvalues] = cvector[(ichunk[i])-1];
        } else buf[n+i*nvalues] = 0.00;
	  }

	  } else {
        if (!(compute->invoked_flag & INVOKED_ARRAY)) {
          compute->compute_array();
          compute->invoked_flag |= INVOKED_ARRAY;
        }
        double **carray = compute->array;
        int icol = argindex[n]-1;
      
        for (int i = 0; i < nlocal; i++) {
          if (mask[i] & groupbit) {
            buf[n+i*nvalues] = carray[(ichunk[i])-1][icol];
          } else buf[n+i*nvalues] = 0.00;
	    }
      }
		
	  /* ........ */	
	  /*for (int i = 0; i < nlocal; i++) {
        if (mask[i] & groupbit) {
          buf[n+i*nvalues] = static_cast<int> (ichunk[i]);
        } else buf[n+i*nvalues] = static_cast<int> (0);
	  }*/
	}
  }
  
}

void ComputeExpandChunkAtom::allocate()
{
  maxchunk = nchunk;
}

/* ----------------------------------------------------------------------
   memory usage of local atom-based array
------------------------------------------------------------------------- */

double ComputeExpandChunkAtom::memory_usage()
{
  double bytes = nmax*nvalues * sizeof(double);
  return bytes;
}