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