2

我正在使用 repa 库（3 和 4）在 haskell 中处理点云。至少我正在努力。

在并行性确实有很大帮助的情况下，我需要大量执行一些操作。其中大多数是对点的（度量）邻域的简单线性代数运算。例如，我需要在每行都是一个点的小矩阵上计算 SVD 的主成分分析。

现在我使用线性包作为向量类型

type Vec3 = V3 Float

以及用于点云的这些向量的一维数组

type Cloud = Array F DIM1 Vec3

所以现在我遇到了使用computeP并行调用hmatrix矩阵分解的函数的问题。我为此尝试了 hmatrix 本身以及 repa -linear-algebra包。我遇到的问题是，对于所有这些调用（无论我如何提供数据，无论我调用什么（svd、特征分解、qr decomp 等）），应用程序总是会因总线错误或段错误而随机崩溃。

我也没有找到任何方法来获得至少可以为我指明正确方向的任何堆栈跟踪。堆栈跟踪通常在 pthread 结束。

此外，我编写了自己的 C 代码，我称之为：

foreign import ccall safe "pca.hpp pca"
    c_pca :: CUInt -> Ptr Float -> Ptr Float -> Ptr Float -> IO ()

{-# INLINE foreignPCA #-}
foreignPCA :: forall r . (Source r Vec3) => Array r DIM1 Vec3 -> ([Vec3], Vec3)
foreignPCA !vs = unsafePerformIO $ do
    n <- return $ Arrays.length vs
    ps <- mallocForeignPtrArray n :: IO (ForeignPtr Vec3) -- point matrix
    computeIntoS ps (delay vs)
    as <- mallocForeignPtrArray 3 :: IO (ForeignPtr Float) -- singular values
    av <- mallocForeignPtrArray 3 :: IO (ForeignPtr Vec3) -- right singular vectors
    withForeignPtr (castForeignPtr ps) $ \pps ->
        withForeignPtr as $ \pas ->
            withForeignPtr (castForeignPtr av) $ \pav -> do
                c_pca (fromIntegral n) pps pas pav
                svalues <- peekArray 3 pas :: IO [Float]
                svecs <- peekArray 3 (castPtr pav :: Ptr Vec3) :: IO [Vec3]
                let [sx, sy, sz] = svalues in
                    return (svecs, (V3 sx sy sz))

这在具有 20 个并行内核的大型点云上运行良好。从未以任何方式坠毁。

现在我非常模糊的想法是 hmatrix 用“安全”调用 C/Fortran 代码，从而允许 pthread 分叉并且实际上不是线程安全的。我无法尝试验证这个假设，因为调试对于 haskell 工具链来说似乎是一个陌生的概念（至少对于我这样一个完整的新手来说）。

总之，我有三个问题：

  • 是否已知 hmatrix 存在并行工作的问题
  • 有没有人致力于这些基本算法的本地实现？
  • 如何防止 FFI 包装的代码在无法访问导入调用的情况下生成 fork() 的线程？
  • 如何调试 hmatrix 之类的东西？

第二个是我特别感兴趣的，因为我发现 hmatrix 非常丑陋（subhask 看起来很有希望，但太不完整以至于不可行）。我的目标是改用haskell，但如果我必须使用我自己的C++代码来处理上述任何琐碎的事情，我可以像现在一样继续用C++编码......

4

0 回答 0