/*
    DFT++ is a density functional package developed by the research group
    of Professor Tomas Arias

    Copyright 1996-2003 Sohrab Ismail-Beigi

    This file is part of DFT++.

    DFT++ is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    DFT++ is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DFT++; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    Please see the file CREDITS for a list of authors.

    For academic users, we request that publications using results obtained with
    this software reference

    "New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
    and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).

    and, if using the wavelet basis, further reference

    "Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
    T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).

    and 

    "Robust ab initio calculation of condensed matter: transparent convergence through
    semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
    Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).

    For your convenience, preprints of the above articles may be obtained from
    http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/

/* $Id: transtest.cpp,v 1.1.2.3 2003/05/29 18:54:36 ivan Exp $ */
/*
 * $Log: transtest.cpp,v $
 * Revision 1.1.2.3  2003/05/29 18:54:36  ivan
 * Update the PRL reference
 *
 * Revision 1.1.2.2  2003/05/23 01:49:17  sahak
 * added the new header
 *
 * Revision 1.1.2.1  2002/09/22 23:36:03  ivan
 * Incorporating development by Prof. Arias:
 *   GGA - both for PW and WL branch
 *   optest/transtest/gradtest for WL (note that for PW are trivial)
 *        updated and incorporated to this code, added option WLtests
 *        in makefile to compile them
 *
 * Revision 1.1  1999/09/14 19:40:35  torkel
 * Initial revision
 *
 * Revision 1.7  1999/08/06 19:25:39  torkel
 * *** empty log message ***
 *
 * Revision 1.6  1999/07/04 20:49:33  torkel
 * Full O implemented.
 *
 * Revision 1.5  1999/06/11 12:56:42  torkel
 * *** empty log message ***
 *
 * Revision 1.2  1999/04/25 21:00:31  muchomas
 * Longer run set for timing
 *
 * Revision 1.1  1999/04/21 02:32:32  muchomas
 * Initial revision
 *
 * Revision 1.2  1999/04/20 15:59:37  muchomas
 * Fully documented code
 *
 * Revision 1.1  1999/04/19 21:09:48  muchomas
 * Initial revision
 * *
 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>
#include "WL_header.h"

#define HFL 3
#define TIMELOOP 50

#define PX 3
#define PY 3
#define PZ 3

int Px,Py,Pz;

/* Computes a multinomial function at a point in space (OUTPUT)
   USES multinomial exponents from GLOBAL variables Px, Py, Pz!!!
   r: real-space point
   return value: value of function
*/
double multi(struct Dvec r) {
  return( pow(r.x,Px)+pow(r.y,Py)+pow(r.z,Pz) );
}

main()
{
  /* Space for grid specification structure */
  struct Gridspec spec[NgridMx];

  /* Pointer to special grids (will be allocated below) */
  struct Grid *ones,*onesc,*nononesc,*zeros,*ghosts,*nonghosts;

  /* Size needed for workspace */
  struct Ivec szmx={0,0,0};

  /* Pointer to "grid"s */
  struct Grid *grid,*grid1,*grid2;

  /* Scale and origin structure for top level */
  struct Dvec scale={1.,1.,1.};
  struct Dvec org={-6.,-6.,-12.};

  /* Work space variables */
  double *work,*work_extra;
  int lwork;

  /* Timer variables */
  struct timeval time1,time2;
  double rtime,mflop;

  /* Misc variables */
  FILE *fp;
  double tmp;
  int i,j,k;

  printf("\n-----\n");
  /* Read grid structure from file */
  if ( (fp=fopen("gridspec","r")) == NULL )
    {
      fprintf(stderr,"Can't open file \"gridspec\" for input.\n");
      fprintf(stderr,"Exiting.\n");
      exit(1);
    }
  readgridspec(fp,spec);
  close(fp);
  /* Set top level grid spacings to something non-trivial */
  spec[0].scale.x=2.; spec[0].scale.y=3.; spec[0].scale.z=4.;
  printf("Seting (dx,dy,dz)=(%lf,%lf,%lf)...\n\n",
	 spec[0].scale.x,spec[0].scale.y,spec[0].scale.z);

  /* Make and allocate special grids */
  zeros=mkgridnew(spec,NULL);
  ones=mkgridnew(spec,NULL);
  onesc=mkgridnew(spec,NULL);
  nononesc=mkgridnew(spec,NULL);
  ghosts=mkgridnew(spec,NULL);
  nonghosts=mkgridnew(spec,NULL);

  /* Get needed work array size and memory */
  getgridsizemx(ones,&szmx);
  lwork=(szmx.x+3*HFL)*(szmx.y+3*HFL)*(szmx.z+3*HFL);
  if ( (work=(double *) malloc(lwork*sizeof(double))) == NULL )
    {
      fprintf(stderr,"malloc failed for %d doubles in main (1)\n",lwork);
      fprintf(stderr,"Exiting.\n");
      exit(1);
    }
  if ( (work_extra=(double *) malloc(lwork*sizeof(double))) == NULL )
    {
      fprintf(stderr,"malloc failed for %d doubles in main (1)\n",lwork);
      fprintf(stderr,"Exiting.\n");
      exit(1);
    }
  
  for(i=0;i<szmx.x;i++)
    for(j=0;j<szmx.y;j++)
      for(k=0;k<szmx.z;k++)
    {
      work[i*szmx.y*szmx.z+j*szmx.z+k]=rand(); 
      work_extra[i*szmx.y*szmx.z+j*szmx.z+k]=rand(); 
    }

  /* Fill in special grids */
  evalgrid(zeros,org,zero);
  evalgrid(ones,org,one);
  evalgrid(onesc,org,one);
  J(onesc,work,work_extra,lwork);
  gridsub(nononesc,ones,onesc);
  gridghost(ghosts);
  gridsub(nonghosts,ones,ghosts);

  /* Make and allocate "grid"s */
  grid=mkgridnew(spec,NULL);
  grid1=mkgridnew(spec,NULL);
  grid2=mkgridnew(spec,NULL);

  /* Time I */
  evalgrid(grid,org,rnd);

  printf("Timing %d iterations... ",TIMELOOP);
  gettimeofday(&time1,NULL);
  for (i=0; i<TIMELOOP; i++)
    I(grid,work,work_extra,lwork);
  gettimeofday(&time2,NULL);
  rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
  for (i=0; spec[i].level!=-1; i++)
    {
      if (spec[i].level==0)
	mflop=0.;
      else
	mflop+=
	  (spec[i].dim.x+2*HFL)*(spec[i].dim.y+2*HFL)*(spec[i].dim.z+2*HFL)*
	  (1.+0.5+0.25)*(3+7)*(TIMELOOP)/1e6;
    }
  printf("   I:  %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
	  mflop,rtime,mflop/rtime);

  evalgrid(grid,org,one); 

  /* Time J */
  evalgrid(grid,org,rnd);

  printf("Timing %d iterations... ",TIMELOOP);
  gettimeofday(&time1,NULL);
  for (i=0; i<TIMELOOP; i++)
    J(grid,work,work_extra,lwork);
  gettimeofday(&time2,NULL);
  rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
  for (i=0; spec[i].level!=-1; i++)
    {
      if (spec[i].level==0)
	mflop=0.;
      else
	mflop+=
	  (spec[i].dim.x+2*HFL)*(spec[i].dim.y+2*HFL)*(spec[i].dim.z+2*HFL)*
	  (1.+0.5+0.25)*(3+7)*(TIMELOOP)/1e6;
    }
  printf("   J:  %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
	  mflop,rtime,mflop/rtime);

  /* Time Idag */
  evalgrid(grid,org,rnd);

  printf("Timing %d iterations... ",TIMELOOP);
  gettimeofday(&time1,NULL);
  for (i=0; i<TIMELOOP; i++)
    Idag(grid,work,work_extra,lwork);
  gettimeofday(&time2,NULL);
  rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
  for (i=0; spec[i].level!=-1; i++)
    {
      if (spec[i].level==0)
	mflop=0.;
      else
	mflop+=
	  (spec[i].dim.x+2*HFL)*(spec[i].dim.y+2*HFL)*(spec[i].dim.z+2*HFL)*
	  (1.+0.5+0.25)*(7+3.5)*(TIMELOOP)/1e6;
    }
  printf("Idag:  %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
	  mflop,rtime,mflop/rtime);

  /* Time Jdag */
  evalgrid(grid,org,rnd);

  printf("Timing %d iterations... ",TIMELOOP);
  gettimeofday(&time1,NULL);
  for (i=0; i<TIMELOOP; i++)
    Jdag(grid,work,work_extra,lwork);
  gettimeofday(&time2,NULL);
  rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
  for (i=0; spec[i].level!=-1; i++)
    {
      if (spec[i].level==0)
	mflop=0.;
      else
	mflop+=
	  (spec[i].dim.x+2*HFL)*(spec[i].dim.y+2*HFL)*(spec[i].dim.z+2*HFL)*
	  (1.+0.5+0.25)*(7+3.5)*(TIMELOOP)/1e6;
    }
  printf("Jdag:  %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n\n",
	  mflop,rtime,mflop/rtime);

  printf("Starting this: \n");

  /* Test linearity of J */
  evalgrid(grid1,org,rnd);
  evalgrid(grid2,org,rnd);
  gridadd(grid,grid1,grid2);

  J(grid,work,work_extra,lwork);
  J(grid1,work,work_extra,lwork);
  J(grid2,work,work_extra,lwork);

  gridsub(grid,grid,grid1);
  gridsub(grid,grid,grid2);
  printf("J is linear: %le %le\n",sqrt(griddot(grid,grid)), 
	 sqrt(griddot(grid2,grid2)) );
  
  /* Test linearity of I (being sure to zero ghost points of input!) */
  evalgrid(grid1,org,rnd); gridptmult(grid1,nonghosts,grid1);
  evalgrid(grid2,org,rnd); gridptmult(grid2,nonghosts,grid2);
  gridadd(grid,grid1,grid2);

  I(grid,work,work_extra,lwork);
  I(grid1,work,work_extra,lwork);
  I(grid2,work,work_extra,lwork);

  gridsub(grid,grid,grid1);
  gridsub(grid,grid,grid2);
  printf("I is linear: %le %le\n\n",sqrt(griddot(grid,grid)),
	 sqrt(griddot(grid2,grid2)) );

  /* Test J */
  evalgrid(grid,org,rnd);
  gridptmult(grid1,ones,grid);

  J(grid,work,work_extra,lwork);
  printf("J is zero on ghosts: %le\n",sqrt(gridckghosts(grid,0)));

  I(grid,work,work_extra,lwork);
  gridsub(grid2,grid,grid1); 
  gridptmult(grid2,grid2,grid2);
  printf("I is inverse of J: %le\n\n",sqrt(gridsum(grid2,0)));

  /* Test I */
  evalgrid(grid,org,one); 
  gridptmult(grid,nonghosts,grid);
  gridptmult(grid1,ones,grid);

  I(grid,work,work_extra,lwork);
  printf("I matches on ghosts: %le\n",sqrt(gridckghosts(grid,1)));

  J(grid,work,work_extra,lwork);

  gridsub(grid2,grid,grid1); gridptmult(grid2,grid2,grid2); 
  printf("J is inverse of I: %le\n\n",sqrt(gridsum(grid2,0)));

  /* Test Jdag */
  evalgrid(grid,org,rnd); gridptmult(grid,nonghosts,grid);
  gridptmult(grid1,ones,grid);

  Jdag(grid,work,work_extra,lwork);
  printf("Jdag matches on ghosts: %le\n",sqrt(gridckghosts(grid,1)));

  Idag(grid,work,work_extra,lwork);
  gridsub(grid2,grid,grid1); 
  gridptmult(grid2,grid2,grid2);
  printf("Idag is inverse of Jdag: %le\n\n",sqrt(gridsum(grid2,0)));

  /* grid1 in coeff space, grid2 in real space */
  evalgrid(grid1,org,rnd); gridptmult(grid1,nonghosts,grid);
  evalgrid(grid2,org,rnd); 

  /* tmp=1'J(2) */
  J(grid2,work,work_extra,lwork); tmp=griddot(grid1,grid2); I(grid2,work,work_extra,lwork);

  /* tmp-=Jdag(1)'2 */
  Jdag(grid1,work,work_extra,lwork); tmp-=griddot(grid1,grid2);
  printf("Jdag is dag of J: %le %le  \n\n",tmp, griddot(grid1,grid2) );


  /* Test Idag */
  evalgrid(grid,org,rnd); 
  gridptmult(grid1,ones,grid);

  Idag(grid,work,work_extra,lwork);
  printf("Idag zero on ghosts: %le\n",sqrt(gridckghosts(grid,0)));

  Jdag(grid,work,work_extra,lwork);
  gridsub(grid2,grid,grid1); gridptmult(grid2,grid2,grid2);
  printf("Jdag is inverse of Idag: %le\n\n",sqrt(gridsum(grid2,0)));

  /* grid1 in real space, grid2 in coef space */
  evalgrid(grid1,org,rnd); 
  evalgrid(grid2,org,rnd); gridptmult(grid2,nonghosts,grid);

  /* tmp=1'I(2) */
  I(grid2,work,work_extra,lwork); tmp=griddot(grid1,grid2); J(grid2,work,work_extra,lwork);

  /* tmp-=Jdag(1)'2 */
  Idag(grid1,work,work_extra,lwork); tmp-=griddot(grid1,grid2);
  printf("Idag is dag of I: %le\n\n",tmp);

  tmp=0.;
  for (Px=0; Px<=PX; Px++)
    for (Py=0; Py<=PY; Py++)
      for (Pz=0; Pz<=PZ; Pz++)
	{
	  evalgrid(grid,org,multi);
	  J(grid,work,work_extra,lwork);
	  gridptmult(grid,grid,nononesc); gridptmult(grid,grid,grid);
	  tmp+=gridsum(grid,0);
	}
  printf("Multinomial tests summed up to (%d,%d,%d) [typically 1e-8]: %le\n",PX,PY,PZ,sqrt(tmp));

	{
	  Px=0; Py=0; Pz=4;
	  evalgrid(grid,org,multi);
	  J(grid,work,work_extra,lwork);
	  gridptmult(grid,grid,nononesc); gridptmult(grid,grid,grid);
	  printf("Test of (%d,%d,%d): %18.12lf\n\n",Px,Py,Pz,sqrt(gridsum(grid,0)));
	}
	  

printf("-----\n\n");
}