/* 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 #include #include #include #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; idim.x; i++) for (j=0; jdim.y; j++) for (k=0; kdim.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 %.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 %.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 %.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 %.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 %.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 %.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 %.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 %.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); } }