/*
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 <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys/time.h>
#include "WL_header.h"
#define TIMELOOP 10
int Px,Py,Pz;
double evalmflopOsame(struct Gridspec spec[])
{
int i;
double mflop=0;
for (i=0; spec[i].level!=-1; i++) {
mflop+=
spec[i].dim.x*spec[i].dim.y*spec[i].dim.z*
(6*(2*HFL_O1+1)
-HFL_O1*(1/((double) spec[i].dim.x)+1
/((double) spec[i].dim.y)+1/((double) spec[i].dim.z)))
*(TIMELOOP)/1e6;
}
return(mflop);
}
double evalmflopOabove(struct Gridspec spec[])
{
int i;
double mflop=0;
double tmp;
for (i=0; spec[i].level!=-1; i++) {
if (spec[i].level==0)
mflop=0.;
else {
mflop+=
(2*HFL_O2+1)
* (3./2. * spec[i].dim.x *spec[i].dim.y *spec[i].dim.z)
* (
(1./4.* (1+3*HFL_O2/((double) spec[i].dim.x)))*
(1+2*HFL_O2/((double) spec[i].dim.y))*
(1+2*HFL_O2/((double) spec[i].dim.z))
+1./2. * (1+3*HFL_O2/((double)spec[i].dim.y))*(1+2*HFL_O2/((double) spec[i].dim.z))
+1 * (1+3*HFL_O2/((double) spec[i].dim.z))
)
*(TIMELOOP)/1e6;
if (spec[i+1].level > spec[i].level)
mflop+=
(2*HFL_T+1)
* (3./2. * spec[i].dim.x *spec[i].dim.y *spec[i].dim.z)
* (
(1./4.* (1+3*HFL_T/((double) spec[i].dim.x)))*
(1+2*HFL_T/((double) spec[i].dim.y))*
(1+2*HFL_T/((double) spec[i].dim.z))
+1./2. * (1+3*HFL_T/((double)spec[i].dim.y))*(1+2*HFL_T/((double) spec[i].dim.z))
+1 * (1+3*HFL_T/((double) spec[i].dim.z))
)
*(TIMELOOP)/1e6;
}
}
return(mflop);
}
double evalmflopObelow(struct Gridspec spec[])
{
int i;
double mflop=0;
double rh;
for (i=0; spec[i].level!=-1; i++) {
if (spec[i].level==0)
mflop=0.;
else {
rh=(5.*HFL_O2*HFL_O2+2.*HFL_O2-2.)/(2.*HFL_O2+1.);
mflop+=(
3./2.* spec[i].dim.x *spec[i].dim.y *spec[i].dim.z * (2*HFL_O2)
* (
1 * (1 + 2*rh/(3.* spec[i].dim.x))
+1./2. * (1 + 2*rh/(3. *spec[i].dim.y)) * (1 + 2*rh/((double) spec[i].dim.x))
+1./4. * (1 + 2*rh/(3. *spec[i].dim.z)) * (1 + 2*rh/((double) spec[i].dim.y))* (1 + 2*rh/((double) spec[i].dim.x)) ) )
*(TIMELOOP)/1e6;
if (spec[i+1].level > spec[i].level)
rh=(5.*HFL_T*HFL_T+2.*HFL_T-2.)/(2.*HFL_T+1.);
mflop+=(
3./2.* spec[i].dim.x *spec[i].dim.y *spec[i].dim.z * (2*HFL_T)
* (
1 * (1 + 2*rh/(3.* spec[i].dim.x))
+1./2. * (1 + 2*rh/(3. *spec[i].dim.y)) * (1 + 2*rh/((double) spec[i].dim.x))
+1./4. * (1 + 2*rh/(3. *spec[i].dim.z)) * (1 + 2*rh/((double) spec[i].dim.y))* (1 + 2*rh/((double) spec[i].dim.x)) ) )
*(TIMELOOP)/1e6;
}
}
return(mflop);
}
double evalmflopO(struct Gridspec spec[])
{
return( evalmflopObelow(spec)
+ evalmflopOsame(spec)
+ evalmflopOabove(spec) );
}
double evalmflopLsame(struct Gridspec spec[])
{
int i;
double mflop=0;
for (i=0; spec[i].level!=-1; i++) {
mflop+=3*
spec[i].dim.x*spec[i].dim.y*spec[i].dim.z*
(6*(2*HFL_O1+1)
-HFL_O1*(1/((double) spec[i].dim.x)+1/((double) spec[i].dim.y)+1/((double) spec[i].dim.z)))
*(TIMELOOP)/1e6;
}
return(mflop);
}
double evalmflopLabove(struct Gridspec spec[])
{
int i;
double mflop=0;
double tmp;
for (i=0; spec[i].level!=-1; i++) {
if (spec[i].level==0)
mflop=0.;
else {
mflop+=3*
(2*HFL_O2+1)
* (3./2. * spec[i].dim.x *spec[i].dim.y *spec[i].dim.z)
* (
(1./4.* (1+3*HFL_O2/((double) spec[i].dim.x)))*
(1+2*HFL_O2/((double) spec[i].dim.y))*
(1+2*HFL_O2/((double) spec[i].dim.z))
+1./2. * (1+3*HFL_O2/((double)spec[i].dim.y))
*(1+2*HFL_O2/((double) spec[i].dim.z))
+1 * (1+3*HFL_O2/((double) spec[i].dim.z))
)
*(TIMELOOP)/1e6;
if (spec[i+1].level > spec[i].level)
mflop+=
(2*HFL_T+1)
* (3./2. * spec[i].dim.x *spec[i].dim.y *spec[i].dim.z)
* (
(1./4.* (1+3*HFL_T/((double) spec[i].dim.x)))*
(1+2*HFL_T/((double) spec[i].dim.y))*
(1+2*HFL_T/((double) spec[i].dim.z))
+1./2. * (1+3*HFL_T/((double)spec[i].dim.y))
*(1+2*HFL_T/((double) spec[i].dim.z))
+1 * (1+3*HFL_T/((double) spec[i].dim.z))
)
*(TIMELOOP)/1e6;
}
}
return(mflop);
}
double evalmflopLbelow(struct Gridspec spec[])
{
int i;
double mflop=0;
double rh;
for (i=0; spec[i].level!=-1; i++) {
if (spec[i].level==0)
mflop=0.;
else {
rh=(5.*HFL_O2*HFL_O2+2.*HFL_O2-2.)/(2.*HFL_O2+1.);
mflop+=3*(
3./2.* spec[i].dim.x *spec[i].dim.y *spec[i].dim.z * (2*HFL_O2)
* (
1 * (1 + 2*rh/(3.* spec[i].dim.x))
+1./2. * (1 + 2*rh/(3. *spec[i].dim.y))
* (1 + 2*rh/((double) spec[i].dim.x))
+1./4. * (1 + 2*rh/(3. *spec[i].dim.z))
* (1 + 2*rh/((double) spec[i].dim.y))
* (1 + 2*rh/((double) spec[i].dim.x)) ) )
*(TIMELOOP)/1e6;
if (spec[i+1].level > spec[i].level)
rh=(5.*HFL_T*HFL_T+2.*HFL_T-2.)/(2.*HFL_T+1.);
mflop+=(
3./2.* spec[i].dim.x *spec[i].dim.y
*spec[i].dim.z * (2*HFL_T)
* (
1 * (1 + 2*rh/(3.* spec[i].dim.x))
+1./2. * (1 + 2*rh/(3. *spec[i].dim.y))
* (1 + 2*rh/((double) spec[i].dim.x))
+1./4. * (1 + 2*rh/(3. *spec[i].dim.z))
* (1 + 2*rh/((double) spec[i].dim.y))
* (1 + 2*rh/((double) spec[i].dim.x)) ) )
*(TIMELOOP)/1e6;
}
}
return(mflop);
}
double evalmflopL(struct Gridspec spec[])
{
return( evalmflopLbelow(spec)
+ evalmflopLsame(spec)
+ evalmflopLabove(spec) );
}
/**************************************************************/
/*Functions used for verifying that all operators are correct */
double compgrids(struct Grid *in1,struct Grid *in2)
{
int i,j,k;
int num;
double tmp=0.;
/*Don't expect things to come out correctly on top level
anyways*/
if (in1->parent)
for (i=0; i<in1->dim.x; i++)
for (j=0; j<in1->dim.y; j++)
for (k=0; k<in1->dim.z; k++)
if ( (i%2 == 1) || (j%2 == 1) || (k%2 == 1) ) {
num=(i*in1->dim.y+j)*in1->dim.z+k;
if ( fabs(in1->dat[num] / (1 + fabs(in2->dat[num])) ) > tmp)
tmp=fabs(in1->dat[num] / (1 + fabs(in2->dat[num])) );
}
if (in1->sibling)
tmp+=compgrids(in1->sibling,in2->sibling);
if (in1->child)
tmp+=compgrids(in1->child,in2->child);
return(tmp);
}
double multiL(struct Dvec r)
{
double tot=0;
if (Px>1)
tot=Px*(Px-1)*pow(r.x,Px-2)*pow(r.y,Py)*pow(r.z,Pz);
if (Py>1)
tot+=Py*(Py-1)*pow(r.x,Px)*pow(r.y,Py-2)*pow(r.z,Pz);
if (Pz>1)
tot+=Pz*(Pz-1)*pow(r.x,Px)*pow(r.y,Py)*pow(r.z,Pz-2);
return(tot);
}
double multi(struct Dvec r)
{
return( pow(r.x,Px)*pow(r.y,Py)*pow(r.z,Pz) );
}
/*
Takes account of the scaling factors from the level of the grid
and physical dimensions. Used in testing
O(J(P))=det(topdims)*2^(-sd)*P(p)
*/
void multfunc( struct Grid *self)
{
int i;
for (i=0;i<(self->dim.x*self->dim.y*self->dim.z);i++)
self->dat[i] /= self->scale.x*self->scale.y*self->scale.z;
if (self->child)
multfunc(self->child);
if (self->sibling)
multfunc(self->sibling);
}
main()
{
struct Ivec NULLIVEC={0,0,0};
/* 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;
/* Pointer to helpgrids (will be allocated below) */
struct Grid *workgrid;
/* Size needed for workspace */
struct Ivec szmx={0,0,0};
/* Pointer to "grid"s */
struct Grid *grid,*grid1,*grid2,*grid3,*grid4,*grid5,*grid6,*grid7,*grid8;
/* Scale and origin structure for top level */
struct Dvec org={-6.,-6.,-12.};
double O1lr[3*(HFL_O1+1)],O1lf[3*(HFL_O1+1)],O1lfull[3*(HFL_O1+1)];
double O2lr[3*(HFL_O2+1)],O2lf[3*(HFL_O2+1)],O2lfull[3*(HFL_O2+1)];
/* Work space variables */
double *workspace1,*workspace2;
int lwork;
int num;
/* Timer variables */
struct timeval time1,time2;
double rtime,mflop;
/* Misc variables */
FILE *fp;
double tmp,tmp2;
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);
workgrid=mkgridnew(spec,NULL);
/* Print grid structure */
printf("\nGrid structure:\n"); printgridstruct(stdout,ones); printf("\n");
/* Get needed work array size and memory */
getgridsizemx(ones,&szmx);
lwork=(szmx.x+3*HFL_O2)*(szmx.y+3*HFL_O2)*(szmx.z+3*HFL_O2);
if ( (workspace1=(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 ( (workspace2=(double *) malloc(lwork*sizeof(double))) == NULL ) {
fprintf(stderr,"malloc failed for %d doubles in main (1)\n",lwork);
fprintf(stderr,"Exiting.\n");
exit(1);
}
/* Fill in special grids */
evalgrid(zeros,org,zero);
evalgrid(ones,org,one);
evalgrid(onesc,org,one); J(onesc,workspace1,workspace2,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);
grid3=mkgridnew(spec,NULL);
grid4=mkgridnew(spec,NULL);
grid5=mkgridnew(spec,NULL);
grid6=mkgridnew(spec,NULL);
/* Timing Oabove */
printf("Timing %d Oabove iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
Oabovenew(grid3,grid1,0,O2lf,O2lr,O2lfull,grid2,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" Oabovenew: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
evalmflopOabove(spec),rtime,evalmflopOabove(spec)/rtime);
/*Timing Osamenew*/
printf("Timing %d Osame iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
Osamenew(grid1,grid2,0,O1lf,O1lr,O1lfull,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" Osamenew: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
evalmflopOsame(spec),rtime,evalmflopOsame(spec)/rtime);
/*Timing Obelownew*/
printf("Timing %d Obelownew iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
Obelownew(grid1,grid2,0,O2lf,O2lr,O2lfull,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" Obelownew: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n\n\n",
evalmflopObelow(spec),rtime,evalmflopObelow(spec)/rtime);
/*Timing O*/
printf("Timing %d O iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
O(grid1,grid2,O1lf,O1lr,O1lfull,O2lf,O2lr,O2lfull,
grid3,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" O: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n\n\n",
evalmflopO(spec),rtime,evalmflopO(spec)/rtime);
/* Timing Labove */
printf("Timing %d Labove iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
Labovenew(grid3,grid1,0,O2lf,O2lr,O2lfull,grid2,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" Labovenew: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
evalmflopLabove(spec),rtime,evalmflopLabove(spec)/rtime);
/*Timing Lsamenew*/
printf("Timing %d Lsame iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
Lsamenew(grid1,grid2,0,O1lf,O1lr,O1lfull,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" Lsamenew: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n",
evalmflopLsame(spec),rtime,evalmflopLsame(spec)/rtime);
/*Timing Lbelownew*/
printf("Timing %d Lbelownew iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
Lbelownew(grid1,grid2,0,O2lf,O2lr,O2lfull,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" Lbelownew: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n\n\n",
evalmflopLbelow(spec),rtime,evalmflopLbelow(spec)/rtime);
/*Timing L*/
printf("Timing %d L iterations... \n",TIMELOOP);
gettimeofday(&time1,NULL);
for (i=0; i<TIMELOOP; i++)
L(grid1,grid2,O1lf,O1lr,O1lfull,O2lf,O2lr,O2lfull,
grid3,workspace1,workspace2);
gettimeofday(&time2,NULL);
rtime=1.*(time2.tv_sec-time1.tv_sec)+(time2.tv_usec-time1.tv_usec)/1.e6;
printf(" L: %.0lf MFLOP in %.3f sec ==> %.0lf MFLOPS.\n\n\n",
evalmflopL(spec),rtime,evalmflopL(spec)/rtime);
/*Check if Osame is Hermitian*/
evalgrid(grid4,org,rnd);
zeroghostrec(grid4);
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
Osamenew(grid1,grid2,0,O1lf,O1lr,O1lfull,workspace1,workspace2);
Osamenew(grid3,grid4,0,O1lf,O1lr,O1lfull,workspace1,workspace2);
tmp=griddot(grid1,grid4);
tmp2=griddot(grid2,grid3);
printf("Osame is dag of Osame to %d digits\n\n",
-(int) log10(fabs((1e-100+tmp-tmp2)/tmp)) );
/*Check if Oabovenew is dagger of Obelownew*/
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
evalgrid(grid4,org,rnd);
zeroghostrec(grid4);
evalgrid(grid5,org,rnd);
Oabovenew(grid1,grid2,0,O2lf,O2lr,O2lfull,grid3,workspace1,workspace2);
Obelownew(grid5,grid4,0,O2lf,O2lr,O2lfull,workspace1,workspace2);
tmp=griddot(grid1,grid4);
tmp2=griddot(grid2,grid5);
printf("Oabovenew is dag of Obelownew to %d digits. \n\n",
-(int) log10(1.e-100+fabs((tmp-tmp2)/tmp)) );
/*Check if Oabovenew is dagger of Obelow*/
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
evalgrid(grid4,org,rnd);
zeroghostrec(grid4);
copygrid(grid5,grid4);
Oabovenew(grid1,grid2,0,O2lf,O2lr,O2lfull,grid3,workspace1,workspace2);
Obelowsimple(grid5,workspace1,workspace2,O2lf,O2lr,O2lfull);
tmp=griddot(grid1,grid4);
tmp2=griddot(grid2,grid5);
printf("Oabovenew is dag of Obelow to %d digits.\n\n",
-(int) log10(fabs((1e-100+tmp-tmp2)/tmp)));
/*Compare Obelow to Obelownew*/
evalgrid(grid3,org,zero);
evalgrid(grid1,org,rnd);
zeroghostrec(grid1);
copygrid(grid2,grid1);
Obelowsimple(grid1,workspace1,workspace2,O2lf,O2lr,O2lfull);
Obelownew(grid3,grid2,0,O2lf,O2lr,O2lfull,workspace1,workspace2);
gridsub(grid4,grid1,grid3);
printf("Obelow is equal to Obelownew to %d digits.\n\n",
(int) -log10(fabs( (1e-100+gridsum(grid4,0))/gridsum(grid1,0))) );
/*Check if O is Hermitian*/
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
evalgrid(grid5,org,one);
zeroghostrec(grid5);
O(grid1,grid2,O1lf,O1lr,O1lfull,O2lf,O2lr,O2lfull,
grid3,workspace1,workspace2);
O(grid4,grid5,O1lf,O1lr,O1lfull,O2lf,O2lr,O2lfull,
grid6,workspace1,workspace2);
tmp=griddot(grid1,grid5);
tmp2=griddot(grid2,grid4);
printf("O is Hermitian to %d digits\n\n",
-(int) log10(fabs((1e-100+tmp-tmp2)/tmp)) );
/*Check if Lsame is Hermitian*/
evalgrid(grid4,org,rnd);
zeroghostrec(grid4);
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
Lsamenew(grid1,grid2,0,O1lf,O1lr,O1lfull,workspace1,workspace2);
Lsamenew(grid3,grid4,0,O1lf,O1lr,O1lfull,workspace1,workspace2);
tmp=griddot(grid1,grid4);
tmp2=griddot(grid2,grid3);
printf("Lsame is Hermitian to %d digits\n\n",
-(int) log10(fabs((1e-100+tmp-tmp2)/tmp)) );
/*Check if Labovenew is dagger of Lbelownew*/
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
evalgrid(grid4,org,rnd);
zeroghostrec(grid4);
evalgrid(grid5,org,rnd);
Labovenew(grid1,grid2,0,O2lf,O2lr,O2lfull,grid3,workspace1,workspace2);
Lbelownew(grid5,grid4,0,O2lf,O2lr,O2lfull,workspace1,workspace2);
tmp=griddot(grid1,grid4);
tmp2=griddot(grid2,grid5);
printf("Labove is dag of Lbelow to %d digits.\n\n",
-(int) log10(1e-300+fabs((tmp-tmp2)/tmp)));
/*Check if L is Hermitian*/
evalgrid(grid2,org,rnd);
zeroghostrec(grid2);
evalgrid(grid5,org,rnd);
zeroghostrec(grid5);
L(grid1,grid2,O1lf,O1lr,O1lfull,O2lf,O2lr,O2lfull,
grid3,workspace1,workspace2);
L(grid4,grid5,O1lf,O1lr,O1lfull,O2lf,O2lr,O2lfull,
grid6,workspace1,workspace2);
tmp=griddot(grid1,grid5);
tmp2=griddot(grid2,grid4);
printf("L is Hermitian to %d digits\n\n",
-(int) log10(fabs((1e-100+tmp-tmp2)/tmp)) );
/*Testing that O(J(P))=det(topdims)*2^(-sd)*P(p) */
tmp=0.;
printf("Multinomial tests for O, zero to following # of digits up to (3,3,3) [typically ~12 digits]\n");
for (Px=0; Px<=3; Px++)
for (Py=0; Py<=3; Py++)
for (Pz=0; Pz<=3; Pz++) {
evalgrid(grid1,org,multi);
evalgrid(grid3,org,multi);
zeroghostrec(grid3);
J(grid1,workspace1,workspace2,lwork);
O(grid2,grid1,O1lf,O1lr,O1lfull,
O2lf,O2lr,O2lfull,grid4,workspace1,workspace2);
multfunc(grid2);
gridsub(grid1,grid2,grid3);
tmp=compgrids(grid1,grid2);
printf(" %d %d %d %d ", Px,Py,Pz, (int) -log10(tmp) );
printf("\n");
}
{
Px=0; Py=0; Pz=4;
evalgrid(grid1,org,multi);
evalgrid(grid3,org,multi);
zeroghostrec(grid3);
J(grid1,workspace1,workspace2,lwork);
O(grid2,grid1,O1lf,O1lr,O1lfull,
O2lf,O2lr,O2lfull,grid4,workspace1,workspace2);
multfunc(grid2);
gridsub(grid1,grid2,grid3);
tmp=compgrids(grid1,grid2);
printf("Test of (%d,%d,%d): %18.12lf\n\n",Px,Py,Pz,tmp);
}
/*Testing that L(J(P))=det(topdims)*2^(-sd)*P(p) */
tmp=0.;
printf("Multinomial tests for L, zero to following # of digits up to (3,3,3) [typically ~7-8 digits]\n");
for (Px=0; Px<=3; Px++)
for (Py=0; Py<=3; Py++)
for (Pz=0; Pz<=3; Pz++)
{
evalgrid(grid1,org,multi);
evalgrid(grid3,org,multiL);
zeroghostrec(grid3);
J(grid1,workspace1,workspace2,lwork);
L(grid2,grid1,O1lf,O1lr,O1lfull,
O2lf,O2lr,O2lfull,grid4,workspace1,workspace2);
multfunc(grid2);
gridsub(grid1,grid2,grid3);
tmp=compgrids(grid1,grid2)+1e-100;
printf(" %d %d %d %d \n ",Px,Py,Pz, (int) -log10(tmp) );
}
printf("\n\n");
{
Px=0; Py=0; Pz=4;
evalgrid(grid1,org,multi);
evalgrid(grid3,org,multiL);
zeroghostrec(grid3);
J(grid1,workspace1,workspace2,lwork);
L(grid2,grid1,O1lf,O1lr,O1lfull,
O2lf,O2lr,O2lfull,grid4,workspace1,workspace2);
multfunc(grid2);
gridsub(grid1,grid2,grid3);
tmp=compgrids(grid1,grid2);
printf("Test of (%d,%d,%d): %18.12lf\n\n",Px,Py,Pz,tmp);
}
}
syntax highlighted by Code2HTML, v. 0.9.1