1try: range = xrange 2except NameError: pass 3 4import sys, petsc4py 5petsc4py.init(sys.argv) 6 7from petsc4py import PETSc 8 9m, n = 16, 32 10 11A = PETSc.Mat().create(PETSc.COMM_WORLD) 12A.setSizes([m*n, m*n]) 13A.setFromOptions() 14A.setUp() 15Istart, Iend = A.getOwnershipRange() 16for I in range(Istart, Iend): 17 A[I,I] = 4 18 i = I//n 19 if i>0 : J = I-n; A[I,J] = -1 20 if i<m-1: J = I+n; A[I,J] = -1 21 j = I-i*n 22 if j>0 : J = I-1; A[I,J] = -1 23 if j<n-1: J = I+1; A[I,J] = -1 24A.assemblyBegin() 25A.assemblyEnd() 26 27x, y = A.createVecs() 28x.set(1) 29A.mult(x,y) 30 31# save 32viewer = PETSc.Viewer().createBinary('matrix-A.dat', 'w') 33viewer(A) 34viewer = PETSc.Viewer().createBinary('vector-x.dat', 'w') 35viewer(x) 36viewer = PETSc.Viewer().createBinary('vector-y.dat', 'w') 37viewer(y) 38 39# load 40viewer = PETSc.Viewer().createBinary('matrix-A.dat', 'r') 41B = PETSc.Mat().load(viewer) 42viewer = PETSc.Viewer().createBinary('vector-x.dat', 'r') 43u = PETSc.Vec().load(viewer) 44viewer = PETSc.Viewer().createBinary('vector-y.dat', 'r') 45v = PETSc.Vec().load(viewer) 46 47# check 48assert B.equal(A) 49assert x.equal(u) 50assert y.equal(v) 51