使用 Eigen 矩阵类型，其中行数和列数在编译时被编码到类型中，使您比 LAPACK 更具优势，其中矩阵大小仅在运行时才知道。这些额外信息允许编译器进行完整或部分循环展开，从而消除大量分支指令。如果您正在考虑使用现有库而不是编写自己的内核，那么拥有可以将矩阵大小作为 C++ 模板参数包含在内的数据类型可能是必不可少的。我所知道的唯一一个这样做的库是blaze，所以这可能值得对 Eigen 进行基准测试。

如果您决定推出自己的实现，您可能会发现 PETSc 对其块 CSR 格式所做的工作是一个有用的示例，尽管 PETSc 本身可能不是您所想的正确工具。他们没有写一个循环，而是明确地写出小矩阵向量乘法的每一个操作（参见他们的存储库中的这个文件）。这保证了没有像循环一样的分支指令。带有 AVX 指令的代码版本是如何实际使用向量扩展的一个很好的例子。例如，此函数使用__m256d数据类型同时对四个双精度进行同时操作。通过使用向量扩展明确写出所有操作，您可以获得可观的性能提升，仅用于 LU 分解而不是矩阵向量乘法。与其实际手动编写 C 代码，不如使用脚本生成它。当您重新排序某些操作以更好地利用指令流水线时，看看是否存在明显的性能差异也可能很有趣。

您还可以从工具STOKE中获得一些好处，它会随机探索可能的程序转换空间以找到更快的版本。

另一个想法可能是使用生成方法（编写程序的程序）。编写一个（元）程序，该程序吐出 C/C++ 指令序列以在 10x10 系统上执行 unpivoted** LU。基本上采用 k/i/j 循环嵌套并将其展平为 O(1000) 左右的行标量算术。然后将该生成的程序提供给任何优化编译器。我认为这里有趣的是，删除循环会暴露每个数据依赖项和冗余子表达式，并为编译器提供最大的机会来重新排序指令，以便它们很好地映射到实际硬件（例如执行单元的数量、危险/停顿，所以在）。

如果您碰巧知道所有矩阵（甚至只是其中的几个），则可以通过调用 SIMD 内在函数/函数 (SSE/AVX) 而不是标量代码来提高吞吐量。在这里，您将利用跨实例的令人尴尬的并行性，而不是在单个实例中追逐任何并行性。例如，您可以使用 AVX256 内在函数同时执行 4 个双精度 LU，方法是“跨”寄存器打包 4 个矩阵并对所有矩阵执行相同的操作**。

** 因此，重点放在未透视的 LU 上。旋转以两种方式破坏了这种方法。首先，它由于选择枢轴而引入了分支，这意味着您的数据依赖关系并不完全为人所知。其次，这意味着不同的 SIMD“插槽”必须做不同的事情，因为实例 A 可能与实例 B 不同。因此，如果您追求其中任何一个，我建议在计算之前静态旋转您的矩阵（置换最大条目每列到对角线）。

您的问题导致了两种不同的考虑。

首先，您需要选择正确的算法。因此，应该考虑矩阵是否具有任何结构的问题。例如，当矩阵对称时，Cholesky 分解比 LU 更有效。当您只需要有限的准确度时，迭代方法可以更快。

其次，您需要有效地实现算法。为此，您需要了解算法的瓶颈。您的实现是否受内存传输速度或计算速度的限制。既然你只考虑$10 \times 10$矩阵，您的矩阵应该完全适合 CPU 缓存。因此，您应该利用处理器的 SIMD 单元（SSE、AVX 等）和内核，在每个周期内进行尽可能多的计算。

总之，您的问题的答案很大程度上取决于您考虑的硬件和矩阵。可能没有明确的答案，您必须尝试一些事情才能找到最佳方法。

我会尝试分块反转。

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen 使用优化的例程来计算 4x4 矩阵的逆矩阵，这可能是您将得到的最好的。尝试尽可能多地使用它。

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

左上角：8x8。右上角：8x2。左下角：2x8。右下角：2x2。使用优化的 4x4 反转代码反转 8x8。剩下的是矩阵产品。

编辑：使用 6x6、6x4、4x6 和 4x4 块已证明比我上面描述的要快一点。

using namespace Eigen;

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
    Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;

    Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
    Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<br_size, br_size>() = DCAB_inv;

    return result;
}

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
    const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
    const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
    const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
    const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();

    return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
    const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
    const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
    const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();

    Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
    Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<2, 2>() = DCAB_inv;

    return result;
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
    const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
    const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
    const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();

    Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
    Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();

    result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<4, 4>() = DCAB_inv;

    return result;
}

Eigen::Matrix<double,10,10>::Random()这是使用一百万个矩阵和Eigen::Matrix<double,10,1>::Random()向量进行的一次基准测试的结果。在我所有的测试中，我的逆总是更快。我的求解例程涉及计算倒数，然后将其乘以向量。有时它比 Eigen 快，有时它不是。我的基准标记方法可能有缺陷（没有禁用涡轮增压等）。此外，Eigen 的随机函数可能不代表真实数据。

• 特征部分枢轴反转：3036 毫秒
• 我的逆 8x8 上块：1638 毫秒
• 我的逆 6x6 上块：1234 毫秒
• 特征部分枢轴求解：1791 毫秒
• 我用 8x8 上块解决：1739 毫秒
• 我用 6x6 上块解决：1286 毫秒

我很想看看是否有人可以进一步优化这一点，因为我有一个有限元应用程序，它可以反转无数个 10x10 矩阵（是的，我确实需要逆的各个系数，所以直接求解线性系统并不总是一种选择） .