The blocking method was made popular by Flyvbjerg and Pedersen (1989)
-and has become one of the standard ways to estimate
-\( V(\widehat{\theta}) \) for exactly one \( \widehat{\theta} \), namely
-\( \widehat{\theta} = \overline{X} \).
+and has become one of the standard ways to estimate the variance
+\( \mathrm{var}(\widehat{\theta}) \) for exactly one estimator \( \widehat{\theta} \), namely
+\( \widehat{\theta} = \overline{X} \), the mean value.
Assume \( n = 2^d \) for some integer \( d>1 \) and \( X_1,X_2,\cdots, X_n \) is a stationary time series to begin with.
@@ -191,6 +687,9 @@
Flyvbjerg and Petersen demonstrated that the sequence
-\( \{e_k\}_{k=0}^{d-1} \) is decreasing, and conjecture that the term
-\( e_k \) can be made as small as we would like by making \( k \) (and hence
-\( d \)) sufficiently large. The sequence is decreasing (Master of Science thesis by Marius Jonsson, UiO 2018).
-It means we can apply blocking transformations until
-\( e_k \) is sufficiently small, and then estimate \( \mathrm{var}(\overline{X}) \) by
-\( \widehat{\sigma}^2_k/n_k \).
-
# 2-electron VMC code for 2dim quantum dot with importance sampling
-# Using gaussian rng for new positions and Metropolis- Hastings
-# Added energy minimization
-frommathimport exp, sqrt
-fromrandomimport random, seed, normalvariate
-importnumpyasnp
-importmatplotlib.pyplotasplt
-frommpl_toolkits.mplot3dimport Axes3D
-frommatplotlibimport cm
-frommatplotlib.tickerimport LinearLocator, FormatStrFormatter
-fromscipy.optimizeimport minimize
-importsys
-importos
-
-# Where to save data files
-PROJECT_ROOT_DIR ="Results"
-DATA_ID ="Results/EnergyMin"
+
More on the blocking method
-ifnot os.path.exists(PROJECT_ROOT_DIR):
- os.mkdir(PROJECT_ROOT_DIR)
-
-ifnot os.path.exists(DATA_ID):
- os.makedirs(DATA_ID)
-
-defdata_path(dat_id):
- return os.path.join(DATA_ID, dat_id)
-
-outfile =open(data_path("Energies.dat"),'w')
-
-
-# Trial wave function for the 2-electron quantum dot in two dims
-defWaveFunction(r,alpha,beta):
- r1 = r[0,0]**2+ r[0,1]**2
- r2 = r[1,0]**2+ r[1,1]**2
- r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
- deno = r12/(1+beta*r12)
- return exp(-0.5*alpha*(r1+r2)+deno)
-
-# Local energy for the 2-electron quantum dot in two dims, using analytical local energy
-defLocalEnergy(r,alpha,beta):
-
- r1 = (r[0,0]**2+ r[0,1]**2)
- r2 = (r[1,0]**2+ r[1,1]**2)
- r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
- deno =1.0/(1+beta*r12)
- deno2 = deno*deno
- return0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha +1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
-
-# Derivate of wave function ansatz as function of variational parameters
-defDerivativeWFansatz(r,alpha,beta):
-
- WfDer = np.zeros((2), np.double)
- r1 = (r[0,0]**2+ r[0,1]**2)
- r2 = (r[1,0]**2+ r[1,1]**2)
- r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
- deno =1.0/(1+beta*r12)
- deno2 = deno*deno
- WfDer[0] =-0.5*(r1+r2)
- WfDer[1] =-r12*r12*deno2
- return WfDer
-
-# Setting up the quantum force for the two-electron quantum dot, recall that it is a vector
-defQuantumForce(r,alpha,beta):
-
- qforce = np.zeros((NumberParticles,Dimension), np.double)
- r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
- deno =1.0/(1+beta*r12)
- qforce[0,:] =-2*r[0,:]*alpha*(r[0,:]-r[1,:])*deno*deno/r12
- qforce[1,:] =-2*r[1,:]*alpha*(r[1,:]-r[0,:])*deno*deno/r12
- return qforce
-
-
-# Computing the derivative of the energy and the energy
-defEnergyDerivative(x0):
-
-
- # Parameters in the Fokker-Planck simulation of the quantum force
- D =0.5
- TimeStep =0.05
- # positions
- PositionOld = np.zeros((NumberParticles,Dimension), np.double)
- PositionNew = np.zeros((NumberParticles,Dimension), np.double)
- # Quantum force
- QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
- QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
-
- energy =0.0
- DeltaE =0.0
- alpha = x0[0]
- beta = x0[1]
- EnergyDer =0.0
- DeltaPsi =0.0
- DerivativePsiE =0.0
- #Initial position
- for i inrange(NumberParticles):
- for j inrange(Dimension):
- PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
- wfold = WaveFunction(PositionOld,alpha,beta)
- QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
-
- #Loop over MC MCcycles
- for MCcycle inrange(NumberMCcycles):
- #Trial position moving one particle at the time
- for i inrange(NumberParticles):
- for j inrange(Dimension):
- PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
- QuantumForceOld[i,j]*TimeStep*D
- wfnew = WaveFunction(PositionNew,alpha,beta)
- QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
- GreensFunction =0.0
- for j inrange(Dimension):
- GreensFunction +=0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
- (D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
- PositionNew[i,j]+PositionOld[i,j])
-
- GreensFunction = exp(GreensFunction)
- ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
- #Metropolis-Hastings test to see whether we accept the move
- if random() <= ProbabilityRatio:
- for j inrange(Dimension):
- PositionOld[i,j] = PositionNew[i,j]
- QuantumForceOld[i,j] = QuantumForceNew[i,j]
- wfold = wfnew
- DeltaE = LocalEnergy(PositionOld,alpha,beta)
- DerPsi = DerivativeWFansatz(PositionOld,alpha,beta)
- DeltaPsi += DerPsi
- energy += DeltaE
- DerivativePsiE += DerPsi*DeltaE
-
- # We calculate mean values
- energy /= NumberMCcycles
- DerivativePsiE /= NumberMCcycles
- DeltaPsi /= NumberMCcycles
- EnergyDer =2*(DerivativePsiE-DeltaPsi*energy)
- return EnergyDer
-
-
-# Computing the expectation value of the local energy
-defEnergy(x0):
- # Parameters in the Fokker-Planck simulation of the quantum force
- D =0.5
- TimeStep =0.05
- # positions
- PositionOld = np.zeros((NumberParticles,Dimension), np.double)
- PositionNew = np.zeros((NumberParticles,Dimension), np.double)
- # Quantum force
- QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
- QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
-
- energy =0.0
- DeltaE =0.0
- alpha = x0[0]
- beta = x0[1]
- #Initial position
- for i inrange(NumberParticles):
- for j inrange(Dimension):
- PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
- wfold = WaveFunction(PositionOld,alpha,beta)
- QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
-
- #Loop over MC MCcycles
- for MCcycle inrange(NumberMCcycles):
- #Trial position moving one particle at the time
- for i inrange(NumberParticles):
- for j inrange(Dimension):
- PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
- QuantumForceOld[i,j]*TimeStep*D
- wfnew = WaveFunction(PositionNew,alpha,beta)
- QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
- GreensFunction =0.0
- for j inrange(Dimension):
- GreensFunction +=0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
- (D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
- PositionNew[i,j]+PositionOld[i,j])
-
- GreensFunction = exp(GreensFunction)
- ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
- #Metropolis-Hastings test to see whether we accept the move
- if random() <= ProbabilityRatio:
- for j inrange(Dimension):
- PositionOld[i,j] = PositionNew[i,j]
- QuantumForceOld[i,j] = QuantumForceNew[i,j]
- wfold = wfnew
- DeltaE = LocalEnergy(PositionOld,alpha,beta)
- energy += DeltaE
- if Printout:
- outfile.write('%f\n'%(energy/(MCcycle+1.0)))
- # We calculate mean values
- energy /= NumberMCcycles
- return energy
-
-#Here starts the main program with variable declarations
-NumberParticles =2
-Dimension =2
-# seed for rng generator
-seed()
-# Monte Carlo cycles for parameter optimization
-Printout =False
-NumberMCcycles=10000
-# guess for variational parameters
-x0 = np.array([0.9,0.2])
-# Using Broydens method to find optimal parameters
-res = minimize(Energy, x0, method='BFGS', jac=EnergyDerivative, options={'gtol': 1e-4,'disp': True})
-x0 = res.x
-# Compute the energy again with the optimal parameters and increased number of Monte Cycles
-NumberMCcycles=2**19
-Printout =True
-FinalEnergy = Energy(x0)
-EResult = np.array([FinalEnergy,FinalEnergy])
-outfile.close()
-#nice printout with Pandas
-importpandasaspd
-frompandasimport DataFrame
-data ={'Optimal Parameters':x0, 'Final Energy':EResult}
-frame = pd.DataFrame(data)
-print(frame)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
Flyvbjerg and Petersen demonstrated that the sequence
+\( \{e_k\}_{k=0}^{d-1} \) is decreasing, and conjecture that the term
+\( e_k \) can be made as small as we would like by making \( k \) (and hence
+\( d \)) sufficiently large. The sequence is decreasing.
+It means we can apply blocking transformations until
+\( e_k \) is sufficiently small, and then estimate \( \mathrm{var}(\overline{X}) \) by
+\( \widehat{\sigma}^2_k/n_k \).
+
The next step is then to use the above data sets and perform a
-resampling analysis using the blocking method
-The blocking code, based on the article of Marius Jonsson is given here
-
-
+
Example code form last week
@@ -164,60 +654,221 @@
Resampling analysis
-
# Common imports
+
# 2-electron VMC code for 2dim quantum dot with importance sampling
+# Using gaussian rng for new positions and Metropolis- Hastings
+# Added energy minimization
+frommathimport exp, sqrt
+fromrandomimport random, seed, normalvariate
+importnumpyasnp
+importmatplotlib.pyplotasplt
+frommpl_toolkits.mplot3dimport Axes3D
+frommatplotlibimport cm
+frommatplotlib.tickerimport LinearLocator, FormatStrFormatter
+fromscipy.optimizeimport minimize
+importsysimportos
-# Where to save the figures and data files
+# Where to save data files
+PROJECT_ROOT_DIR ="Results"
DATA_ID ="Results/EnergyMin"
+ifnot os.path.exists(PROJECT_ROOT_DIR):
+ os.mkdir(PROJECT_ROOT_DIR)
+
+ifnot os.path.exists(DATA_ID):
+ os.makedirs(DATA_ID)
+
defdata_path(dat_id):
return os.path.join(DATA_ID, dat_id)
-infile =open(data_path("Energies.dat"),'r')
-
-fromnumpyimport log2, zeros, mean, var, sum, loadtxt, arange, array, cumsum, dot, transpose, diagonal, sqrt
-fromnumpy.linalgimport inv
-
-defblock(x):
- # preliminaries
- n =len(x)
- d =int(log2(n))
- s, gamma = zeros(d), zeros(d)
- mu = mean(x)
-
- # estimate the auto-covariance and variances
- # for each blocking transformation
- for i in arange(0,d):
- n =len(x)
- # estimate autocovariance of x
- gamma[i] = (n)**(-1)*sum( (x[0:(n-1)]-mu)*(x[1:n]-mu) )
- # estimate variance of x
- s[i] = var(x)
- # perform blocking transformation
- x =0.5*(x[0::2] + x[1::2])
-
- # generate the test observator M_k from the theorem
- M = (cumsum( ((gamma/s)**2*2**arange(1,d+1)[::-1])[::-1] ) )[::-1]
-
- # we need a list of magic numbers
- q =array([6.634897,9.210340, 11.344867, 13.276704, 15.086272, 16.811894, 18.475307, 20.090235, 21.665994, 23.209251, 24.724970, 26.216967, 27.688250, 29.141238, 30.577914, 31.999927, 33.408664, 34.805306, 36.190869, 37.566235, 38.932173, 40.289360, 41.638398, 42.979820, 44.314105, 45.641683, 46.962942, 48.278236, 49.587884, 50.892181])
-
- # use magic to determine when we should have stopped blocking
- for k in arange(0,d):
- if(M[k] < q[k]):
- break
- if (k >= d-1):
- print("Warning: Use more data")
- return mu, s[k]/2**(d-k)
-
-
-x = loadtxt(infile)
-(mean, var) = block(x)
-std = sqrt(var)
+outfile =open(data_path("Energies.dat"),'w')
+
+
+# Trial wave function for the 2-electron quantum dot in two dims
+defWaveFunction(r,alpha,beta):
+ r1 = r[0,0]**2+ r[0,1]**2
+ r2 = r[1,0]**2+ r[1,1]**2
+ r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
+ deno = r12/(1+beta*r12)
+ return exp(-0.5*alpha*(r1+r2)+deno)
+
+# Local energy for the 2-electron quantum dot in two dims, using analytical local energy
+defLocalEnergy(r,alpha,beta):
+
+ r1 = (r[0,0]**2+ r[0,1]**2)
+ r2 = (r[1,0]**2+ r[1,1]**2)
+ r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
+ deno =1.0/(1+beta*r12)
+ deno2 = deno*deno
+ return0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha +1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
+
+# Derivate of wave function ansatz as function of variational parameters
+defDerivativeWFansatz(r,alpha,beta):
+
+ WfDer = np.zeros((2), np.double)
+ r1 = (r[0,0]**2+ r[0,1]**2)
+ r2 = (r[1,0]**2+ r[1,1]**2)
+ r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
+ deno =1.0/(1+beta*r12)
+ deno2 = deno*deno
+ WfDer[0] =-0.5*(r1+r2)
+ WfDer[1] =-r12*r12*deno2
+ return WfDer
+
+# Setting up the quantum force for the two-electron quantum dot, recall that it is a vector
+defQuantumForce(r,alpha,beta):
+
+ qforce = np.zeros((NumberParticles,Dimension), np.double)
+ r12 = sqrt((r[0,0]-r[1,0])**2+ (r[0,1]-r[1,1])**2)
+ deno =1.0/(1+beta*r12)
+ qforce[0,:] =-2*r[0,:]*alpha*(r[0,:]-r[1,:])*deno*deno/r12
+ qforce[1,:] =-2*r[1,:]*alpha*(r[1,:]-r[0,:])*deno*deno/r12
+ return qforce
+
+
+# Computing the derivative of the energy and the energy
+defEnergyDerivative(x0):
+
+
+ # Parameters in the Fokker-Planck simulation of the quantum force
+ D =0.5
+ TimeStep =0.05
+ # positions
+ PositionOld = np.zeros((NumberParticles,Dimension), np.double)
+ PositionNew = np.zeros((NumberParticles,Dimension), np.double)
+ # Quantum force
+ QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
+ QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
+
+ energy =0.0
+ DeltaE =0.0
+ alpha = x0[0]
+ beta = x0[1]
+ EnergyDer =0.0
+ DeltaPsi =0.0
+ DerivativePsiE =0.0
+ #Initial position
+ for i inrange(NumberParticles):
+ for j inrange(Dimension):
+ PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
+ wfold = WaveFunction(PositionOld,alpha,beta)
+ QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
+
+ #Loop over MC MCcycles
+ for MCcycle inrange(NumberMCcycles):
+ #Trial position moving one particle at the time
+ for i inrange(NumberParticles):
+ for j inrange(Dimension):
+ PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
+ QuantumForceOld[i,j]*TimeStep*D
+ wfnew = WaveFunction(PositionNew,alpha,beta)
+ QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
+ GreensFunction =0.0
+ for j inrange(Dimension):
+ GreensFunction +=0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
+ (D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
+ PositionNew[i,j]+PositionOld[i,j])
+
+ GreensFunction = exp(GreensFunction)
+ ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
+ #Metropolis-Hastings test to see whether we accept the move
+ if random() <= ProbabilityRatio:
+ for j inrange(Dimension):
+ PositionOld[i,j] = PositionNew[i,j]
+ QuantumForceOld[i,j] = QuantumForceNew[i,j]
+ wfold = wfnew
+ DeltaE = LocalEnergy(PositionOld,alpha,beta)
+ DerPsi = DerivativeWFansatz(PositionOld,alpha,beta)
+ DeltaPsi += DerPsi
+ energy += DeltaE
+ DerivativePsiE += DerPsi*DeltaE
+
+ # We calculate mean values
+ energy /= NumberMCcycles
+ DerivativePsiE /= NumberMCcycles
+ DeltaPsi /= NumberMCcycles
+ EnergyDer =2*(DerivativePsiE-DeltaPsi*energy)
+ return EnergyDer
+
+
+# Computing the expectation value of the local energy
+defEnergy(x0):
+ # Parameters in the Fokker-Planck simulation of the quantum force
+ D =0.5
+ TimeStep =0.05
+ # positions
+ PositionOld = np.zeros((NumberParticles,Dimension), np.double)
+ PositionNew = np.zeros((NumberParticles,Dimension), np.double)
+ # Quantum force
+ QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
+ QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
+
+ energy =0.0
+ DeltaE =0.0
+ alpha = x0[0]
+ beta = x0[1]
+ #Initial position
+ for i inrange(NumberParticles):
+ for j inrange(Dimension):
+ PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
+ wfold = WaveFunction(PositionOld,alpha,beta)
+ QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
+
+ #Loop over MC MCcycles
+ for MCcycle inrange(NumberMCcycles):
+ #Trial position moving one particle at the time
+ for i inrange(NumberParticles):
+ for j inrange(Dimension):
+ PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
+ QuantumForceOld[i,j]*TimeStep*D
+ wfnew = WaveFunction(PositionNew,alpha,beta)
+ QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
+ GreensFunction =0.0
+ for j inrange(Dimension):
+ GreensFunction +=0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
+ (D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
+ PositionNew[i,j]+PositionOld[i,j])
+
+ GreensFunction = exp(GreensFunction)
+ ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
+ #Metropolis-Hastings test to see whether we accept the move
+ if random() <= ProbabilityRatio:
+ for j inrange(Dimension):
+ PositionOld[i,j] = PositionNew[i,j]
+ QuantumForceOld[i,j] = QuantumForceNew[i,j]
+ wfold = wfnew
+ DeltaE = LocalEnergy(PositionOld,alpha,beta)
+ energy += DeltaE
+ if Printout:
+ outfile.write('%f\n'%(energy/(MCcycle+1.0)))
+ # We calculate mean values
+ energy /= NumberMCcycles
+ return energy
+
+#Here starts the main program with variable declarations
+NumberParticles =2
+Dimension =2
+# seed for rng generator
+seed()
+# Monte Carlo cycles for parameter optimization
+Printout =False
+NumberMCcycles=10000
+# guess for variational parameters
+x0 = np.array([0.9,0.2])
+# Using Broydens method to find optimal parameters
+res = minimize(Energy, x0, method='BFGS', jac=EnergyDerivative, options={'gtol': 1e-4,'disp': True})
+x0 = res.x
+# Compute the energy again with the optimal parameters and increased number of Monte Cycles
+NumberMCcycles=2**19
+Printout =True
+FinalEnergy = Energy(x0)
+EResult = np.array([FinalEnergy,FinalEnergy])
+outfile.close()
+#nice printout with Pandasimportpandasaspdfrompandasimport DataFrame
-data ={'Mean':[mean], 'STDev':[std]}
-frame = pd.DataFrame(data,index=['Values'])
+data ={'Optimal Parameters':x0, 'Final Energy':EResult}
+frame = pd.DataFrame(data)
print(frame)
The blocking method was made popular by Flyvbjerg and Pedersen (1989)
-and has become one of the standard ways to estimate
-\( V(\widehat{\theta}) \) for exactly one \( \widehat{\theta} \), namely
-\( \widehat{\theta} = \overline{X} \).
+
The next step is then to use the above data sets and perform a
+resampling analysis using the blocking method
+The blocking code, based on the article of Marius Jonsson is given here
-
Assume \( n = 2^d \) for some integer \( d>1 \) and \( X_1,X_2,\cdots, X_n \) is a stationary time series to begin with.
-Moreover, assume that the time series is asymptotically uncorrelated. We switch to vector notation by arranging \( X_1,X_2,\cdots,X_n \) in an \( n \)-tuple. Define:
-
The strength of the blocking method is when the number of
-observations, \( n \) is large. For large \( n \), the complexity of dependent
-bootstrapping scales poorly, but the blocking method does not,
-moreover, it becomes more accurate the larger \( n \) is.
-
+
+
+
+
+
+
+
+
# Common imports
+importos
+
+# Where to save the figures and data files
+DATA_ID ="Results/EnergyMin"
+
+defdata_path(dat_id):
+ return os.path.join(DATA_ID, dat_id)
+
+infile =open(data_path("Energies.dat"),'r')
+
+fromnumpyimport log2, zeros, mean, var, sum, loadtxt, arange, array, cumsum, dot, transpose, diagonal, sqrt
+fromnumpy.linalgimport inv
+
+defblock(x):
+ # preliminaries
+ n =len(x)
+ d =int(log2(n))
+ s, gamma = zeros(d), zeros(d)
+ mu = mean(x)
+
+ # estimate the auto-covariance and variances
+ # for each blocking transformation
+ for i in arange(0,d):
+ n =len(x)
+ # estimate autocovariance of x
+ gamma[i] = (n)**(-1)*sum( (x[0:(n-1)]-mu)*(x[1:n]-mu) )
+ # estimate variance of x
+ s[i] = var(x)
+ # perform blocking transformation
+ x =0.5*(x[0::2] + x[1::2])
+
+ # generate the test observator M_k from the theorem
+ M = (cumsum( ((gamma/s)**2*2**arange(1,d+1)[::-1])[::-1] ) )[::-1]
+
+ # we need a list of magic numbers
+ q =array([6.634897,9.210340, 11.344867, 13.276704, 15.086272, 16.811894, 18.475307, 20.090235, 21.665994, 23.209251, 24.724970, 26.216967, 27.688250, 29.141238, 30.577914, 31.999927, 33.408664, 34.805306, 36.190869, 37.566235, 38.932173, 40.289360, 41.638398, 42.979820, 44.314105, 45.641683, 46.962942, 48.278236, 49.587884, 50.892181])
+
+ # use magic to determine when we should have stopped blocking
+ for k in arange(0,d):
+ if(M[k] < q[k]):
+ break
+ if (k >= d-1):
+ print("Warning: Use more data")
+ return mu, s[k]/2**(d-k)
+
+
+x = loadtxt(infile)
+(mean, var) = block(x)
+std = sqrt(var)
+importpandasaspd
+frompandasimport DataFrame
+data ={'Mean':[mean], 'STDev':[std]}
+frame = pd.DataFrame(data,index=['Values'])
+print(frame)
+
We now define
-blocking transformations. The idea is to take the mean of subsequent
-pair of elements from \( \vec{X} \) and form a new vector
-\( \vec{X}_1 \). Continuing in the same way by taking the mean of
-subsequent pairs of elements of \( \vec{X}_1 \) we obtain \( \vec{X}_2 \), and
-so on.
-Define \( \vec{X}_i \) recursively by:
-
The quantity \( \vec{X}_k \) is
-subject to \( k \) blocking transformations. We now have \( d \) vectors
-\( \vec{X}_0, \vec{X}_1,\cdots,\vec X_{d-1} \) containing the subsequent
-averages of observations. It turns out that if the components of
-\( \vec{X} \) is a stationary time series, then the components of
-\( \vec{X}_i \) is a stationary time series for all \( 0 \leq i \leq d-1 \)
-
-
-
We can then compute the autocovariance, the variance, sample mean, and
-number of observations for each \( i \).
-Let \( \gamma_i, \sigma_i^2,
-\overline{X}_i \) denote the autocovariance, variance and average of the
-elements of \( \vec{X}_i \) and let \( n_i \) be the number of elements of
-\( \vec{X}_i \). It follows by induction that \( n_i = n/2^i \).
-
Till now we have not paid much attention to speed and possible optimization possibilities
+inherent in the various compilers. We have compiled and linked as
+
+
+
+
+
+
+
+
+
c++-cmycode.cpp
+c++-omycode.exemycode.o
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
Using the
-definition of the blocking transformation and the distributive
-property of the covariance, it is clear that since \( h =|i-j| \)
-we can define
+
For Fortran replace with for example gfortran or ifort.
+This is what we call a flat compiler option and should be used when we develop the code.
+It produces normally a very large and slow code when translated to machine instructions.
+We use this option for debugging and for establishing the correct program output because
+every operation is done precisely as the user specified it.
The quantity \( \hat{X} \) is asymptotic uncorrelated by assumption, \( \hat{X}_k \) is also asymptotic uncorrelated. Let's turn our attention to the variance of the sample mean \( V(\overline{X}) \).
+
It is instructive to look up the compiler manual for further instructions by writing
We have additional compiler options for optimization. These may include procedure inlining where
+performance may be improved, moving constants inside loops outside the loop,
+identify potential parallelism, include automatic vectorization or replace a division with a reciprocal
+and a multiplication if this speeds up the code.
+
-
The term \( e_k \) is called the truncation error:
After you have run the code you can obtain the profiling information via
-
By repeated use of this equation we get \( V(\overline{X}_i) = V(\overline{X}_0) = V(\overline{X}) \) for all \( 0 \leq i \leq d-1 \). This has the consequence that
Flyvbjerg and Petersen demonstrated that the sequence
-\( \{e_k\}_{k=0}^{d-1} \) is decreasing, and conjecture that the term
-\( e_k \) can be made as small as we would like by making \( k \) (and hence
-\( d \)) sufficiently large. The sequence is decreasing (Master of Science thesis by Marius Jonsson, UiO 2018).
-It means we can apply blocking transformations until
-\( e_k \) is sufficiently small, and then estimate \( V(\overline{X}) \) by
-\( \widehat{\sigma}^2_k/n_k \).
+
When you have profiled properly your code, you must take out this option as it
+slows down performance.
+For memory tests use valgrind. An excellent environment for all these aspects, and much more, is Qt creator.
This option generates debugging information allowing you to trace for example if an array is properly allocated. Some compilers work best with the no optimization option -O0.
+
+
+
+
+
+
+
Depending on the compiler, one can add flags which generate code that catches integer overflow errors.
+The flag -ftrapv does this for the CLANG compiler on OS X operating systems.
+
The next step is then to use the above data sets and perform a
-resampling analysis using the blocking method
-The blocking code, based on the article of Marius Jonsson is given here
-
-
+
Other hints
+
+
+
+
In general, irrespective of compiler options, it is useful to
+
+
avoid if tests or call to functions inside loops, if possible.
+
avoid multiplication with constants inside loops if possible
+
+
Here is an example of a part of a program where specific operations lead to a slower code
-
+
-
# Common imports
-importos
-
-# Where to save the figures and data files
-DATA_ID ="Results/EnergyMin"
-
-defdata_path(dat_id):
- return os.path.join(DATA_ID, dat_id)
-
-infile =open(data_path("Energies.dat"),'r')
-
-fromnumpyimport log2, zeros, mean, var, sum, loadtxt, arange, array, cumsum, dot, transpose, diagonal, sqrt
-fromnumpy.linalgimport inv
-
-defblock(x):
- # preliminaries
- n =len(x)
- d =int(log2(n))
- s, gamma = zeros(d), zeros(d)
- mu = mean(x)
-
- # estimate the auto-covariance and variances
- # for each blocking transformation
- for i in arange(0,d):
- n =len(x)
- # estimate autocovariance of x
- gamma[i] = (n)**(-1)*sum( (x[0:(n-1)]-mu)*(x[1:n]-mu) )
- # estimate variance of x
- s[i] = var(x)
- # perform blocking transformation
- x =0.5*(x[0::2] + x[1::2])
-
- # generate the test observator M_k from the theorem
- M = (cumsum( ((gamma/s)**2*2**arange(1,d+1)[::-1])[::-1] ) )[::-1]
-
- # we need a list of magic numbers
- q =array([6.634897,9.210340, 11.344867, 13.276704, 15.086272, 16.811894, 18.475307, 20.090235, 21.665994, 23.209251, 24.724970, 26.216967, 27.688250, 29.141238, 30.577914, 31.999927, 33.408664, 34.805306, 36.190869, 37.566235, 38.932173, 40.289360, 41.638398, 42.979820, 44.314105, 45.641683, 46.962942, 48.278236, 49.587884, 50.892181])
-
- # use magic to determine when we should have stopped blocking
- for k in arange(0,d):
- if(M[k] < q[k]):
- break
- if (k >= d-1):
- print("Warning: Use more data")
- return mu, s[k]/2**(d-k)
+
Here we avoid a repeated multiplication inside a loop.
+Most compilers, depending on compiler flags, identify and optimize such bottlenecks on their own, without requiring any particular action by the programmer. However, it is always useful to single out and avoid code examples like the first one discussed here.
+
Vectorization and the basic idea behind parallel computing
+
+
+
+
Present CPUs are highly parallel processors with varying levels of parallelism. The typical situation can be described via the following three statements.
+
+
Pursuit of shorter computation time and larger simulation size gives rise to parallel computing.
+
Multiple processors are involved to solve a global problem.
+
The essence is to divide the entire computation evenly among collaborative processors. Divide and conquer.
+
+
Before we proceed with a more detailed discussion of topics like vectorization and parallelization, we need to remind ourselves about some basic features of different hardware models.
+
+
-
Automatic vectorization and vectorization inhibitors, memory stride
-
-
-For C++ programmers it is also worth keeping in mind that an array notation is preferred to the more compact use of pointers to access array elements. The compiler can often not tell if it is safe to vectorize the code.
-
-
-When dealing with arrays, you should also avoid memory stride, since this slows down considerably vectorization. When you access array element, write for example the inner loop to vectorize using unit stride, that is, access successively the next array element in memory, as shown here
-
-
-
for (int i =0; i < n; i++) {
- for (int j =0; j < n; j++) {
- a[i][j] += b[i][j];
- }
- }
-
Cache memory contains a copy of the main memory data
-
Cache is faster but consumes more space and power. It is normally assumed to be much faster than main memory
-
Registers contain working data only
-
-
-
Modern CPUs perform most or all operations only on data in register
-
-
-
Multiple Cache memories contain a copy of the main memory data
+
A rough classification of hardware models
+
+
+
-
Cache items accessed by their address in main memory
-
L1 cache is the fastest but has the least capacity
-
L2, L3 provide intermediate performance/size tradeoffs
+
Conventional single-processor computers are named SISD (single-instruction-single-data) machines.
+
SIMD (single-instruction-multiple-data) machines incorporate the idea of parallel processing, using a large number of processing units to execute the same instruction on different data.
+
Modern parallel computers are so-called MIMD (multiple-instruction-multiple-data) machines and can execute different instruction streams in parallel on different data.
+
+
-
-
-Loads and stores to memory can be as important as floating point operations when we measure performance.
-
-
-How do we measure performance? What is wrong with this code to time a loop?
-
+
+
Task parallelism: the work of a global problem can be divided into a number of independent tasks, which rarely need to synchronize. Monte Carlo simulations represent a typical situation. Integration is another. However this paradigm is of limited use.
+
Data parallelism: use of multiple threads (e.g. one or more threads per processor) to dissect loops over arrays etc. Communication and synchronization between processors are often hidden, thus easy to program. However, the user surrenders much control to a specialized compiler. Examples of data parallelism are compiler-based parallelization and OpenMP directives.
Message passing: all involved processors have an independent memory address space. The user is responsible for partitioning the data/work of a global problem and distributing the subproblems to the processors. Collaboration between processors is achieved by explicit message passing, which is used for data transfer plus synchronization.
+
This paradigm is the most general one where the user has full control. Better parallel efficiency is usually achieved by explicit message passing. However, message-passing programming is more difficult.
+
+
+
-
-
Timers are not infinitely accurate
-
All clocks have a granularity, the minimum time that they can measure
-
The error in a time measurement, even if everything is perfect, may be the size of this granularity (sometimes called a clock tick)
-
Always know what your clock granularity is
-
Ensure that your measurement is for a long enough duration (say 100 times the tick)
Vectorization is a special
+case of Single Instructions Multiple Data (SIMD) to denote a single
+instruction stream capable of operating on multiple data elements in
+parallel.
+We can think of vectorization as the unrolling of loops accompanied with SIMD instructions.
+
-
Problems with cold start
+
Vectorization is the process of converting an algorithm that performs scalar operations
+(typically one operation at the time) to vector operations where a single operation can refer to many simultaneous operations.
+Consider the following example
+
-
-What happens when the code is executed? The assumption is that the code is ready to
-execute. But
+
+
+
+
+
+
+
for(i=0;i<n;i++){
+a[i]=b[i]+c[i];
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
Code may still be on disk, and not even read into memory.
-
Data may be in slow memory rather than fast (which may be wrong or right for what you are measuring)
-
Multiple tests often necessary to ensure that cold start effects are not present
-
Special effort often required to ensure data in the intended part of the memory hierarchy.
-
+
If the code is not vectorized, the compiler will simply start with the first element and
+then perform subsequent additions operating on one address in memory at the time.
+
A SIMD instruction can operate on multiple data elements in one single instruction.
+It uses the so-called 128-bit SIMD floating-point register.
+In this sense, vectorization adds some form of parallelism since one instruction is applied
+to many parts of say a vector.
+
-
-
If the result of the computation is not used, the compiler may eliminate the code
-
Performance will look impossibly fantastic
-
Even worse, eliminate some of the code so the performance looks plausible
-
Ensure that the results are (or may be) used.
-
+
The number of elements which can be operated on in parallel
+range from four single-precision floating point data elements in so-called
+Streaming SIMD Extensions and two double-precision floating-point data
+elements in Streaming SIMD Extensions 2 to sixteen byte operations in
+a 128-bit register in Streaming SIMD Extensions 2. Thus, vector-length
+ranges from 2 to 16, depending on the instruction extensions used and
+on the data type.
+
+
IN summary, our instructions operate on 128 bit (16 byte) operands
Some parts of the hardware do not always perform with exactly the same performance
-
-
-
Make multiple tests and report
-
Easy choices include
+
+
Number of elements that can acted upon, examples
+
We start with the simple scalar operations given by
+
+
+
+
+
+
+
+
for(i=0;i<n;i++){
+a[i]=b[i]+c[i];
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
Average tests represent what users might observe over time
-
+
If the code is not vectorized and we have a 128-bit register to store a 32 bits floating point number,
+it means that we have \( 3\times 32 \) bits that are not used.
+
-
+
We have thus unused space in our SIMD registers. These registers could hold three additional integers.
-The next step is the backward substitution step. The last row is multiplied by \( c_{N-3}/b_{N-2} \) and subtracted from the second to last row, thus eliminating \( c_{N-3} \) from the last row. The general backward substitution procedure is
-$$
- c_i = 0,
-$$
-
-and
-$$
- f_{i-1} = f_{i-1} - \frac{c_{i-1}}{b_i}f_i
-$$
-
-All that ramains to be computed is the solution, which is the very straight forward process of
-$$
-x_i = \frac{f_i}{b_i}
-$$
+
The speedup depends on the size of the vectors. In the example here we have run with \( 10^7 \) elements.
+The example here was run on an IMac17.1 with OSX El Capitan (10.11.4) as operating system and an Intel i5 3.3 GHz CPU.
+
This particular C++ compiler speeds up the above loop operations with a factor of 1.5
+Performing the same operations for \( 10^9 \) elements results in a smaller speedup since reading from main memory is required. The non-vectorized code is seemingly faster.
+
-
// Forward substitution
-// Note that we can simplify by precalculating a[i-1]/b[i-1]
- for (int i=1; i < n; i++) {
- b[i] = b[i] - (a[i-1]*c[i-1])/b[i-1];
- f[i] = g[i] - (a[i-1]*f[i-1])/b[i-1];
- }
- x[n-1] = f[n-1] / b[n-1];
- // Backwards substitution
- for (int i = n-2; i >=0; i--) {
- f[i] = f[i] - c[i]*f[i+1]/b[i+1];
- x[i] = f[i]/b[i];
- }
-
We can also add vectorization analysis, see for example
-using namespace std; // note use of namespace
-int main (int argc, char* argv[])
-{
- // read in dimension of square matrix
- int n = atoi(argv[1]);
- double **A, **B;
- // Allocate space for the two matrices
- A = new double*[n]; B = new double*[n];
- for (int i = 0; i < n; i++){
- A[i] = new double[n];
- B[i] = new double[n];
- }
- // Set up values for matrix A
- for (int i = 0; i < n; i++){
- for (int j = 0; j < n; j++) {
- A[i][j] = cos(i*1.0)*sin(j*3.0);
- }
- }
- clock_t start, finish;
- start = clock();
- // Then compute the transpose
- for (int i = 0; i < n; i++){
- for (int j = 0; j < n; j++) {
- B[i][j]= A[j][i];
- }
- }
+
+
The variable \( n \) does need to be known at compile time. However, this variable must stay the same for the entire duration of the loop. It implies that an exit statement inside the loop cannot be data dependent.
-using namespace std; // note use of namespace
-int main (int argc, char* argv[])
-{
- // read in dimension of square matrix
- int n = atoi(argv[1]);
- double s = 1.0/sqrt( (double) n);
- double **A, **B, **C;
- // Start timing
- clock_t start, finish;
- start = clock();
- // Allocate space for the two matrices
- A = new double*[n]; B = new double*[n]; C = new double*[n];
- for (int i = 0; i < n; i++){
- A[i] = new double[n];
- B[i] = new double[n];
- C[i] = new double[n];
- }
- // Set up values for matrix A and B and zero matrix C
- for (int i = 0; i < n; i++){
- for (int j = 0; j < n; j++) {
- double angle = 2.0*M_PI*i*j/ (( double ) n);
- A[i][j] = s * ( sin ( angle ) + cos ( angle ) );
- B[j][i] = A[i][j];
- }
- }
- // Then perform the matrix-matrix multiplication
- for (int i = 0; i < n; i++){
- for (int j = 0; j < n; j++) {
- double sum = 0.0;
- for (int k = 0; k < n; k++) {
- sum += B[i][k]*A[k][j];
- }
- C[i][j] = sum;
- }
- }
- // Compute now the Frobenius norm
- double Fsum = 0.0;
- for (int i = 0; i < n; i++){
- for (int j = 0; j < n; j++) {
- Fsum += C[i][j]*C[i][j];
- }
- }
- Fsum = sqrt(Fsum);
- finish = clock();
- double timeused = (double) (finish - start)/(CLOCKS_PER_SEC );
- cout << setiosflags(ios::showpoint | ios::uppercase);
- cout << setprecision(10) << setw(20) << "Time used for matrix-matrix multiplication=" << timeused << endl;
- cout << " Frobenius norm = " << Fsum << endl;
- // Free up space
- for (int i = 0; i < n; i++){
- delete[] A[i];
- delete[] B[i];
- delete[] C[i];
- }
- delete[] A;
- delete[] B;
- delete[] C;
- return 0;
-}
-
Automatic vectorization and vectorization inhibitors, exit criteria
-
How do we define speedup? Simplest form
-
-
-
+
An exit statement should in general be avoided.
+If the exit statement contains data-dependent conditions, the loop cannot be vectorized.
+The following is an example of a non-vectorizable loop
+
-
-
Speedup measures the ratio of performance between two objects
-
Versions of same code, with different number of processors
-
Serial and vector versions
-
Try different programing languages, c++ and Fortran
-The key is choosing the correct baseline for comparison
-
-
-
For our serial vs. vectorization examples, using compiler-provided vectorization, the baseline is simple; the same code, with vectorization turned off
-
-
-
For parallel applications, this is much harder:
-
-
-
Choice of algorithm, decomposition, performance of baseline case etc.
-
-
-
-
-
+
Automatic vectorization and vectorization inhibitors, straight-line code
+
+
SIMD instructions perform the same type of operations multiple times.
+A switch statement leads thus to a non-vectorizable loop since different statemens cannot branch.
+The following code can however be vectorized since the if statement is implemented as a masked assignment.
+
These operations can be performed for all data elements but only those elements which the mask evaluates as true are stored. In general, one should avoid branches such as switch, go to, or return statements or if constructs that cannot be treated as masked assignments.
Automatic vectorization and vectorization inhibitors, function calls
-
Speedup and memory
-
-
-
-The speedup on \( p \) processors can
-be greater than \( p \) if memory usage is optimal!
-Consider the case of a memorybound computation with \( M \) words of memory
+
Calls to programmer defined functions ruin vectorization. However, calls to intrinsic functions like
+\( \sin{x} \), \( \cos{x} \), \( \exp{x} \) etc are allowed since they are normally efficiently vectorized.
+The following example is fully vectorizable
+
-
-
If \( M/p \) fits into cache while \( M \) does not, the time to access memory will be different in the two cases:
-
\( T_1 \) uses the main memory bandwidth
-
\( T_p \) uses the appropriate cache bandwidth
-
+
+
+
+
+
+
+
for(inti=0;i<n;i++){
+a[i]=log10(i)*cos(i);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
Similarly, inline functions defined by the programmer, allow for vectorization since the function statements are glued into the actual place where the function is called.
Automatic vectorization and vectorization inhibitors, data dependencies
+
+
One has to keep in mind that vectorization changes the order of operations inside a loop. A so-called
+read-after-write statement with an explicit flow dependency cannot be vectorized. The following code
+
-Assume that almost all parts of a code are perfectly
-parallelizable (fraction \( f \)). The remainder,
-fraction \( (1-f) \) cannot be parallelized at all.
-
-
-That is, there is work that takes time \( W \) on one process; a fraction \( f \) of that work will take
-time \( Wf/p \) on \( p \) processors.
-
-
-
What is the maximum possible speedup as a function of \( f \)?
-
+
is an example of flow dependency and results in wrong numerical results if vectorized. For a scalar operation, the value \( a[i-1] \) computed during the iteration is loaded into the right-hand side and the results are fine. In vector mode however, with a vector length of four, the values \( a[0] \), \( a[1] \), \( a[2] \) and \( a[3] \) from the previous loop will be loaded into the right-hand side and produce wrong results. That is, we have
and if the two first iterations are executed at the same by the SIMD instruction, the value of say \( a[1] \) could be used by the second iteration before it has been calculated by the first iteration, leading thereby to wrong results.
-On one processor we have
-$$
-T_1 = (1-f)W + fW = W
-$$
-
-On \( p \) processors we have
-$$
-T_p = (1-f)W + \frac{fW}{p},
-$$
-
-resulting in a speedup of
-$$
-\frac{T_1}{T_p} = \frac{W}{(1-f)W+fW/p}
-$$
-
-
-As \( p \) goes to infinity, \( fW/p \) goes to zero, and the maximum speedup is
-$$
-\frac{1}{1-f},
-$$
-
-meaning that if
-if \( f = 0.99 \) (all but \( 1\% \) parallelizable), the maximum speedup
-is \( 1/(1-.99)=100 \)!
+
Automatic vectorization and vectorization inhibitors, more data dependencies
+
+
On the other hand, a so-called
+write-after-read statement can be vectorized. The following code
+
is an example of flow dependency that can be vectorized since no iteration with a higher value of \( i \)
+can complete before an iteration with a lower value of \( i \). However, such code leads to problems with parallelization.
+
-If any non-parallel code slips into the
-application, the parallel
-performance is limited.
-
-
-In many simulations, however, the fraction of non-parallelizable work
-is \( 10^{-6} \) or less due to large arrays or objects that are perfectly parallelizable.
-
-
+
Automatic vectorization and vectorization inhibitors, memory stride
+
+
For C++ programmers it is also worth keeping in mind that an array notation is preferred to the more compact use of pointers to access array elements. The compiler can often not tell if it is safe to vectorize the code.
+
+
When dealing with arrays, you should also avoid memory stride, since this slows down considerably vectorization. When you access array element, write for example the inner loop to vectorize using unit stride, that is, access successively the next array element in memory, as shown here
Cache memory contains a copy of the main memory data
+
Cache is faster but consumes more space and power. It is normally assumed to be much faster than main memory
+
Registers contain working data only
-
Distributed memory is the dominant hardware configuration. There is a large diversity in these machines, from MPP (massively parallel processing) systems to clusters of off-the-shelf PCs, which are very cost-effective.
-
Message-passing is a mature programming paradigm and widely accepted. It often provides an efficient match to the hardware. It is primarily used for the distributed memory systems, but can also be used on shared memory systems.
-
Modern nodes have nowadays several cores, which makes it interesting to use both shared memory (the given node) and distributed memory (several nodes with communication). This leads often to codes which use both MPI and OpenMP.
+
Modern CPUs perform most or all operations only on data in register
+
Multiple Cache memories contain a copy of the main memory data
+
+
Cache items accessed by their address in main memory
+
L1 cache is the fastest but has the least capacity
+
L2, L3 provide intermediate performance/size tradeoffs
+
+
+
Loads and stores to memory can be as important as floating point operations when we measure performance.
-Our lectures will focus on both MPI and OpenMP.
-
-
Develop codes locally, run with some few processes and test your codes. Do benchmarking, timing and so forth on local nodes, for example your laptop or PC.
-
When you are convinced that your codes run correctly, you can start your production runs on available supercomputers.
-
-
-
-
-
-
+
Problems with measuring time
+
+
Timers are not infinitely accurate
+
All clocks have a granularity, the minimum time that they can measure
+
The error in a time measurement, even if everything is perfect, may be the size of this granularity (sometimes called a clock tick)
+
Always know what your clock granularity is
+
Ensure that your measurement is for a long enough duration (say 100 times the tick)
-To install MPI is rather easy on hardware running unix/linux as operating systems, follow simply the instructions from the OpenMPI website. See also subsequent slides.
-When you have made sure you have installed MPI on your PC/laptop,
-
-
-
Compile with mpicxx/mpic++ or mpif90
-
-
-
-
-
-
# Compile and link
- mpic++-O3 -o nameofprog.x nameofprog.cpp
- # run code with for example 8 processes using mpirun/mpiexec
- mpiexec -n 8 ./nameofprog.x
-
-
-
-
-
-
-
+
Problems with cold start
+
+
What happens when the code is executed? The assumption is that the code is ready to
+execute. But
+
+
+
Code may still be on disk, and not even read into memory.
+
Data may be in slow memory rather than fast (which may be wrong or right for what you are measuring)
+
Multiple tests often necessary to ensure that cold start effects are not present
+
Special effort often required to ensure data in the intended part of the memory hierarchy.
Can I do it on my own PC/laptop? OpenMP installation
-
-
-
-If you wish to install MPI and OpenMP
-on your laptop/PC, we recommend the following:
-
-
-
For OpenMP, the compile option -fopenmp is included automatically in recent versions of the C++ compiler and Fortran compilers. For users of different Linux distributions, simply use the available C++ or Fortran compilers and add the above compiler instructions, see also code examples below.
-
For OS X users however, install libomp
-
-
-
-
-
-
brew install libomp
-
-
-and compile and link as
-
-
-
-
c++-o <name executable><name program.cpp>-lomp
-
-
-
-
-
-
-
+
Problems with smart compilers
+
+
+
If the result of the computation is not used, the compiler may eliminate the code
+
Performance will look impossibly fantastic
+
Even worse, eliminate some of the code so the performance looks plausible