/*
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);
}
syntax highlighted by Code2HTML, v. 0.9.1