/* Libraries needed by various routines */
#include <stdio.h>
#include <stdlib.h>
/* Include standard C++ complex library */
/* And, make short-hand "Complex" for complex<double> */
#include <complex.h>
typedef complex<double> Complex;
/* Make C functions in nrutil.h accessible to C++ */
extern "C"
{
void fournp480(double data[], unsigned long nn[], int ndim, int isign);
#include "nrutil.h"
}
/* Global definition of pi */
const double pi=4.*atan(1.);
/* Include your C3tensor routines between here ... */
/* ... and here! */
/* Simple three-vector classes */
class i3vec {
public:
int x,y,z;
// Constructor
i3vec() {};
// Copy constructor
// Copy assignment
// Destructor
// Initializer
i3vec(int i,int j, int k) {
x=i; y=j; z=k;
}
};
class d3vec {
public:
double x,y,z;
// Constructor
d3vec() {};
// Copy constructor
// Copy assignment
// Destructor
// Initializer
d3vec(double i,double j, double k) {
x=i; y=j; z=k;
}
};
/* The pvector class */
class pvector {
struct i3vec N;
struct d3vec d;
public:
Complex ***dat;
// Constructor
pvector(i3vec iN, struct d3vec id) {
int i,j,k;
N.x=iN.x; d.x=id.x;
N.y=iN.y; d.y=id.y;
N.z=iN.z; d.z=id.z;
dat=C3tensor(0,N.x-1,0,N.y-1,0,N.z-1);
}
// Copy constructor
pvector(const pvector& in) {
int i,j,k;
N.x=in.N.x; d.x=in.d.x;
N.y=in.N.y; d.y=in.d.y;
N.z=in.N.z; d.z=in.d.z;
dat=C3tensor(0,N.x-1,0,N.y-1,0,N.z-1);
for (i=0; i<N.x; i++)
for (j=0; j<N.y; j++)
for (k=0; k<N.z; k++)
dat[i][j][k]=in.dat[i][j][k];
}
// Copy assignment: cleanup and copy
pvector& operator=(const pvector& in) {
int i,j,k;
for (i=0; i<N.x; i++)
for (j=0; j<N.y; j++)
for (k=0; k<N.z; k++)
dat[i][j][k]=in.dat[i][j][k];
}
// Destructor
~pvector() {
free_C3tensor(dat,0,N.x-1,0,N.y-1,0,N.z-1);
}
// Print out as 3 dimensional array
print3d(char filename[]) {
FILE* ostream;
if ( (ostream=fopen(filename,"w"))==NULL ) {
printf("Unable to open file '%s'!\n",filename);
exit(1);
}
int i,j,k;
for (i=0; i<N.x; i++) {
fprintf(ostream,"%4d ----------------------------------------------n",i);
for (j=0; j<N.y; j++) {
for (k=0; k<N.z; k++)
fprintf(ostream,"(%9f,%9f) ",
real(dat[i][j][k]),imag(dat[i][j][k]));
fprintf(ostream,"\n");
}
}
fprintf(ostream,"\n");
}
// Print out as radial data
printr(char filename[]) {
FILE* ostream;
if ( (ostream=fopen(filename,"w"))==NULL ) {
printf("Unable to open file '%s'!\n",filename);
exit(1);
}
// Loop over entire array
for (int i=0; i<N.x; i++)
for (int j=0; j<N.y; j++)
for (int k=0; k<N.z; k++) {
// Wrap box around so end of array is "near" origin
// (box dimensions are in an i3vec called "N")
int ii = ( i<N.x/2 ? i : i-N.x);
int jj = ( j<N.y/2 ? j : j-N.y);
int kk = ( k<N.z/2 ? k : k-N.z);
// Print distance to origin and data
fprintf(ostream,"%e %e %e\n",
sqrt( pow(ii*d.x,2)+pow(jj*d.y,2)+pow(kk*d.z,2) ),
dat[i][j][k]);
}
}
// vector Hermitian inner product
friend Complex operator^(const pvector& in1,const pvector& in2) {
Complex out;
// Initialize sum
out=Complex(0.,0.);
// Loop over all components, doing sum for dagger operation
// (box dimensions are in both i3vec's, "in1.N" or "in2.N")
for (int i=0; i<in2.N.x; i++)
for (int j=0; j<in2.N.y; j++)
for (int k=0; k<in2.N.z; k++)
out+=conj(in1.dat[i][j][k])*in2.dat[i][j][k];
return out;
}
// vector addition
friend pvector operator+(const pvector& in1,const pvector& in2) {
// Create a pvector called "out" which is a copy of pvector "in1"
pvector out(in1);
// Loop over all components, adding in2 onto out
// (which already contains in1)
for (int i=0; i<in2.N.x; i++)
for (int j=0; j<in2.N.y; j++)
for (int k=0; k<in2.N.z; k++)
out.dat[i][j][k]+=in2.dat[i][j][k];
return out;
}
// Initialize according to a function
eval(Complex (*func)()) {
for (int i=0; i<N.x; i++)
for (int j=0; j<N.y; j++)
for (int k=0; k<N.z; k++)
dat[i][j][k]=func();
}
// I (forward) transform
/* friend pvector I(const pvector& in) { */
/* pvector out(in); */
/* long unsigned int dim[4]; */
/* dim[1]=in.N.x; */
/* dim[2]=in.N.y; */
/* dim[3]=in.N.z; */
/* fournp480((double *)&out.dat[0][0][0]-1, dim, 3, 1); */
/* return out; */
/* } */
/* // L */
/* friend pvector L(const pvector& in) { */
/* // Get nornalization constants for G-vectors */
/* d3vec g=d3vec( */
/* 2.*pi/(in.N.x*in.d.x), */
/* 2.*pi/(in.N.y*in.d.y), */
/* 2.*pi/(in.N.z*in.d.z) */
/* ); */
/* d3vec g2=d3vec( */
/* g.x*g.x, */
/* g.y*g.y, */
/* g.z*g.z */
/* ); */
/* // Creat space for output and copy input into it */
/* pvector out(in); */
/* // Calculate volume of unit cell */
/* double V=(in.N.x*in.d.x)*(in.N.y*in.d.y)*(in.N.z*in.d.z); */
/* // Loop over entire box */
/* for (int i=0; i<in.N.x; i++) */
/* for (int j=0; j<in.N.y; j++) */
/* for (int k=0; k<in.N.z; k++) { */
/* // Treat aliasing properly!! */
/* int ii = ( i<in.N.x/2 ? i : i-in.N.x); */
/* int jj = ( j<in.N.y/2 ? j : j-in.N.y); */
/* int kk = ( k<in.N.z/2 ? k : k-in.N.z); */
/* out.dat[i][j][k]*=-V*(g2.x*(ii*ii)+g2.y*(jj*jj)+g2.z*(kk*kk)); */
/* } */
/* return out; */
/* } */
};
// Function to return random, complex numbers
Complex rnd() {
return Complex(
((double)rand())/RAND_MAX-0.5,
((double)rand())/RAND_MAX-0.5);
}
main() {
int n; /* Dimensions of FFT box */
double d=1.; /* Grid spacing in real space */
// Scan problem size and store in N
printf("n: "); scanf("%d",&n);
printf("Testing on %dx%dx%d grid ...\n",n,n,n);
// Initialize
pvector v1(i3vec(n,n,n),d3vec(d,d,d));
pvector v2(i3vec(n,n,n),d3vec(d,d,d));
pvector v3(i3vec(n,n,n),d3vec(d,d,d));
// Fill in with random data
v1.eval(rnd);
v2.eval(rnd);
v3.eval(rnd);
// Test that random filling is working
printf("Ok to be real, but better not be all zeros!\n");
printf(" |v1|^2=(%f,%f)\n",v1^v1);
printf(" |v2|^2=(%f,%f)\n",v2^v2);
printf(" |v3|^2=(%f,%f)\n\n",v3^v3);
// "Random" complex number for testing
Complex c1=Complex(2.,3.);
printf("These SHOULD be ALL zeros!!!\n");
// Test (v1^v2)=conj(v2^v1))
printf(" Hermitian inner product is Hermitian: (%15.12f,%15.12f)\n",
(v1^v2)-conj(v2^v1));
/* // Test sum(v1+v2)=sum(v1)+sum(v2) */
/* printf(" sum() is a linear operation: (%15.12f,%15.12f)\n", */
/* sum(v1+v2)-sum(v1)-sum(v2)); */
/* // Test pvector subtraction */
/* printf(" subtraction test: (%15.12f,%15.12f)\n", */
/* sum(v1-v2)-(sum(v1)-sum(v2))); */
/* // Test pvector products */
/* printf(" scalar-pvector product test: (%15.12f,%15.12f)\n", */
/* sum(c1*v1)-c1*sum(v1)); */
/* printf(" pvector-scalar product test: (%15.12f,%15.12f)\n", */
/* sum(v1*c1)-c1*sum(v1)); */
/* // Test algebra of inner product */
/* printf(" Hermitian inner product is linear: (%15.12f,%15.12f)\n", */
/* (v1^(v2+c1*v3))-((v1^v2)+c1*(v1^v3)) */
/* ); */
/* printf(" Hermitian inner product is co-linear: (%15.12f,%15.12f)\n\n", */
/* ((v2+c1*v3)^v1)-((v2^v1)+conj(c1)*(v3^v1)) */
/* ); */
/* // Simple test of I */
/* v3=I(v1+c1*v2)-(I(v1)+c1*I(v2)); */
/* printf(" I operator is linear: (%15.12f,%15.12f)\n", */
/* v3^v3); */
/* // Tests of J and I */
/* v3=J(I(v1))-v1; */
/* printf(" J is inverse of I: (%15.12f,%15.12f)\n", */
/* v3^v3); */
/* v3=I(J(v1))-v1; */
/* printf(" I is inverse of J: (%15.12f,%15.12f)\n", */
/* v3^v3); */
/* printf(" Unitarity: (%15.12f,%15.12f)\n", */
/* (v1^I(v2))/conj(v2^J(v1))-n*n*n); */
/* // Test of O */
/* v3=O(v1)-(pow(n*d,3)*v1); */
/* printf(" Test of O: (%15.12f,%15.12f)\n", */
/* v3^v3); */
/* // Test of Linv */
/* v1.dat[0][0][0]=Complex(0.,0.); /* Zero out G=0 component! */ */
/* v3=L(Linv(v1))-v1; */
/* printf(" Test of Linv: (%15.12f,%15.12f)\n", */
/* v3^v3); */
}