如何构造petsc矩阵？答案

【问题标题】：How to construct petsc matrices?如何构造petsc矩阵？
【发布时间】：2016-09-12 22:46:59
【问题描述】：

我正在使用 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

【问题讨论】：

格式 MATMPIDENSE 对应于切片矩阵：每个进程将整个矩阵的一组给定行存储为 MATSEQDENSE 矩阵。函数MatGetValues() 只允许访问本地数据：这就是打印错误消息的原因。同样MatDenseGetArray_MPIDense() 允许访问本地数据。见mpidense.c。调度数据所需的 MPI 通信留给您！我们需要mcve 为您提供答案！

标签： petsc

【解决方案1】：

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

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

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

【讨论】：