attachment:laplace.py of PerformancePython - SciPy wiki dump (original) (raw)

1 2 15 18 19 20 21 26 31 34 35 36 42 44 46 48 50 51 52 56 59 60 61 62 63 67 69 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 89 90 92 98 101 102 104 105 106 107 109 110 112 113 114 115 116 117 118 119 121 124 125 126 127 128 129 133 137 139 141 142 143 144 149 151 153 155 157 158 159 160 163 164 165 166 167 168 169 170 171 172 173 174 175 176 179 180 181 182 184 185 186 187 188 189 190 191 192 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 216 217 218 220 221 222 223 224 225 226 227 228 230 231 232 233 234 236 237 239 241 242 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 262 263 264 265 266 267 269 271 272 273 274 275 277 282 287 288 289 290 291 292 293 294 295 296 297 298 299 303 304 305 306 307 308 312 313 314 315 316 317 319 320 321 322 323 325 326 327 328 329 330 331 332 333 334 335 336 341 3 """ 4 This script compares different ways of implementing an iterative 5 procedure to solve Laplace's equation. These provide a general 6 guideline to using Python for high-performance computing and also 7 provide a simple means to compare the computational time taken by the 8 different approaches. The script compares functions implemented in 9 pure Python, Numeric, weave.blitz, weave.inline, fortran (via f2py) 10 and Pyrex. The function main(), additionally accelerates the pure 11 Python version using Psyco and provides some numbers on how well that 12 works. To compare all the options you need to have Numeric, weave, 13 f2py, Pyrex and Psyco installed. If Psyco is not installed the script 14 will print a warning but will perform all other tests. 16 The fortran and pyrex modules are compiled using the setup.py script 17 that is provided with this file. You can build them like so: python setup.py build_ext --inplace 22 Author: Prabhu Ramachandran 23 License: BSD 24 Last modified: Sep. 18, 2004 25 """ 27 import numpy 28 from scipy import weave 29 from scipy.weave import converters 30 import sys, time 32 msg = """************************************************** 33 Please build the fortran and Pyrex modules like so: python setup.py build_ext --inplace 37 You will require f2py and Pyrex. 38 ************************************************** 39 """ 40 build = 0 41 try: import flaplace 43 except ImportError: build = 1 45 try: import pyx_lap 47 except ImportError: build = 1 49 if build: print msg 53 class Grid: 54
55 """A simple grid class that stores the details and solution of the computational grid.""" 57
58 def init(self, nx=10, ny=10, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0): self.xmin, self.xmax, self.ymin, self.ymax = xmin, xmax, ymin, ymax self.dx = float(xmax-xmin)/(nx-1) self.dy = float(ymax-ymin)/(ny-1) self.u = numpy.zeros((nx, ny), 'd') 64
65 self.old_u = self.u.copy()
66 def setBC(self, l, r, b, t):
68 """Sets the boundary condition given the left, right, bottom and top values (or arrays)"""
70 self.u[0, :] = l self.u[-1, :] = r self.u[:, 0] = b self.u[:,-1] = t self.old_u = self.u.copy() def setBCFunc(self, func): """Sets the BC given a function of two variables.""" xmin, ymin = self.xmin, self.ymin xmax, ymax = self.xmax, self.ymax x = numpy.arange(xmin, xmax + self.dx0.5, self.dx) y = numpy.arange(ymin, ymax + self.dy0.5, self.dy) self.u[0 ,:] = func(xmin,y) self.u[-1,:] = func(xmax,y) self.u[:, 0] = func(x,ymin) self.u[:,-1] = func(x,ymax) def computeError(self):
88 """Computes absolute error using an L2 norm for the solution. This requires that self.u and self.old_u must be appropriately setup."""
91 v = (self.u - self.old_u).flat return numpy.sqrt(numpy.dot(v,v)) 93
94 95 class LaplaceSolver: 96
97 """A simple Laplacian solver that can use different schemes to solve the problem.""" 99
100 def init(self, grid, stepper='numeric'): self.grid = grid self.setTimeStepper(stepper) 103 def slowTimeStep(self, dt=0.0): """Takes a time step using straight forward Python loops.""" g = self.grid nx, ny = g.u.shape
108 dx2, dy2 = g.dx2, g.dy2 dnr_inv = 0.5/(dx2 + dy2) u = g.u 111 err = 0.0 for i in range(1, nx-1): for j in range(1, ny-1): tmp = u[i,j] u[i,j] = ((u[i-1, j] + u[i+1, j])*dy2 + (u[i, j-1] + u[i, j+1])*dx2)dnr_inv diff = u[i,j] - tmp err += diffdiff 120 return numpy.sqrt(err) 122
123 def numericTimeStep(self, dt=0.0): """Takes a time step using a numeric expressions.""" g = self.grid dx2, dy2 = g.dx2, g.dy2 dnr_inv = 0.5/(dx2 + dy2) u = g.u g.old_u = u.copy() 130 131
132 u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 + (u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv 134
135 return g.computeError() 136 def blitzTimeStep(self, dt=0.0):
138 """Takes a time step using a numeric expression that has been blitzed using weave."""
140 g = self.grid dx2, dy2 = g.dx2, g.dy2 dnr_inv = 0.5/(dx2 + dy2) u = g.u g.old_u = u.copy() 145 146
147 expr = "u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 + "
148 "(u[1:-1,0:-2] + u[1:-1, 2:])dx2)dnr_inv" weave.blitz(expr, check_size=0) 150 return g.computeError() 152 def inlineTimeStep(self, dt=0.0):
154 """Takes a time step using inlined C code -- this version uses blitz arrays."""
156 g = self.grid nx, ny = g.u.shape dx2, dy2 = g.dx2, g.dy2 dnr_inv = 0.5/(dx2 + dy2) u = g.u 161
162 code = """ #line 120 "laplace.py" double tmp, err, diff; err = 0.0; for (int i=1; i<nx-1; ++i) { for (int j=1; j<ny-1; ++j) { tmp = u(i,j); u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2 + (u(i,j-1) + u(i,j+1))*dx2)*dnr_inv; diff = u(i,j) - tmp; err += diff*diff; } } return_val = sqrt(err); """ 177
178 err = weave.inline(code, ['u', 'dx2', 'dy2', 'dnr_inv', 'nx','ny'], type_converters = converters.blitz, compiler = 'gcc') return err 183 def fastInlineTimeStep(self, dt=0.0): """Takes a time step using inlined C code -- this version is faster, dirtier and manipulates the numeric array in C. This code was contributed by Eric Jones. """ g = self.grid nx, ny = g.u.shape dx2, dy2 = g.dx2, g.dy**2 dnr_inv = 0.5/(dx2 + dy2) u = g.u 193
194 code = """ #line 151 "laplace.py" double tmp, err, diff; double *uc, *uu, *ud, *ul, *ur; err = 0.0; for (int i=1; i<nx-1; ++i) { uc = u+i*ny+1; ur = u+i*ny+2; ul = u+i*ny; ud = u+(i+1)*ny+1; uu = u+(i-1)*ny+1; for (int j=1; j<ny-1; ++j) { tmp = *uc; *uc = ((*ul + *ur)*dy2 + (*uu + *ud)*dx2)*dnr_inv; diff = *uc - tmp; err += diff*diff; uc++;ur++;ul++;ud++;uu++; } } return_val = sqrt(err); """ 214
215 err = weave.inline(code, ['u', 'dx2', 'dy2', 'dnr_inv', 'nx','ny'], compiler='gcc') return err 219 def fortranTimeStep(self, dt=0.0): """Takes a time step using a simple fortran module that implements the loop in fortran. Use f2py to compile flaplace.f like so: f2py -c flaplace.c -m flaplace. You need the latest f2py version for this to work. This Fortran example was contributed by Pearu Peterson. """ g = self.grid g.u, err = flaplace.timestep(g.u, g.dx, g.dy) return err 229 def pyrexTimeStep(self, dt=0.0): """Takes a time step using a function written in Pyrex. Use the given setup.py to build the extension using the command python setup.py build_ext --inplace. You will need Pyrex installed to run this."""
235 g = self.grid err = pyx_lap.pyrexTimeStep(g.u, g.dx, g.dy) return err 238 def setTimeStepper(self, stepper='numeric'):
240 """Sets the time step scheme to be used while solving given a string which should be one of ['slow', 'numeric', 'blitz', 'inline', 'fastinline', 'fortran']."""
243 if stepper == 'slow': self.timeStep = self.slowTimeStep elif stepper == 'numeric': self.timeStep = self.numericTimeStep elif stepper == 'blitz': self.timeStep = self.blitzTimeStep elif stepper == 'inline': self.timeStep = self.inlineTimeStep elif stepper.lower() == 'fastinline': self.timeStep = self.fastInlineTimeStep elif stepper == 'fortran': self.timeStep = self.fortranTimeStep elif stepper == 'pyrex': self.timeStep = self.pyrexTimeStep else: self.timeStep = self.numericTimeStep
259
260 def solve(self, n_iter=0, eps=1.0e-16):
261 """Solves the equation given an error precision -- eps. If n_iter=0 the solving is stopped only on the eps condition. If n_iter is finite then solution stops in that many iterations or when the error is less than eps whichever is earlier. Returns the error if the loop breaks on the n_iter condition and returns the iterations if the loop breaks on the error condition."""
268 err = self.timeStep() count = 1 270 while err > eps: if n_iter and count >= n_iter: return err err = self.timeStep() count = count + 1 276 return count 278 279 280 def BC(x, y):
281 """Used to set the boundary condition for the grid of points. Change this as you feel fit."""
283 return (x2 - y2) 284 285 286 def test(nmin=5, nmax=30, dn=5, eps=1.0e-16, n_iter=0, stepper='numeric'): iters = [] n_grd = numpy.arange(nmin, nmax, dn) times = [] for i in n_grd: g = Grid(nx=i, ny=i) g.setBCFunc(BC) s = LaplaceSolver(g, stepper) t1 = time.clock() iters.append(s.solve(n_iter=n_iter, eps=eps)) dt = time.clock() - t1 times.append(dt) print "Solution for nx = ny = %d, took %f seconds"%(i, dt) return (n_grd2, iters, times) 300 301 302 def time_test(nx=500, ny=500, eps=1.0e-16, n_iter=100, stepper='numeric'): g = Grid(nx, ny) g.setBCFunc(BC) s = LaplaceSolver(g, stepper) t = time.clock() s.solve(n_iter=n_iter, eps=eps) return time.clock() - t 309
310 311 def main(n=500, n_iter=100): print "Doing %d iterations on a %dx%d grid"%(n_iter, n, n) for i in ['numeric', 'blitz', 'inline', 'fastinline', 'fortran', 'pyrex']: print i, sys.stdout.flush() print "took", time_test(n, n, stepper=i, n_iter=n_iter), "seconds" 318 print "slow (1 iteration)", sys.stdout.flush() s = time_test(n, n, stepper='slow', n_iter=1) print "took", s, "seconds" print "%d iterations should take about %f seconds"%(n_iter, sn_iter) 324 try: import psyco except ImportError: print "You don't have Psyco installed!" else: psyco.bind(LaplaceSolver) psyco.bind(Grid) print "slow with Psyco (1 iteration)", sys.stdout.flush() s = time_test(n, n, stepper='slow', n_iter=1) print "took", s, "seconds" print "%d iterations should take about %f seconds"%
337 (n_iter, sn_iter) 338 339 340 if name == "main": main()

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the

link, since this is subject to change and can break easily.