这个过程并不完全是微不足道的，但您可能有兴趣在此处观看视频讲座 21.6 和 21.65： https ://www.math.colostate.edu/~bangerth/videos.html

我为结构力学问题做了类似的事情，但我没有删除任何行或列。我通过修改我想要施加一些 BC (Dirichlet) 的自由度的对角项来使用惩罚方法。

我不确切知道您是如何组装全局稀疏刚度矩阵的，但我会跟踪全局刚度矩阵中属于第 i 个自由度的条目的位置。所以基本上是第 i 个自由度的对角线条目落在哪里的查找表，我可以准确地修改该值而不会干扰矩阵的其余部分。希望有帮助。

这并不难，并且有一个 fortran 代码示例可以说明这一点。它来自 John Burkardt 的代码示例集合。C 和 C++ 版本也可通过在同一集合中的简单搜索获得。

这是坐标格式，或者他称之为稀疏三元组。通过简单的修改来存储指向开始行的元素的指针ia而不是行号，您可以获得 CSR 版本。

该过程是将对角元素的位置存储在大小为 的稀疏矩阵元素数组中nnz。如果它对应于 Dirichlet 边界处的节点/自由度，则访问它并将其设置为 1，而右侧向量中的相应元素获取 Dirichlet BC 的值。

subroutine dirichlet_apply_dsp ( node_num, node_xy, node_condition, &
  nz_num, ia, ja, a, f )

!*****************************************************************************80
!
!! DIRICHLET_APPLY_DSP accounts for Dirichlet boundary conditions.
!
!  Discussion:
!
!    It is assumed that the matrix A and right hand side F have already been
!    set up as though there were no boundary conditions.  This routine
!    then modifies A and F, essentially replacing the finite element equation
!    at a boundary node NODE by a trivial equation of the form
!
!      A(NODE,NODE) * U(NODE) = NODE_BC(NODE)
!
!    where A(NODE,NODE) = 1.
!
!    This routine assumes that the coefficient matrix is stored in a
!    sparse triplet format.
!
!    This routine implicitly assumes that the sparse matrix has a storage
!    location for every diagonal element...or at least for those diagonal
!    elements corresponding to boundary nodes.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 July 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) NODE_NUM, the number of nodes.
!
!    Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of nodes.
!
!    Input, integer ( kind = 4 ) NODE_CONDITION(NODE_NUM), reports the
!    condition used to set the unknown associated with the node.
!    0, unknown.
!    1, finite element equation.
!    2, Dirichlet condition;
!    3, Neumann condition.
!
!    Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero entries.
!
!    Input, integer ( kind = 4 ) IA(NZ_NUM), JA(NZ_NUM), the row and column
!    indices of the nonzero entries.
!
!    Input/output, real ( kind = 8 ) A(NZ_NUM), the coefficient matrix,
!    stored in sparse triplet format; on output, the matrix has been adjusted
!    for Dirichlet boundary conditions.
!
!    Input/output, real ( kind = 8 ) F(NODE_NUM), the right hand side.
!    On output, the right hand side has been adjusted for Dirichlet
!    boundary conditions.
!
  implicit none

  integer ( kind = 4 ) node_num
  integer ( kind = 4 ) nz_num

  real ( kind = 8 ), dimension(nz_num) :: a
  integer ( kind = 4 ) column
  integer ( kind = 4 ), parameter :: DIRICHLET = 2
  real ( kind = 8 ), dimension(node_num) :: f
  integer ( kind = 4 ) ia(nz_num)
  integer ( kind = 4 ) ja(nz_num)
  integer ( kind = 4 ) node
  real ( kind = 8 ), dimension ( node_num ) :: node_bc
  integer ( kind = 4 ) node_condition(node_num)
  real ( kind = 8 ), dimension(2,node_num) :: node_xy
  integer ( kind = 4 ) nz
!
!  Retrieve the Dirichlet boundary condition value at every node. 
!
  call dirichlet_condition ( node_num, node_xy, node_bc ) ! SUPPLY THIS FUNCTION
!
!  Consider every matrix entry, NZ.
!
!  If the row I corresponds to a boundary node, then
!  zero out all off diagonal matrix entries, set the diagonal to 1,
!  and the right hand side to the Dirichlet boundary condition value.
!
  do nz = 1, nz_num

    node = ia(nz)

    if ( node_condition(node) == DIRICHLET ) then

      column = ja(nz)

      if ( column == node ) then
        a(nz) = 1.0D+00
        f(node) = 0.0D0 ! node_bc(node) FIXME
      else
        a(nz) = 0.0D+00
      end if

    end if

  end do

  return
end
```