您跟踪您发现的 i,j 交互的想法是可行的，我认为这就是您和 Stefano M 所指的“CS 技巧”。这相当于以列表格式构造您的稀疏矩阵。

不确定您有多少 CS，所以如果您已经知道这一点，我深表歉意：在链表数据结构中，每个条目都存储一个指向它之后的条目和之前的条目的指针。从中添加和删除条目很便宜，但在其中查找项目并不那么简单——您可能必须查看所有条目。

因此，对于每个节点 i，您存储一个链表。然后遍历所有元素；如果你发现两个节点 i 和 j 是连接的，你去查看 i 的链表。如果 j 不存在，则将其添加到列表中，同样将 i 添加到 j 的列表中。如果按顺序添加它们是最简单的。

填充列表列表后，您现在知道矩阵每行中非零条目的数量：它是该节点列表的长度。这些信息正是您在 PETSc 的矩阵数据结构中预分配稀疏矩阵所需的信息。然后您可以释放您的列表列表，因为您不再需要它。

但是，这种方法假定您所拥有的只是每个元素包含的节点的列表。

一些网格生成包（例如Triangle）不仅可以输出元素列表及其包含的节点，还可以输出三角剖分中每条边的列表。在这种情况下，您不会冒高估非零条目数量的风险：对于分段线性元素，每条边都会给您正好 2 个刚度矩阵条目。您使用的是分段二次，因此每条边都计为 4 个条目，但您明白了。在这种情况下，您可以使用普通数组通过边缘列表找到每行的非零条目数。

使用这种方法，您必须从硬盘读取一个额外的大文件，如果您的实际计算不是那么大，这实际上可能比使用元素列表要慢。尽管如此，我认为它更简单。

如果您将网格指定为 DMPlex，将数据布局指定为 PetscSection，则 DMCreateMatrix() 将自动为您提供正确的预分配矩阵。以下是泊松问题和斯托克斯问题的 PETSc 示例。

赞成

我个人不知道有什么便宜的方法可以做到这一点，所以我只是高估了这个数字，即对所有行使用一个相当大的值。

例如，对于由线性 8 节点十六进制元素组成的完美结构网格，对角线块和非对角线块中每行的 nnzs 都是自由度*27。对于大多数完全非结构化的自动生成的六角网格，该数量很少超过 dof*54。对于线性 tets，我从来不需要超过 dof*30。对于某些形状非常糟糕/纵横比元素非常低的网格，您可能必须使用稍大的值。

代价是本地（按等级）内存消耗介于 2x-5x 之间，因此您可能不得不在集群上使用比平时更多的计算节点。

顺便说一句，我确实尝试使用可搜索列表，但确定稀疏结构所花费的时间超过了组装/求解。但是我的实现非常简单，并且没有使用有关边缘的信息。

另一种选择是使用 DMMeshCreateExodus 之类的例程，如本例所示。

您正在寻找枚举所有唯一 (gi,gj) 连接，这建议将它们全部放入（非重复）关联容器中，然后计算其基数 - 在 C++ 中，这将是 std::set < std::pair < 整数，整数 > >。在您的伪代码中，您将“nnz[i]++”替换为“s.insert[pair(gi,gj)]”，然后非零的最终数量是 s.size()。它应该在 O(n-log-n) 时间内运行，其中 n 是非零的数量。

由于您可能已经知道可能的 gi 的范围，因此您可以通过 gi 索引“展开”表格以提高性能。这将用 std::vector < std::set < int >> 替换您的集合。你用“v[gi].insert(gj)”填充它，然后非零的总数来自对所有 gi 求和的 v[gi].size()。这应该在 O(n-log-k) 时间内运行，其中 k 是每个元素的未知数（对您来说是六个 - 对于大多数 pde 代码来说基本上是一个常数，除非您在谈论 hp 方法）。

（注意 - 希望这是对所选答案的评论，但太长了 - 抱歉！）