diff --git a/Introductory/mandelbroat/compare_languages/display.py b/Introductory/mandelbroat/compare_languages/display.py new file mode 100644 index 0000000..d201a41 --- /dev/null +++ b/Introductory/mandelbroat/compare_languages/display.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python +import sys +import numpy as np +import matplotlib.cm as cm +from matplotlib.colors import LogNorm +import pylab as pl + +filename = 'mand.dat' +if len(sys.argv)>1: + filename = sys.argv[1] + +X,Y,Z = np.loadtxt(filename,unpack=True) +ext=[X[0],X[-1],Y[0],Y[-1]] +N = int(np.sqrt(len(Z))) +dat = Z.reshape(N,N) +print(ext) +pl.imshow(dat.T,extent=ext,cmap=cm.hot,norm=LogNorm()) # pylab's function for displaying 2D image +pl.show() diff --git a/Introductory/mandelbroat/compare_languages/makefile b/Introductory/mandelbroat/compare_languages/makefile index 3770d43..3c08e7d 100644 --- a/Introductory/mandelbroat/compare_languages/makefile +++ b/Introductory/mandelbroat/compare_languages/makefile @@ -1,13 +1,14 @@ -CC = g++-10 +CC = g++-12 F90 = gfortran +CFLAGS = -O3 -fopenmp all : mandc mandf mandc : mandc.cc - $(CC) -O3 -fopenmp -o $@ $< + $(CC) $(CFLAGS) -o $@ $< mandf : mandf.f90 - $(F90) -O3 -fopenmp -o $@ $< + $(F90) $(CFLAGS) -o $@ $< clean : rm -f mandc mandf diff --git a/Introductory/mandelbroat/compare_languages/mandc.cc b/Introductory/mandelbroat/compare_languages/mandc.cc index fdedfe8..c170476 100644 --- a/Introductory/mandelbroat/compare_languages/mandc.cc +++ b/Introductory/mandelbroat/compare_languages/mandc.cc @@ -1,14 +1,15 @@ -#include -#include -#include -#include -#include -using namespace std; +#include // for printing to stdout +#include // for complex numbers +#include // for measuring time +#include // for using array/vector +#include // openMP and measuing time with openMP +using namespace std;// everything from standard template library will be in global scope +// std::cout and std::endl can be called cout, endl, etc... int Mandelb(const complex& z0, int max_steps) { - complex z=0; - for (int i=0; i z=0; // specify type where we start using it. (still need to specify type not like python with auto-typing) + for (int i=0; i2.) return i; z = z*z + z0; } @@ -17,35 +18,38 @@ int Mandelb(const complex& z0, int max_steps) int main() { - const int Nx = 1000; + const int Nx = 1000; // constants are named cost rather than parameter const int Ny = 1000; - int max_steps = 1000; - double ext[]={-2,1,-1,1}; + int max_steps = 1000; // allowed to change later. + double ext[]={-2,1,-1,1}; // this specifies fixed array of four numbers and initializes it. - vector mand(Nx*Ny); - clock_t startTimec = clock(); - double start = omp_get_wtime(); + vector mand(Nx*Ny); // allocating one dimensional array (vector) of size Nx*Ny + // multidimensional dynamic arrays are not standard in C++. One can use pointers to allocate/deallocate memory, + // and write one's own class interface for that. Or one has to use extension of C++ (blitz++ is excellent). + // Unfortunately the standard C++ still does not support standard class for that. + clock_t startTimec = clock(); // cpu time at the start + double start = omp_get_wtime(); // wall time at the start #pragma omp parallel for for (int i=0; i(x,y), max_steps); + double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.); // x in the interval ext[0]...ext[1] + double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.); // y in the interval ext[2]...ext[3] + mand[i*Ny+j] = Mandelb(complex(x,y), max_steps); // storing values in 2D array using 1D array } } clock_t endTimec = clock(); - double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC; - double diff = omp_get_wtime()-start; + double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC; // how to get seconds from cpu time + double diff = omp_get_wtime()-start; // openMP time is already in seconds - clog<<"clock time : "<2.) then - Mandelb = i-1 + if (abs(z)>2.) then ! |f(z)|>2 we know it will explode for large n, hence not part of the set. + Mandelb = i-1 ! return how many iterations it took to know this is not part of the set. return end if z = z*z + z0 @@ -31,7 +31,7 @@ program mand INTEGER, parameter :: max_steps = 1000 REAL*8 :: ext(4) = (/-2., 1., -1., 1./) ! The limits of plotting REAL*8 :: mande(Nx,Ny) - REAL :: start, finish, startw, finishw + REAL*8 :: start, finish, startw, finishw call cpu_time(start) startw = OMP_get_wtime() @@ -39,8 +39,8 @@ program mand !$OMP PARALLEL DO PRIVATE(j,x,y,z0) do i=1,Nx do j=1,Ny - x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.) - y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.) + x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.) ! x \in [ext(1)...ext(2)] + y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.) ! y \in [ext(3)...ext(4)] z0 = dcmplx(x,y) mande(i,j) = Mandelb(z0, max_steps) enddo diff --git a/Introductory/mandelbroat/compare_languages/manp.py b/Introductory/mandelbroat/compare_languages/manp.py index 7f86d16..bcd712d 100755 --- a/Introductory/mandelbroat/compare_languages/manp.py +++ b/Introductory/mandelbroat/compare_languages/manp.py @@ -1,11 +1,12 @@ #!/usr/bin/env python - -from scipy import * -from pylab import * -import time - +from scipy import * # scientific python +from pylab import * # plotting library +import time # timeing +from numba import jit # This is the new line with numba + +@jit(nopython=True) # This is the second new line with numba def Mand(z0, max_steps): - z = 0j + z = 0j # no need to specify type. To initialize to complex number, just assign 0j==i*0 for itr in range(max_steps): if abs(z)>2.: return itr @@ -18,17 +19,18 @@ def Mand(z0, max_steps): Ny = 1000 max_steps = 1000 #50 - ext = [-2,1,-1,1] + ext = [-2,1,-1,1] # no need to specify this is a python list, just initialize it t0 = time.time() - data = zeros( (Nx,Ny) ) + data = zeros( (Nx,Ny) ) # initialize a 2D dynamic array (something C++ is still lacking) for i in range(Nx): for j in range(Ny): x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.) y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.) - data[i,j] = Mand(x + y*1j, max_steps) + data[i,j] = Mand(x + y*1j, max_steps) # creating complex number of the fly + print ('clock time: '+str( time.time()-t0) ) - imshow(transpose(1./data), extent=ext) - show() + imshow(transpose(1./data), extent=ext) # pylab's function for displaying 2D image + show() # showing the image.