[petsc-users] How to create a local to global mapping and construct matrix correctly

[Barry Smith]

   "Gives wrong answers" is not very informative. What answer do you expect and what answer do you get?

  Note that each process is looping over mSizes?

for i in range(len(mSizes)):
    for j in range(len(mSizes)):

  Is this what you want? It doesn't seem likely that you want all processes to generate all information in the matrix. Each process should be doing a subset of the generation.

   Barry

> On Sep 16, 2016, at 11:03 AM, Ji Zhang <gotofd at gmail.com> wrote:
> 
> Dear all, 
> 
> I have a number of small 'mpidense' matrices mij, and I want to construct them to a big 'mpidense' matrix M like this:
>      [  m11  m12  m13  ]
> M =  |  m21  m22  m23  |   ,
>      [  m31  m32  m33  ]
> 
> And a short demo is below. I'm using python, but their grammar are similar. 
> import numpy as np
> from petsc4py import PETSc
> import sys, petsc4py
> 
> 
> petsc4py.init(sys.argv)
> mSizes = (2, 2)
> mij = []
> 
> # create sub-matrices mij
> for i in range(len(mSizes)):
>     for j in range(len(mSizes)):
>         temp_m = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
>         temp_m.setSizes(((None, mSizes[i]), (None, mSizes[j])))
>         temp_m.setType('mpidense')
>         temp_m.setFromOptions()
>         temp_m.setUp()
>         temp_m[:, :] = np.random.random_sample((mSizes[i], mSizes[j]))
>         temp_m.assemble()
>         temp_m.view()
>         mij.append(temp_m)
> 
> # Now we have four sub-matrices. I would like to construct them into a big matrix M.
> M = PETSc.Mat().create(comm=PETSc.COMM_WORLD)
> M.setSizes(((None, np.sum(mSizes)), (None, np.sum(mSizes))))
> M.setType('mpidense')
> M.setFromOptions()
> M.setUp()
> mLocations = np.insert(np.cumsum(mSizes), 0, 0)    # mLocations = [0, mSizes]
> for i in range(len(mSizes)):
>     for j in range(len(mSizes)):
>         temp_m = mij[i*len(mSizes)+j].getDenseArray()
>         for k in range(temp_m.shape[0]):
>             M.setValues(mLocations[i]+k, np.arange(mLocations[j],mLocations[j+1],dtype='int32'), temp_m[k, :])
> M.assemble()
> M.view()
> The code works well in a single cup, but give wrong answer for 2 and more cores. 
> 
> Thanks. 
> 2016-09-17
> Best, 
> Regards,
> Zhang Ji 
> Beijing Computational Science Research Center 
> E-mail: gotofd at gmail.com
> 
>