问题本身的数值成本线性依赖于行数。如果您单独解决这些块，您将达到这些最低成本。不利用系统内部结构的常用矩阵求解器会导致数值成本随行数的三次方增长。所以，它们更贵。

但是，有一些通用求解器可以检测矩阵结构，然后以线性成本求解系统。矩阵结构的检测具有行数二次方的成本，但影响非常低，因为测试一个矩阵元素所需的时间非常低。我们正在为SimulationX使用这样的求解器。使用此求解器，求解线性系统的时间很少成为模拟中的限制因素。

如果只有一百万个，您应该能够在一台像样的笔记本电脑或台式机上在大约 1 秒内轻松解决所有问题，因此性能甚至不应该成为问题（除非 1 秒对于您的预期应用程序来说仍然太慢）。

为了说明这一点，请考虑以下 Mathematica 代码，它解决了 100 万个随机系统：$3\times 3$

First@Timing@Table[LinearSolve[RandomReal[{0, 1}, {3, 3}]], {1000000}]

在我使用单个 CPU 内核的 7 年前的廉价笔记本电脑上，这需要 30 秒。在现代台式计算机上使用 4-6 个并行 CPU 内核，这大约需要 1 秒。在具有 GPU 并行性的 CUDA 中，您可能可以在几毫秒内完成。

如果您对矩阵有所了解，那么使用它几乎总是值得的。在您的情况下，使用您知道它由许多块组成的事实并单独解决它们 - 例如，您可以并行执行。如果您只是将一个大块对角矩阵提供给迭代或稀疏直接求解器，那么他们能做的最好的事情就是分析矩阵并完全按照您自己的方式做，但是这个分析步骤（i）通常不做，至少在迭代求解器中，（ii）昂贵。$3\times 3$

我认为最好解决单独的 3 x 3 系统。做这么大的矩阵，包含数百万个 3 x 3 系统的问题：

1. 组装的计算成本：组装如此大的矩阵可能需要大量时间（即使使用稀疏求解器）。
2. 并行化开销：如果需要并行化代码，则必须将大矩阵分布在许多处理器上，并且在处理器之间传递信息会产生通信开销。另一方面，对于单独的 3 x 3 系统，系数可以很容易地分布在各个处理器上，之后将是“令人尴尬的并行”操作（即处理器之间不需要通信）。此外，使用 OpenMP 或 MPI 进行并行化也不会太困难。另一方面，对于大型矩阵，使用库（例如 PETSc：例如ex2.c示例是一个不错的开始）可能比编写自己的代码更容易。

最后，我想指出，对于 3 x 3 系统，使用计算机代数系统（或使用 Wolfram Alpha）很容易得出分析公式。作为说明，这里是如何使用 Python 包 sympy 完成的：

from sympy import *
x1,x2,x3=symbols('x1 x2 x3')
b1,b2,b3=symbols('b1 b2 b3')
a11,a12,a13=symbols('a11 a12 a13')
a21,a22,a23=symbols('a21 a22 a23')
a31,a32,a33=symbols('a31 a32 a33')
eq1=Eq(a11*x1+a12*x2+a13*x3,b1)
eq2=Eq(a21*x1+a22*x2+a23*x3,b2)
eq3=Eq(a31*x1+a32*x2+a33*x3,b3)
x=solve([eq1,eq2,eq3], [x1,x2,x3])
x.get(x1)
x.get(x2)
x.get(x3)

只需将 x1、x2、x3、b1、b2、b3 等替换为现有代码中的实际变量即可。如果需要生成 C 代码：

[(c_name, c_code), (h_name, c_header)] = codegen(("x1",x.get(x1)), "C", "test", header=False, empty=False)
print(c_code)
[(c_name, c_code), (h_name, c_header)] = codegen(("x2",x.get(x2)), "C", "test", header=False, empty=False)
print(c_code)
[(c_name, c_code), (h_name, c_header)] = codegen(("x3",x.get(x3)), "C", "test", header=False, empty=False)
print(c_code)

示例输出（仅 x1）：

>>> print(c_code)
#include "test.h"
#include <math.h>
double x1(double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33, double b1, double b2, double b3) {
   return (a12*a23*b3 - a12*a33*b2 - a13*a22*b3 + a13*a32*b2 + a22*a33*b1 - a23*a32*b1)/(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
}

在这里，预先计算行列式以进一步降低计算成本可能是有意义的，即：

determinant = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);