您应该能够使用 LAPACK/BLAS、使用 QR 分解 [geqrf]、正交乘法 [ormqr] 和三角求解 [trsm] 例程有效地做到这一点。LAPACK 和 BLAS 非常重视存储经济性，因此这些算法中的许多算法就地运行，用它们的输出覆盖它们的输入（geqrf、ormqr 和 trsm 都是这样工作的）。

你应该从 geqrf 开始，它会用\和。的输出只是直接列在上三角形中，但以“反射器形式”紧凑地表示。这种表示需要少量/矢量大小的额外工作区（），您提供给 geqrf 的空存储，它将在退出时填充。您需要挂在做任何其他事情。$\mathbf B$$\mathbf Q$$\mathbf R$$\mathbf R$$\mathbf Q$$\mathbf \tau$$\mathbf \tau$$\mathbf Q$

因为以特殊/反射器形式存储，所以不能只使用 BLAS/gemm 来应用它。相反，您应该使用 LAPACK 例程 ormqr，它可以从左侧或右侧将应用到并应用转置。它也操作就地，用（例如）或。这可能会使推理有点棘手，因为（例如）和的大小不一定相同（$\mathbf Q$$\mathbf Q$$\mathbf A$$\mathbf A$$\mathbf Q^T \mathbf A$$\mathbf A \mathbf Q$$\mathbf Q^T \mathbf A$$\mathbf A$$\mathbf Q$高/瘦）。尽管如此，您应该能够使用对 ormqr 的两次调用（从 'L'eft 使用 'T'transpose，然后从 'R'ight 使用 'N'otranspose）来获得的左上角。$\mathbf C = \mathbf Q^T\mathbf A\mathbf Q$$\mathbf A$

最后一步，形成可以使用 BLAS 例程 trsm 完成。它应用三角系统的逆，用（例如）输出 mathbf C 。它可以从左或右你需要两个调用，一个来自'L'eft 的'N'otranspose，另一个来自'R'right 的'T'ranspose。在这两种情况下，您都需要使用 uplo='U'pper 因为那是 geqrf 存储的地方。$\mathbf R ^{-1} \mathbf C \mathbf R^{-T}$$\mathbf C$$\mathbf R^{-1} \mathbf C$$\mathbf R$$\mathbf R$

不幸的是，我不是 C# 专家，但我想有一些常用的方法来调用本机/非托管“extern C”库，这就是你想要开始的地方。