pythonpetsc

Convert a list of Vector to a matrix with petsc4py

I have a list of vectors that I want to gather in a single matrix Z. Here is the code I tried to use :

import sys, slepc4py
slepc4py.init(sys.argv)
from petsc4py import PETSc
from slepc4py import SLEPc

n = 5
d = 10

l = []   # creation of the list of vectors
for i in range(d):
    v = PETSc.Vec().create()
    v.setSizes(n)
    v.setFromOptions()
    v.set(i)
    l.append(v)

Z = PETSc.Mat().create()
Z.setSizes([n, d])
Z.setFromOptions()
Z.setUp()
for i, v in enumerate(l):
    # Z.setValuesBlocked(range(n), i, v)
    Z.setValues(range(n), i, v)
Z.assemble()
Z.view()    # to display the matrix in the terminal

I tried setValues and setValuesBlocked. When I run this script on one core, that works well. The result printed is :

row 0: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 1: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 2: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 3: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 4: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.)

but in parallel I have issues :

Traceback (most recent call last):
  File "t.py", line 26, in <module>
    Z.setValues(range(n), i, v)
  File "PETSc/Mat.pyx", line 885, in petsc4py.PETSc.Mat.setValues
  File "PETSc/petscmat.pxi", line 826, in petsc4py.PETSc.matsetvalues
ValueError: incompatible array sizes: ni=5, nj=1, nv=3
Traceback (most recent call last):
  File "t.py", line 26, in <module>
    Z.setValues(range(n), i, v)
  File "PETSc/Mat.pyx", line 885, in petsc4py.PETSc.Mat.setValues
  File "PETSc/petscmat.pxi", line 826, in petsc4py.PETSc.matsetvalues
ValueError: incompatible array sizes: ni=5, nj=1, nv=2

I think it's because PetSc waits for the local indices in the function setValues, but I don't know how to access it (local_size gives the local size of the matrix, if I put Z.local_size[0] in setValues, the result in parallel is

Mat Object: 2 MPI processes
  type: mpiaij
row 0: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 1: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 2: (0, 0.)  (1, 1.)  (2, 2.)  (3, 3.)  (4, 4.)  (5, 5.)  (6, 6.)  (7, 7.)  (8, 8.)  (9, 9.) 
row 3:
row 4:

)

Which function or what argument should I use to fulfill this ?

EDIT :

I tried a more tedious way to fill the matrix Z, using the following double loop :

for i, v in enumerate(l):
    for j in range(n):
        Z.setValue(j, i, v[j])

once again in sequential, that works fine, but if I run it over 2 cores here is the result, totally wrong !

Mat Object: 2 MPI processes
  type: mpiaij
row 0: (0, 6.93107e-310)  (1, 6.93107e-310)  (2, 6.93107e-310)  (3, 0.)  (4, 6.93107e-310)  (5, 6.93107e-310)  (6, 6.93107e-310)  (7, 6.93107e-310)  (8, 6.93107e-310)  (9, 6.93107e-310) 
row 1: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, -1.)  (7, -1.)  (8, 0.)  (9, 0.) 
row 2: (0, 1.63042e-322)  (1, 1.63042e-322)  (2, 1.63042e-322)  (3, 1.63042e-322)  (4, 1.63042e-322)  (5, 1.63042e-322)  (6, 1.63042e-322)  (7, 1.63042e-322)  (8, 1.63042e-322)  (9, 1.63042e-322) 
row 3: (0, 0.)  (1, 0.)  (2, 0.)  (3, 0.)  (4, 0.)  (5, 0.)  (6, 0.)  (7, 0.)  (8, 0.)  (9, 0.) 
row 4: (0, 0.)  (1, -1.)  (2, 0.)  (3, 0.)  (4, -1.)  (5, 4.24399e-314)  (6, 0.)  (7, -1.)  (8, -1.)  (9, 0.)
Solution

  • Thank to the help of the team developping PetSc, I managed to do it. Here is the code :

    from petsc4py import PETSc

n = 10
d = 5

l = []   # creation of the list of vectors
for i in range(d):
    v = PETSc.Vec().create(comm=PETSc.COMM_WORLD)
    v.setSizes(n)
    v.setFromOptions()
    v.setUp()
    v.set(i)
    l.append(v)

Z = PETSc.Mat().createDense((n,d),comm=PETSc.COMM_WORLD)
Z.setFromOptions()
Z.setUp()
for i, v in enumerate(l):
    r = Z.getDenseColumnVec(i)
    v.copy(r)
    Z.restoreDenseColumnVec(i)
Z.assemble()
Z.view()

    We need to use a dense matrix, and the function getDenseColumnVec allows to access to a column of the matrix, and modify it at will.

    That works both in sequential and parallel.

    Mat Object: 6 MPI processes
  type: mpidense
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00 
0.0000000000000000e+00 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00