PETSc from python

[Marek Wojciechowski]

On Wed, 16 Aug 2006 11:07:42 -0000, Matthew Knepley <knepley at gmail.com>  
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 at 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