/*
    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.
*/

#include "header.h"

/* Set dir1 or dir2 to be a negative number if you only want to do a derivative
   in one direction. */
#define X 0
#define Y 1
#define Z 2

void Dabove( struct Grid *out, 
		struct Grid *self, 
		int sign,
		double forwardfilter[], 
		double reversefilter[],
		double fullfilter[],
		struct Grid *workgrid, 
		double *workspace1,
		double *workspace2, int dir1, int dir2 )
{ 
  int i,j;

  /* Return immediately if there is nothing to do */
  /* No excutables BEFORE this statement, please! */
  if (self==NULL) return;

    /*Make sure output on toplevel is zero if we are not
    doing a cumulative convolution*/
  if ( (sign==0) && (self->level == 0) )
      for (i=0;i<(out->dim.x*out->dim.y*out->dim.z);i++)
          out->dat[i]=0;
  
  
  if ( ( self->level > 1 )   && (self->child) )
    /* Take I from workgrid->parent onto workgrid. We also
       add the values of in onto the target workgrid, to ensure
       that we actually take I of in onto workgrid without
       changing in itself*/
      aconvP3dadd( workgrid->parent->dat+
                   (workgrid->org.x * self->parent->dim.y 
                    + workgrid->org.y) 
                   * self->parent->dim.z
                   + workgrid->org.z,
                   workgrid->dat,
                   self->dat,
                   tr, tf, ts,
                   workgrid->parent->dim.x, 
                   workgrid->parent->dim.y,
                   workgrid->parent->dim.z,
                   workgrid->dim.x, 
                   workgrid->dim.y,
                   workgrid->dim.z,
                   HFL_T, workspace1,workspace2 );
  else if ( ( self->level == 1 )   && (self->child) )
      /* Take I from in->parent onto workgrid */
      aconvP3dadd( self->parent->dat+
                   (workgrid->org.x * self->parent->dim.y 
                    + workgrid->org.y) 
                   * self->parent->dim.z
                   + workgrid->org.z,
                   workgrid->dat,
                   self->dat,
                   tr, tf, ts,
                   self->parent->dim.x, 
                   self->parent->dim.y,
                   self->parent->dim.z,
                   workgrid->dim.x, 
                   workgrid->dim.y,
                   workgrid->dim.z,
                   HFL_T, workspace1,workspace2 );

    /* We need to set the 'dir' part of the filter to the derivative 
       filter in order to computer the proper {x,y,z} derivative.  The simplest
       way to do this is to set all components to the O filter, and then
       reset the proper part to D.  This is wasteful, but it's only a few 
       flops.  -- mhe 04/12/02 */
  double Oscale=pow(2,self->level);
  for (j=0;j<=HFL_cD2;j++) 
    {
      reversefilter[0*(1+HFL_cD2)+j] = (O2l_r[j])/Oscale;
      forwardfilter[0*(1+HFL_cD2)+j] = (O2l_f[j]);
      fullfilter[0*(1+HFL_cD2)+j] = (O2l_s[j])/Oscale;
      
      reversefilter[1*(1+HFL_cD2)+j] = (O2l_r[j])/Oscale;
      forwardfilter[1*(1+HFL_cD2)+j] = (O2l_f[j]);
      fullfilter[1*(1+HFL_cD2)+j] = (O2l_s[j])/Oscale;
      
      reversefilter[2*(1+HFL_cD2)+j] = (O2l_r[j])/Oscale;
      forwardfilter[2*(1+HFL_cD2)+j] = (O2l_f[j]);
      fullfilter[2*(1+HFL_cD2)+j] = (O2l_s[j])/Oscale;
    }
  double Lscale=Oscale;
  for(j=0;j<=HFL_cD2;j++)
    {
      /* If dir1 == dir2, then we're really doing L. */
      if(dir1 == dir2)
	{
/*  	  switch(dir1) */
/*  	    { */
/*  	    case 0: */
/*  	      Lscale=self->scale.x; */
/*  	      break; */
/*  	    case 1: */
/*  	      Lscale=self->scale.y; */
/*  	      break; */
/*  	    case 2: */
/*  	      Lscale=self->scale.z; */
/*  	      break; */
/*  	    } */
	  reversefilter[dir1*(1+HFL_L2)+j]=(L2l_r[j])*Oscale;
	  forwardfilter[dir1*(1+HFL_L2)+j]=(L2l_f[j]);
	  fullfilter[dir1*(1+HFL_L2)+j]=(L2l_s[j])*Oscale;
	}
      if((dir1 >= 0) && (dir1 <= 2) && (dir1 != dir2))
	{
	  reversefilter[dir1*(1+HFL_cD2)+j] = (cD2l_r[j]);
	  forwardfilter[dir1*(1+HFL_cD2)+j] = (cD2l_f[j]);
	  fullfilter[dir1*(1+HFL_cD2)+j] = (cD2l_s[j]);
	}
      if((dir2 >= 0) && (dir2 <=2) && (dir1 != dir2))
	{
	  reversefilter[dir2*(1+HFL_cD2)+j] = (cD2l_r[j]);
	  forwardfilter[dir2*(1+HFL_cD2)+j] = (cD2l_f[j]);
	  fullfilter[dir2*(1+HFL_cD2)+j] = (cD2l_s[j]);
	}
    }
  
  if (workgrid->level > 1)
      /*D2level from workgrid->parent to out*/
    aconvP3d( workgrid->parent->dat+
	      +(workgrid->org.x * workgrid->parent->dim.y + 
		workgrid->org.y) * workgrid->parent->dim.z
	      + workgrid->org.z,
	      out->dat,
	      sign, reversefilter, forwardfilter, fullfilter,
	      workgrid->parent->dim.x,
	      workgrid->parent->dim.y,
	      workgrid->parent->dim.z,
	      workgrid->dim.x, 
	      workgrid->dim.y, 
	      workgrid->dim.z,
	      HFL_cD2, workspace1,workspace2 );
  else if (workgrid->level == 1)
    /*D2level directly from in->parent to out*/
    aconvP3d( self->parent->dat+
	      +(out->org.x * self->parent->dim.y + 
		out->org.y) * self->parent->dim.z
	      + out->org.z,
	      out->dat,
	      sign, reversefilter, forwardfilter, fullfilter,
	      self->parent->dim.x,
	      self->parent->dim.y,
	      self->parent->dim.z,
	      out->dim.x, 
	      out->dim.y, 
	      out->dim.z,
	      HFL_cD2, workspace1, workspace2 );

  /* Recursive call to Dabove for children. */
  if (self->child)
    Dabove(out->child,self->child,sign,
	      forwardfilter, reversefilter, fullfilter,
	      workgrid->child, workspace1, workspace2, dir1, dir2 );
  
  /* Recursive call to Dabove for siblings.*/
  if (self->sibling)
    Dabove(out->sibling,self->sibling,sign,
	      forwardfilter,reversefilter, fullfilter,
	      workgrid->sibling,workspace1, workspace2, dir1, dir2 );
  
}


void Dsame(struct Grid *out, struct Grid *self, int sign,
	      double forwardfilter[], 
	      double reversefilter[],
	      double fullfilter[],
	      double *workspace1,
	      double *workspace2, int dir1, int dir2)
{
  int j;

    /* We need to set the 'dir' part of the filter to the derivative 
       filter in order to computer the proper {x,y,z} derivative.  The simplest
       way to do this is to set all components to the O filter, and then
       reset the proper part to D.  This is wasteful, but it's only a few 
       flops.  -- mhe 04/12/02 */
  double Oscale=pow(2,self->level);
  for (j=0;j<=HFL_cD1;j++) 
    {
      reversefilter[0*(1+HFL_cD1)+j] = (O1l_r[j])/Oscale;
      forwardfilter[0*(1+HFL_cD1)+j] = (O1l_f[j]);
      fullfilter[0*(1+HFL_cD1)+j] = (O1l_s[j])/Oscale;
      
      reversefilter[1*(1+HFL_cD1)+j] = (O1l_r[j])/Oscale;
      forwardfilter[1*(1+HFL_cD1)+j] = (O1l_f[j]);
      fullfilter[1*(1+HFL_cD1)+j] = (O1l_s[j])/Oscale;
      
      reversefilter[2*(1+HFL_cD1)+j] = (O1l_r[j])/Oscale;
      forwardfilter[2*(1+HFL_cD1)+j] = (O1l_f[j]);
      fullfilter[2*(1+HFL_cD1)+j] = (O1l_s[j])/Oscale;
    }
  double Lscale=Oscale;
  for(j=0;j<=HFL_cD1;j++)
    {
      /* If dir1 == dir2, then we're really doing L. */
      if(dir1 == dir2)
	{
/*  	  switch(dir1) */
/*  	    { */
/*  	    case 0: */
/*  	      Lscale=self->scale.x; */
/*  	      break; */
/*  	    case 1: */
/*  	      Lscale=self->scale.y; */
/*  	      break; */
/*  	    case 2: */
/*  	      Lscale=self->scale.z; */
/*  	      break; */
/*  	    } */
	  reversefilter[dir1*(1+HFL_L1)+j]=(L1l_r[j])*Lscale;
	  forwardfilter[dir1*(1+HFL_L1)+j]=(L1l_f[j]);
	  fullfilter[dir1*(1+HFL_L1)+j]=(L1l_s[j])*Lscale;
	}
      if((dir1 >= 0) && (dir1 <=2) && (dir1 != dir2))
	{
	  reversefilter[dir1*(1+HFL_cD1)+j] = (cD1l_r[j]);
	  forwardfilter[dir1*(1+HFL_cD1)+j] = (cD1l_f[j]);
	  fullfilter[dir1*(1+HFL_cD1)+j] = (cD1l_s[j]);
	}
      if((dir2 >= 0) && (dir2 <=2) && (dir1 != dir2))
	{
	  reversefilter[dir2*(1+HFL_cD1)+j] = (cD1l_r[j]);
	  forwardfilter[dir2*(1+HFL_cD1)+j] = (cD1l_f[j]);
	  fullfilter[dir2*(1+HFL_cD1)+j] = (cD1l_s[j]);
	}
    }

  /*Take Dsame on self. Need to make it possible to distinguish
   between periodic boundary conditions and non periodic b.c.*/  
  /* IPD: Now the grids know about the boundary conditions */

  /* Do convolution */
  conv3dnew(self->dat, out->dat, sign, self->ifperiodic,
            reversefilter,forwardfilter, fullfilter,
            self->dim.x, self->dim.y, self->dim.z,
            HFL_cD1, workspace1, workspace2);


  if (self->sibling) 
    Dsame(out->sibling,self->sibling,sign,
	     forwardfilter,reversefilter,fullfilter,
	     workspace1, workspace2, dir1, dir2);
	       
  if (self->child) 
    Dsame(out->child,self->child,sign,
	     forwardfilter,reversefilter, fullfilter,
	     workspace1, workspace2, dir1, dir2);
}

void 
Dbelow( struct Grid *out, struct Grid *self, int sign,
		double forwardfilter[],
                double reversefilter[],
                double fullfilter[],
		double *workspace1, double *workspace2, int dir1, int dir2)
{
  int i,j,k;
  
  /*Same process on siblings*/
  if (self->sibling) 
    Dbelow(out->sibling,self->sibling,sign,
	      forwardfilter,reversefilter,fullfilter,
	      workspace1,workspace2,dir1,dir2);
    
  /*Same process on children*/
  if (self->child) 
    Dbelow(out->child,self->child,sign,
	      forwardfilter,reversefilter,fullfilter,
	      workspace1,workspace2,dir1,dir2); 


  /*Have to zero out every gridpoint on bottom level,
    will not be necessary in final version of Obelow, need
    to change that eventually*/
  if ( (sign == 0) && (self->child==NULL) )
     for (i=0;i<(out->dim.x*out->dim.y*out->dim.z);i++)
      out->dat[i]=0;

  if ( (self->level > 0)  ) 
    {
      /* We need to set the 'dir' part of the filter to the derivative 
	 filter in order to computer the proper {x,y,z} derivative.  The 
	 simplest way to do this is to set all components to the O filter, and
	 then reset the proper part to D.  This is wasteful, but it's only a 
	 few flops.  -- mhe 04/12/02 */
      double Oscale=pow(2,self->level);
      for (j=0;j<=HFL_cD2;j++) 
	{
	  reversefilter[0*(1+HFL_cD2)+j] = (O2l_r[j])/Oscale;
	  forwardfilter[0*(1+HFL_cD2)+j] = (O2l_f[j]);
	  fullfilter[0*(1+HFL_cD2)+j] = (O2l_s[j])/Oscale;
	  
          reversefilter[1*(1+HFL_cD2)+j] = (O2l_r[j])/Oscale;
	  forwardfilter[1*(1+HFL_cD2)+j] = (O2l_f[j]);
	  fullfilter[1*(1+HFL_cD2)+j] = (O2l_s[j])/Oscale;

	  reversefilter[2*(1+HFL_cD2)+j] = (O2l_r[j])/Oscale;
	  forwardfilter[2*(1+HFL_cD2)+j] = (O2l_f[j]);
	  fullfilter[2*(1+HFL_cD2)+j] = (O2l_s[j])/Oscale;
	}
      double Lscale=Oscale;
      for(j=0;j<=HFL_cD2;j++)
	{
	  /* If dir1 == dir2, then we're really doing L. */
	  if(dir1 == dir2)
  	    { 
/*  	      switch(dir1) */
/*  		{ */
/*  		case 0: */
/*  		  Lscale=self->scale.x; */
/*  		  break; */
/*  		case 1: */
/*  		  Lscale=self->scale.y; */
/*  		  break; */
/*  		case 2: */
/*  		  Lscale=self->scale.z; */
/*  		  break; */
/*  		} */
	      reversefilter[dir1*(1+HFL_L2)+j]=(L2l_r[j])*Lscale;
	      forwardfilter[dir1*(1+HFL_L2)+j]=(L2l_f[j]);
	      fullfilter[dir1*(1+HFL_L2)+j]=(L2l_s[j])*Lscale;
	    }
	  if((dir1 >= 0) && (dir1 <=2) && (dir1 != dir2))
	    {
	      reversefilter[dir1*(1+HFL_cD2)+j] = (cD2l_r[j]);
	      forwardfilter[dir1*(1+HFL_cD2)+j] = (cD2l_f[j]);
	      fullfilter[dir1*(1+HFL_cD2)+j] = (cD2l_s[j]);
	    }
	  if((dir2 >= 0) && (dir2 <=2) && (dir1 != dir2))
	    {
	      reversefilter[dir2*(1+HFL_cD2)+j] = (cD2l_r[j]);
	      forwardfilter[dir2*(1+HFL_cD2)+j] = (cD2l_f[j]);
	      fullfilter[dir2*(1+HFL_cD2)+j] = (cD2l_s[j]);
	    }
	}

      if (sign == 0)
	if (self->sibling == NULL)
	  /*Have to zero out every gridpoint on out->parent*/
	  for (i=0;
	       i<(out->parent->dim.x*out->parent->dim.y*out->parent->dim.z);
	       i++)
	    out->parent->dat[i]=0;
      
      /*  D2lconvolution from self to out->parent */
      aPconv3d( out->parent->dat+
		( (self->org.x) * out->parent->dim.y 
		  + (self->org.y)) * out->parent->dim.z
		+ (self->org.z),
		self->dat,
		1, reversefilter, forwardfilter, fullfilter,
		out->parent->dim.x, 
		out->parent->dim.y, 
		out->parent->dim.z, 
		self->dim.x,
		self->dim.y,
		self->dim.z,
		HFL_cD2, workspace1, workspace2 );

      /*This is a hack which could be improved by doing the addition inside
	of the aPconv3d function call below, but it only costs about 2% in 
	performance, so I will not modify the function to make 
	that possible for now..*/
      for (i=0; i<out->dim.x; i+=2)
	for (j=0; j<out->dim.y; j+=2)
	  for (k=0; k<out->dim.z; k+=2)
	    *( out->parent->dat + 
	       ( (out->org.x+i/2) * out->parent->dim.y + 
		 (out->org.y+j/2) ) * out->parent->dim.z +
	       (out->org.z+k/2) )
	      +=
	      *( out->dat+(i*out->dim.y+j)*out->dim.z+k );

      /*  Idag convolution from out to out->parent
	  This part of the function ensures that we
	  get the contribution to O from ALL levels
	  below, not just the level right below
	  Need to make condition a little more rigorous than
	  this*/
      if (out->child) 
	aPconv3d( out->parent->dat+
		  ( (out->org.x) * out->parent->dim.y 
		    + (out->org.y)) * out->parent->dim.z
		  + (out->org.z),
		  out->dat,
		  1, tr, tf, ts,
		  out->parent->dim.x, 
		  out->parent->dim.y, 
		  out->parent->dim.z, 
		  out->dim.x,
		  out->dim.y,
		  out->dim.z,
		  HFL_T, workspace1, workspace2 );
    }
}

/* Multi-level first derivative operator, for x, y, or z directions.
   workgrid:  a grid onto which we place the results of
              some intermediate calculations
   workspace: additional workspace which is used inside 
              the convolutions	     */
void
cD( struct Grid *out, struct Grid *in,
   double D1lf[], double D1lr[], double D1lfull[],
   double D2lf[], double D2lr[], double D2lfull[],
   struct Grid *workgrid,double *workspace1, double *workspace2, 
   int dir1, int dir2 )
{
  /* Sanity check for dir1 and dir2. */
  if((dir1 > 2) || (dir2 > 2))
    {
      printf("Error in D:  directions must be 0 (x), 1 (y), or 2 (z).\n");
      printf("Exiting.\n");
      exit(1);
    }
  
  if ( (HFL_T>HFL_O2) || (HFL_O1>HFL_O2) || (HFL_T>HFL_O1) ) {
    printf("Warning: our code assumes HFL_O2 > HFL_O1 > HFL_T\n");
    printf("while you have entered the lengths %d %d %d\n\n",
	   HFL_O2,HFL_O1,HFL_T);
    printf("This creates problems with the size of workspace\n");
    printf("so you must change the parameters determining the workspace.\n");
    exit(1);
  }

  Dbelow(out,in,0,D2lf,D2lr,D2lfull,workspace1,workspace2,dir1,dir2);

  Dabove(out,in,1,D2lf,D2lr,D2lfull,workgrid,workspace1,workspace2,dir1,dir2);

  Dsame(out,in,1,D1lf,D1lr,D1lfull,workspace1,workspace2,dir1,dir2);

  zeroghostrec(out);
}