#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <sys/resource.h>
#define NMAX 8192

void get_cputime(double *laptime, double *splittime)
{
  struct timeval tval;
  struct timezone *dum=0;
  gettimeofday(&tval,dum);
  *laptime = tval.tv_sec + tval.tv_usec * 1.0e-6 - *splittime;
  *splittime = tval.tv_sec + tval.tv_usec * 1.0e-6 ;
}

#pragma goose func
double rsqrt(double r2)
{
  return 1.0/sqrt(r2);
}

int
main(int argc, char **argv)
{
  static double m[NMAX];
  static double x[NMAX][3];
  static double xi[NMAX][3];
  static double a[NMAX][3];
  static double pot[NMAX];
  double r2, mrinv, mr3inv, rinv, eps2, dx[3];
  int i, j, k, n, nstep, step;
  double lt, st, total_time = 0.0;

  if(argc == 3){
    n = atoi(argv[1]);
    nstep = atoi(argv[2]);
  }else if(argc == 2){
    n = atoi(argv[1]);
    nstep = 1;
  }else{
    n = 18;
    nstep = 1;
  }
  printf("n:%d  nstep:%d\n", n, nstep);

  for(step=0;step<nstep;step++){
    eps2 = (1.0/128.0)*(1.0/128.0);
    for(i=0;i<n;i++){
      x[i][0] = drand48()-0.5;
      x[i][1] = drand48()-0.5;
      x[i][2] = drand48()-0.5;
      m[i] = 1.0/((double)n);
    }

    get_cputime(&lt,&st);

#pragma omp parallel for private(j,dx,r2,rinv,mrinv,mr3inv,k)
#pragma goose parallel for precision ("double") loopcounter(i, j) \
        result(a[i][0..2], pot[i])
    for(i=0;i<n;i++) {
      for(k=0;k<3;k++) a[i][k] = 0.0;
      pot[i] = 0.0;
      for(j=0;j<n;j++) {
        for(k=0;k<3;k++) dx[k] = x[j][k] - x[i][k];
        r2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] + eps2;
        rinv = rsqrt(r2);
        mrinv = m[j]*rinv;
        mr3inv = mrinv*rinv*rinv;
        for(k=0;k<3;k++) a[i][k] += mr3inv * dx[k];
        pot[i] -= mrinv;
      }
    }
    get_cputime(&lt,&st);
    total_time += lt;

    for(i=0;i<n;i++) pot[i] += m[i]/sqrt(eps2);
  }

  for(i=0;i<n;i++){
    printf("i %2d    a %22.14e %22.14e %22.14e    p %22.14e\n",
           i, a[i][0], a[i][1], a[i][2], pot[i]);
  }
  printf("total time: %g\n",total_time);
  exit(0);
}