ипичный случай использования для систем уравнений 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(vals[0, 0, 0] + vals[1, 1, 2])
I.append(nodes_edge_cells[0, 0])  # == nodes_edge_cells[1, 2]
J.append(nodes_edge_cells[0, 0])  # == nodes_edge_cells[1, 2]

Это делаетI, J, V короче и таким образом ускоряетtocsr.

Прямо сейчас, края "на ячейку". Я мог бы идентифицировать равные края друг с другом, используяnumpy.uniqueэффективно экономя около половиныI, J, V, Однако я обнаружил, что это тоже занимает некоторое время. (Неудивительно.)

Еще одна мысль, которая у меня была, заключалась в том, что я мог бы заменить диагональV, I, J простымnumpy.add.at если бы былcsr_matrix-подобная структура данных, где главная диагональ хранится отдельно. Я знаю, что это существует в некоторых других программных пакетах, но не смог найти его в scipy. Правильный?

Возможно, есть разумный способ построить КСО напрямую?