如何构造petsc矩阵？

Question

我正在使用petsc4py，现在遇到了一些问题。我有一些小的petsc密集矩阵mij，我想把它们构造成一个像这样的大矩阵M：

     [  m11  m12  m13  ]
M =  |  m21  m22  m23  |   ,
     [  m31  m32  m33  ]

下面显示了一段mcve代码，我使用的是PETSc的python包装，但是它们的语法是相似的。

import numpy as np
from petsc4py import PETSc

mSizes = (5, 8, 6)
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()
        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)):
        M[mLocations[i]:mLocations[i+1], mLocations[j]:mLocations[j+1]] = \ 
            mij[i*len(mSizes)+j][:, :]
M.assemble()

它们会报告这样的错误消息：

Traceback (most recent call last):
  File "tryMatConstuct.py", line 29, in <module>
    mij[i*len(mSizes)+j][:, :]
  File "PETSc/Mat.pyx", line 227, in petsc4py.PETSc.Mat.__getitem__ (src/petsc4py.PETSc.c:110477)
  File "PETSc/petscmat.pxi", line 997, in petsc4py.PETSc.mat_getitem (src/petsc4py.PETSc.c:30481)
  File "PETSc/petscmat.pxi", line 917, in petsc4py.PETSc.matgetvalues (src/petsc4py.PETSc.c:29242)
petsc4py.PETSc.Error: error code 56
[1] MatGetValues() line 1818 in /home/zhangji/PycharmProjects/petsc-petsc-31a1859eaff6/src/mat/interface/matrix.c
[1] MatGetValues_MPIDense() line 154 in /home/zhangji/PycharmProjects/petsc-petsc-31a1859eaff6/src/mat/impls/dense/mpi/mpidense.c
[1] No support for this operation for this object type
[1] Only local values currently supported

Answer 1

我必须使用二维数组而不是petsc矩阵来缩短子矩阵mij。

在petsc4py中，函数是PETSc.Mat().getDenseArray()，即

temp_m = mij[i*len(mSizes)+j].getDenseArray()