[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: PETSc from python

On Wed, 16 Aug 2006 11:07:42 -0000, Matthew Knepley <knepley@xxxxxxxxx> wrote:

This is very interesting. I also know Lisandro who wrote bwpython. I just had a user have the opposite experience with these packages, but I guess that
is why we have multiple packages. I am very happy you got this going.

Yes. I got this going. And I'm quite satisfied. To prove it, that's an examplary session:


mwojc@evo ~ $ mpirun -np 2 `which bwpython` `which ipython` -nobanner

In [1]: from petscinit import * # This imports for me PETSc extensions, initializes and creates stdviewer

In [2]: from mpi4py import MPI

In [3]:

In [3]: v=Vec()

In [4]: size = (MPI.rank + 1) * 2

In [5]: totsize = MPI.COMM_WORLD.Allreduce(size)

In [6]: print 'Process [%d]: size=%d' %(MPI.rank, size)
Process [0]: size=2
Process [1]: size=4

In [7]: print 'Process [%d]: totsize=%d' %(MPI.rank, totsize)
Process [0]: totsize=6
Process [1]: totsize=6

In [8]: v.setSizes(size, PETSC_DECIDE)

In [9]: v.setFromOptions()

In [10]: if MPI.rank == 0: v.setValues(xrange(totsize), [1]*totsize, INSERT_VALUES)
....:

In [11]: v.assemblyBegin()

In [12]: v.assemblyEnd()

In [13]: v.view(stdviewer)
Process [0]
0: 1
1: 1
Process [1]
2: 1
3: 1
4: 1
5: 1

In [14]:

In [14]: A=Mat()

In [15]: A.setSizes(size, size, PETSC_DECIDE, PETSC_DECIDE)

In [16]: A.setFromOptions()

In [17]: Range = A.getOwnershipRange()

In [18]: print Range
(0, 2)
(2, 6)

In [19]: rows=xrange(Range[0], Range[1])

In [20]: cols=xrange(totsize)

In [21]: import random

In [22]: values=[random.uniform(-1,1) for i in xrange(size*totsize)]

In [23]: A.setValues(rows, cols, values, INSERT_VALUES)

In [24]: A.assemblyBegin(MAT_FINAL_ASSEMBLY)

In [25]: A.assemblyEnd(MAT_FINAL_ASSEMBLY)

In [26]: A.view(stdviewer)
row 0: (0, 0.630031) (1, 0.673476) (2, -0.734869) (3, 0.105727) (4, 0.538428) (5, 0.12576)
row 1: (0, -0.857206) (1, -0.0761736) (2, -0.143492) (3, -0.938166) (4, 0.41378) (5, -0.210328)
row 2: (0, 0.50173) (1, -0.214067) (2, 0.59921) (3, 0.848044) (4, -0.819785) (5, -0.436404)
row 3: (0, 0.30529) (1, 0.968145) (2, 0.377928) (3, -0.656585) (4, 0.882831) (5, 0.850657)
row 4: (0, -0.304465) (1, 0.496273) (2, -0.277161) (3, -0.81206) (4, 0.63498) (5, 0.58123)
row 5: (0, 0.538759) (1, -0.654964) (2, -0.256906) (3, -0.335948) (4, 0.748973) (5, 0.813876)

In [27]:

In [27]: b = v.duplicate() #right hand vector

In [28]: A.mult(v, b)

In [29]: x=b.duplicate() #unknowns

In [30]:

In [30]: solver = KSP()

In [31]: solver.setOperators(A,A,DIFFERENT_NONZERO_PATTERN)

In [32]: solver.setFromOptions()

In [33]: solver.solve(b,x)

In [34]: x.view(stdviewer)    # IS SOLUTION CORRECT? (SHOULD BE ONES?)
Process [0]
0: 1
1: 1
Process [1]
2: 1
3: 1
4: 1
5: 1

In [35]: #AND SO ON...


Isn't it nice?
But I have also a question. Suppose my matrix A is "dense" and I would like to get local arrays:


In [36]: A.setType("dense")

In [37]: A.setValues(rows, cols, values, INSERT_VALUES)

In [38]: A.assemblyBegin(MAT_FINAL_ASSEMBLY)

In [39]: A.assemblyEnd(MAT_FINAL_ASSEMBLY)

In [40]: A.view(stdviewer)
6.3003e-01 6.7348e-01 -7.3487e-01 1.0573e-01 5.3843e-01 1.2576e-01
-8.5721e-01 -7.6174e-02 -1.4349e-01 -9.3817e-01 4.1378e-01 -2.1033e-01
5.0173e-01 -2.1407e-01 5.9921e-01 8.4804e-01 -8.1978e-01 -4.3640e-01
3.0529e-01 9.6815e-01 3.7793e-01 -6.5659e-01 8.8283e-01 8.5066e-01
-3.0446e-01 4.9627e-01 -2.7716e-01 -8.1206e-01 6.3498e-01 5.8123e-01
5.3876e-01 -6.5496e-01 -2.5691e-01 -3.3595e-01 7.4897e-01 8.1388e-01

In [41]: print A.getArray()
array([0.63003134676816508, -0.85720600630413402, 0.67347596051594882, -0.076173622638120664], 'd')
array([0.5017303230541359, 0.30529049321568058, -0.3044646591900968, 0.53875941791034943, -0.21406665507357259, 0.96814544907801547, 0.49627269833290644, -0.65496405056370244, 0.59920997972823065, 0.3779276741558788, -0.2771608364332232, -0.25690564212998068, 0.8480435858237616, -0.65658527037718994, -0.81206017169215183, -0.33594803008233565], 'd')


Why the obtained arrays does not represent what is stored on each process. I thought there should be 2 first whole rows in the first array and the rest in the second or I'm missing something.


I also observed that PETSc objects are not picklable. Is there any good reason for that?


Greetings

--
Marek Wojciechowski