Unity中MathNET求解大型复稀疏矩阵的硬件加速性能问题

背景 #

我在Unity 6项目中使用MathNET库，通过NuGetForUnity工具成功安装了MathNET V5且运行正常。

目标 #

利用硬件加速（GPU或多核CPU）求解规模约为10k*10k的大型复稀疏矩阵方程组Ax=b。

当前问题 #

使用直接或迭代求解器求解稀疏矩阵时，无法利用硬件加速。

已尝试方案 #

• CUDA与openBLAS：未能在Unity中正确安装，导致MathNET无法检测到；
• MKL：安装成功并被MathNET识别，但仅对稠密矩阵（DenseMatrix）有效，对稀疏矩阵无效。

待解答问题 #

1. MathNET是否支持复稀疏矩阵的硬件加速？若支持，对应哪些原生提供商？
2. 如何在Unity中为MathNET安装CUDA和openBLAS提供商？

附加环境 #

• Unity版本：6
• 已安装NuGet包：
  • MathNet.Numerics.MKL.Win-x64.3.0.0（需勾选Load On Startup才能被MathNET识别）
  • MathNet.Numerics.Providers.CUDA.5.0.0
  • MathNet.Numerics.Providers.MKL.5.0.0
  • MathNet.Numerics.Providers.OpenBLAS.5.0.0

性能测试示例代码 #

using UnityEngine;
using MathNet.Numerics.LinearAlgebra.Complex;
using MathNet.Numerics;
using System.Diagnostics;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Statistics;
using MathNet.Numerics.LinearAlgebra.Solvers;
using MathNet.Numerics.LinearAlgebra.Complex.Solvers;
using System.Threading;

using Debug = UnityEngine.Debug;
using Complex = System.Numerics.Complex;

public class PerformanceTest : MonoBehaviour
{
    private Thread _t;
    void Start()
    {
        // 降低Unity渲染负载
        Application.targetFrameRate = 30;

        // 配置MathNet控制参数
        Control.ConfigureAuto();
        Debug.Log($"Native MKL : {Control.TryUseNativeMKL()}");
        Debug.Log(Control.Describe());

        // 在独立线程中执行计算，避免Unity长时间卡顿
        Debug.Log("启动性能测试任务。");
        _t = new Thread( () => {
            Thread.CurrentThread.IsBackground = false;
            PerfTest(); 
        });
        _t.Start();
    }

    private void PerfTest()
    {
        // 配置计时器用于精确测量
        Stopwatch sw = new Stopwatch();

        // 生成N*N矩阵M和稠密向量b，求解Mx=b
        sw.Start();
        int N = 100000;
        //DenseMatrix M = new DenseMatrix(N, N);
        SparseMatrix M = new SparseMatrix(N, N);
        DenseVector b = new DenseVector(N);
        for (int k = 0; k < N; k++)
        {
            // 生成一个示例用的带状矩阵
            M[k, k] = new Complex(k / N, -k / N);
            b[k] = new Complex(2 * k / N, -2 * k / N);
            if (k > 0) M[k, k - 1] = new Complex(-5f, 0f);
            if (k < N - 1) M[k, k + 1] = new Complex(-5f, 0f);
        }
        Debug.Log($"生成{ N }*{ N }矩阵耗时：{ TickToMS(sw.ElapsedTicks) } ms");
        sw.Stop();

        // 尝试用LU分解求解并验证结果
        sw.Restart();
        var x = M.LU().Solve(b);
        Debug.Log($"LU分解求解耗时：{ TickToMS(sw.ElapsedTicks) } ms");
        var err = M * x - b;
        err.MapInplace(c => c.Magnitude, Zeros.AllowSkip);
        double err_max = err.Enumerate().MaximumMagnitudePhase().Real;
        Debug.Log($"最大误差：{ err_max }");

        // 尝试用迭代求解器
        sw.Restart();
        var iterationCountStopCriterion = new IterationCountStopCriterion<Complex>(1000);
        var residualStopCriterion = new ResidualStopCriterion<Complex>(1e-10);
        var monitor = new Iterator<Complex>(iterationCountStopCriterion, residualStopCriterion);
        var solver = new BiCgStab(); // 可在此更换迭代求解器类型
        x = M.SolveIterative(b, solver, monitor);
        Debug.Log($"迭代求解器耗时：{ TickToMS(sw.ElapsedTicks) } ms");
        err = M * x - b;
        err.MapInplace(c => c.Magnitude, Zeros.AllowSkip);
        err_max = err.Enumerate().MaximumMagnitudePhase().Real;
        Debug.Log($"最大误差：{ err_max }");
    }

    private float TickToMS(long tick)
    {
        return (1000f * tick)/(float)Stopwatch.Frequency;
    }

    private void OnDisable()
    {
        if (_t != null && _t.IsAlive)
            _t.Abort();
    }
}

内容的提问来源于stack exchange，提问作者Eledwin