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