高效构造有限元/FVM矩阵

Question

这是FEM/FVM方程系统的一个典型用例，因此可能会引起更广泛的兴趣。从三角形网格

​

​

我想创建一个scipy.sparse.csr_matrix。矩阵行/列表示网格节点上的值。该矩阵在主对角线和两个节点通过边缘连接的地方都有条目。

下面是一个MWE，它首先构建节点->边缘->单元关系，然后构建矩阵：

import numpy
import meshzoo
from scipy import sparse

nx = 1600
ny = 1000
verts, cells = meshzoo.rectangle(0.0, 1.61, 0.0, 1.0, nx, ny)

n = len(verts)

nds = cells.T
nodes_edge_cells = numpy.stack([nds[[1, 2]], nds[[2, 0]],nds[[0, 1]]], axis=1)

# assign values to each edge (per cell)
alpha = numpy.random.rand(3, len(cells))
vals = numpy.array([
    [alpha**2, -alpha],
    [-alpha, alpha**2],
    ])

# Build I, J, V entries for COO matrix
I = []
J = []
V = []
#
V.append(vals[0][0])
V.append(vals[0][1])
V.append(vals[1][0])
V.append(vals[1][1])
#
I.append(nodes_edge_cells[0])
I.append(nodes_edge_cells[0])
I.append(nodes_edge_cells[1])
I.append(nodes_edge_cells[1])
#
J.append(nodes_edge_cells[0])
J.append(nodes_edge_cells[1])
J.append(nodes_edge_cells[0])
J.append(nodes_edge_cells[1])
# Create suitable data for coo_matrix
I = numpy.concatenate(I).flat
J = numpy.concatenate(J).flat
V = numpy.concatenate(V).flat

matrix = sparse.coo_matrix((V, (I, J)), shape=(n, n))
matrix = matrix.tocsr()

使用

python -m cProfile -o profile.prof main.py
snakeviz profile.prof

您可以创建和查看上述的概要文件：

​

​

方法tocsr()在这里获得了运行时的最大份额，但是当构建alpha更复杂时，这也是正确的。因此，我正在寻找加速这一进程的方法。

我已经发现了：

• 由于数据的结构，矩阵对角线上的值可以预先求和，即， V.append(vals0，0，0+ vals1，1，2) I.append(nodes_edge_cells0，0) # == nodes_edge_cells1，2 J.append(nodes_edge_cells0，0) # == nodes_edge_cells1，2 这使得I、J、V更短，从而加快了tocsr的速度。
• 现在，边缘是“每个细胞”。我可以使用numpy.unique识别彼此相等的边，有效地节省了大约一半的I，J，V。然而，我发现这也需要一些时间。(不足为奇)

我的另一个想法是，如果存在一个主对角线单独保存的V -like数据结构，则可以用简单的numpy.add.at替换对角线numpy.add.at、I、J。我知道在其他软件包中也存在这种情况，但在find中找不到它。对，是这样?

也许有一种明智的方法直接构建企业社会责任？

Answer 1

因此，到头来，这是首席运营官和CSR的sum_duplicates之间的区别(就像@hpaulj怀疑的那样)。多亏这里的每一个人的努力(特别是@保罗-帕泽)，公关正在为tocsr提供一个巨大的加速。

SciPy的tocsr在(I, J)上做了一个lexsort，所以它可以帮助组织索引，这样(I, J)就会得到很好的排序。

对于nx=4，上面示例中的ny=2，I和J是

[1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7]
[1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7 1 6 3 5 2 7 5 5 7 4 5 6 0 2 2 0 1 2 5 5 7 4 5 6 0 2 2 0 1 2 1 6 3 5 2 7]

首先对cells的每一行进行排序，然后按照第一列对行进行排序，如下

cells = numpy.sort(cells, axis=1)
cells = cells[cells[:, 0].argsort()]

产生

[1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6]
[1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6 1 4 2 5 3 6 5 5 5 6 7 7 0 0 1 2 2 2 5 5 5 6 7 7 0 0 1 2 2 2 1 4 2 5 3 6]

对于原始帖子中的数字，排序将运行时从大约40秒减少到8秒。

如果首先对节点进行更适当的编号，也许可以实现更好的排序。我在想卡蒂尔-麦基和朋友。

Answer 2

我将尝试直接创建csr结构，特别是如果您使用np.unique，因为这会给您排序的键，这是工作完成的一半。

我假设您所处的位置是这样的：i, j按字典顺序排序，重叠v在np.unique的可选inverse输出上使用np.add.at相加。

那么v和j已经采用了csr格式。剩下要做的就是创建indptr，您只需通过np.searchsorted(i, np.arange(M+1))获得它，其中M是列的长度。您可以直接将它们传递给sparse.csr_matrix构造函数。

好的，让代码说：

import numpy as np
from scipy import sparse
from timeit import timeit


def tocsr(I, J, E, N):
    n = len(I)
    K = np.empty((n,), dtype=np.int64)
    K.view(np.int32).reshape(n, 2).T[...] = J, I  
    S = np.argsort(K)
    KS = K[S]
    steps = np.flatnonzero(np.r_[1, np.diff(KS)])
    ED = np.add.reduceat(E[S], steps)
    JD, ID = KS[steps].view(np.int32).reshape(-1, 2).T
    ID = np.searchsorted(ID, np.arange(N+1))
    return sparse.csr_matrix((ED, np.array(JD, dtype=int), ID), (N, N))

def viacoo(I, J, E, N):
    return sparse.coo_matrix((E, (I, J)), (N, N)).tocsr()


#testing and timing

# correctness
N = 1000
A = np.random.random((N, N)) < 0.001
I, J = np.where(A)

E = np.random.random((2, len(I)))
D = np.zeros((2,) + A.shape)
D[:, I, J] = E
D2 = tocsr(np.r_[I, I], np.r_[J, J], E.ravel(), N).A

print('correct:', np.allclose(D.sum(axis=0), D2))

# speed
N = 100000
K = 10

I, J = np.random.randint(0, N, (2, K*N))
E = np.random.random((2 * len(I),))
I, J, E = np.r_[I, I, J, J], np.r_[J, J, I, I], np.r_[E, E]

print('N:', N, ' --  nnz (with duplicates):', len(E))
print('direct: ', timeit('f(a,b,c,d)', number=10, globals={'f': tocsr, 'a': I, 'b': J, 'c': E, 'd': N}), 'secs for 10 iterations')
print('via coo:', timeit('f(a,b,c,d)', number=10, globals={'f': viacoo, 'a': I, 'b': J, 'c': E, 'd': N}), 'secs for 10 iterations')

指纹：

correct: True
N: 100000  --  nnz (with duplicates): 4000000
direct:  7.702431229001377 secs for 10 iterations
via coo: 41.813509466010146 secs for 10 iterations

加速比:5倍