现有串行代码矩阵分解计算多线程使用OpenMP案例

Question

我遇到了一个代码，该代码使用低性能系列 for 循环来计算所谓的“Crout 分解”。如果您需要了解更多，我认为这个概念类似于 described here。乍一看，代码是一个 for 循环，可以通过 omp 指令简单地变为并行：

SparseMtrx *Skyline :: factorized()
{
    // Returns the receiver in  U(transp).D.U  Crout factorization form.
        ...
        // this loop is very expensive and only uses single thread:
        for ( int k = 2; k <= n; k++ ) {
            int ack = adr.at(k);
            int ack1 = adr.at(k + 1);
            int acrk = k - ( ack1 - ack ) + 1;
            for ( int i = acrk + 1; i < k; i++ ) {
                int aci = adr.at(i);
                int aci1 = adr.at(i + 1);
                int acri = i - ( aci1 - aci ) + 1;
                int ac;
                if ( acri < acrk ) {
                    ac = acrk;
                } else {
                    ac = acri;
                }
    
                int acj = k - ac + ack;
                int acj1 = k - i + ack;
                int acs = i - ac + aci;
                double s = 0.0;
                for ( int j = acj; j > acj1; j-- ) {
                    s += mtrx [ j ] * mtrx [ acs ];
                    acs--;
                }
    
                mtrx [ acj1 ] -= s; //mtrx  here is the matrix values (shared)
            }
            double s = 0.0;
            for ( int i = ack1 - 1; i > ack; i-- ) {
                double g = mtrx [ i ];
                int acs = adr.at(acrk);
                acrk++;
                mtrx [ i ] /= mtrx [ acs ];
                s += mtrx [ i ] * g;
            }
    
            mtrx [ ack ] -= s;
        }
    ...

但问题是似乎每一步都使用矩阵值的更新值 mtrx 来计算下一个分解值。

我的问题是有什么方法可以使用多线程来做这个计算吗？

完整代码可通过 this link 在 github 中获得。

Answer 1

Crout 分解是高斯消去法的一种变体。您可以通过 1. 它们的 end-product 2. 它们的处理方式来描述此类算法。 Crout 分解计算 LU，其中 U 的对角线是恒等式。其他因式分解对 L 的对角线进行归一化，或者它们计算 LDU 时 L,U 都被归一化，或者它们计算 LU 时 L 和 U 的对角线相等。所有这些都与并行化无关。

接下来，您可以通过计算方式来描述高斯消除算法的特征。这些都是基本三重嵌套循环在数学上等价的re-organizations。由于它们在数学上是等价的，您可以选择一个或另一个以获得有利的计算特性。

要解决一件事：高斯消元法在计算枢轴时具有本质上的顺序成分，因此您无法使其完全并行。（好吧，有些关于行列式的嘀咕，但那是不现实的。）外循环是连续的，步数等于矩阵大小。问题是你是否可以使内部循环并行。

Crout 算法的特点是，在分解的步骤“k”中，它计算 U 的行“k”和 L 的列“k”。由于两者中的元素不递归相关，因此这些更新可以在并行循环中完成，它为您提供了一个顺序的“k”循环，并且对于每个 k 值，有两个并行的单循环。这就剩下第三个循环级别：那个循环级别来自这样一个事实，即每个独立的迭代都涉及一个总和，因此是一个缩减。

这听起来不太适合并行化：如果内部是缩减，则不能 collapse 两个循环，因此您必须选择要并行化的级别。也许您应该使用不同的高斯消去公式。例如，普通算法基于“rank-k 更新”，并且这些算法非常并行。该算法有一个顺序循环，每个步骤都有一个 collapse(2) 并行循环。

现有串行代码矩阵分解计算多线程使用OpenMP案例

Case of using OpenMP for multi-threading of a matrix factorization calculation of an existing serial code

c++

openmp

sparse-matrix