我将 SciPysparse.csr_matrix格式用于有限元代码。在应用基本边界条件时，我在右侧向量中设置所需的值，并将相应的行和列设置为标识（即对角线上的 1 和其他地方的零。）
我使用的代码类似于：
K = csr_matrix((data,(rows,cols)))
for idx in boundary:
K[:,idx] = 0.0
K[idx,:] = 0.0
K[idx,idx] = 1.0
此代码导致为 idx 行和列中的所有条目插入值，这非常慢。
我希望能够将对角线条目设置为 1（这不应该改变稀疏性），并且仅将idx行和列中尚未为零的条目设置为零。
有没有办法做到这一点？
我认为直接处理 CSR 矩阵的数据是最有效的。看看这个CSR 数据结构的解释——它是德文的，但是这个例子可以帮助你理解。
下面是一个 python 脚本，它的函数设置 CSR 矩阵的指定行以匹配身份：
In [13]: run t_csrdiag.py
sparse mat A was before
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]
[12 13 14 15]]
index of the row that is modified: 2
sparse mat A is now
[[ 0 1 2 3]
[ 4 5 6 7]
[ 0 0 1 0]
[12 13 14 15]]
随意循环或矢量化它。
这是代码：
import scipy.sparse as sps
import numpy as np
def main():
# index of the row to become the identity row
rowind = 2
A = np.arange(16).reshape((4, 4))
A = sps.csr_matrix(A)
print 'sparse mat A was before'
print A.todense()
print '\nindex of the row that is modified: {0}'.format(rowind)
setcsrrow2id(A, rowind)
print '\nsparse mat A is now'
print A.todense()
def setcsrrow2id(amat, rowind):
indptr = amat.indptr
values = amat.data
indxs = amat.indices
# get the range of the data that is changed
rowpa = indptr[rowind]
rowpb = indptr[rowind+1]
# new value and its new rowindex
values[rowpa] = 1.0
indxs[rowpa] = rowind
# number of new zero values
diffvals = rowpb - rowpa - 1
# filter the data and indices and adjust the range
values = np.r_[values[:rowpa+1], values[rowpb:]]
indxs = np.r_[indxs[:rowpa+1], indxs[rowpb:]]
indptr = np.r_[indptr[:rowind+1], indptr[rowind+1:]-diffvals]
# hard set the new sparse data
amat.indptr = indptr
amat.data = values
amat.indices = indxs
if __name__ == '__main__':
main()
您可以使用该方法nonzero返回给定行或列的所有非零条目的数组，例如：
for i in boundary:
_, J = A[i,:].nonzero() # ignored the first return value, which is just [i,...,i]
for j in J:
if j!=i:
A[i,j] = 0.0
A[i,i] = 1.0
当然有一个摆脱条件的聪明方法，但你明白了。
另外，我认为您只想将i与边界中的索引对应的行设置为标识，而不是列。如果节点重新排序以使边界节点位于开头，则矩阵应如下所示
其中表示内部节点之间的耦合和边界条件对内部节点的影响。
这是我非常简单的功能，可以消除bc_id矩阵行中的所有条目，只保留对角线条目。这就是我通常设置边界条件的方式。我密集地使用它。我知道它不是最优的，它不是很优雅，但是每当我分析我的代码时，这个函数从来都不是问题。可能最好包含一个断言来检查矩阵是否是平方的。
from scipy import sparse
import numpy as np
import matplotlib.pyplot as plt
def set_diag(A,bc_id):
ndofs = A.shape[0]
diago = A.diagonal()
uno = np.ones((ndofs))
uno[bc_id] = 0
uno = sparse.dia_matrix((uno,0), shape = (ndofs,ndofs))
# up to here I delete the rows
# I multiply A by an identity matrix
# where i set to zero the rows I want
# to delete
A = uno.dot(A)
new_diag_entries = np.zeros((ndofs))
# here I set the diagonal entries
# equals to the value on the diagonal.
# Use the second line if you want to
# to set the diagonal entry to one
new_diag_entries[bc_id] = diago[bc_id]
#new_diag_entries[bc_id] = uno[bc_id]
uno = sparse.dia_matrix((new_diag_entries,0), shape = (ndofs,ndofs))
A = A+uno # here I set the diagonal entry
return A
A = sparse.csr_matrix([[1, 2, 3, 10], [4,5, 6,11], [7, 8, 9,12],[13,14,15,16]])
bc_id = np.array([2,3])
A = set_diag(A,bc_id)
print A
您可以根据 csr 矩阵乘法和矩阵加法在纯线性代数级别对此进行公式化，这很快。我发现这是我笔记本电脑上 scipy 中最快和最节省内存的方法。
具体来说，让是原始矩阵，并且是修改后的版本，其中与边界节点关联的行和列清零，除了放置值 1 的对角线条目。然后
在哪里是对角矩阵，其位置对应于内部自由度，其他位置为零，反之亦然是对角矩阵，对应于边界自由度的位置为 1，其他位置为 0。
因此，总体策略是：
形成对角矩阵和.
兑换和为 CSR 格式。
形成修改后的矩阵.
这是我使用的代码：
def dirichlet_zero_matrix_modification(matrix, dirichlet_zero_dofs):
"""Enforces dirichlet zero B.C.'s by zeroing-out rows and columns
associated with dirichlet dofs, and putting 1 on the diagonal there."""
N = matrix.shape[1]
chi_interior = np.ones(N)
chi_interior[dirichlet_zero_dofs] = 0.0
I_interior = sps.spdiags(chi_interior, [0], N, N).tocsr()
chi_boundary = np.zeros(N)
chi_boundary[dirichlet_zero_dofs] = 1.0
I_boundary = sps.spdiags(chi_boundary, [0], N, N).tocsr()
matrix_modified = I_interior * matrix * I_interior + I_boundary
return matrix_modified
根据我使用 scipy 的经验，与实际操作矩阵条目或切片的其他方法相比，用更高级别的线性代数运算（在可能的情况下）表达矩阵操作几乎总是产生更好的性能。在其他语言和框架中，情况可能并非如此。