使用 GMP/ARB 矩阵减少 OpenMP

Question

我想并行化我编写的一个程序，该程序计算涉及矩阵和向量乘积的序列，结果是一个向量。由于参数变得非常小和大，我使用 ARB（基于 GMP、MPFR 和 flint）来防止失去意义。此外，由于系列元素是独立的，矩阵维度并不大，但系列需要评估多达 50k 个元素，因此让多个线程每个计算系列的几个元素是有意义的，即 5 个线程每个可以计算 10k并行元素，然后将结果向量相加。

现在的问题是，用于将向量和矩阵相加的 ARB 函数不是可以在 openmp 缩减中轻松使用的标准操作。当我天真地尝试编写自定义缩减时，g++ 抱怨 void 类型，因为 ARB 中的操作没有 return 值：

void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec)¶

将计算并设置 z 为 z=x+y，精度为 prec 位，但 arb_add 本身是一个空函数。

举个例子：一个类似问题的随机for循环看起来像这样（当然我的实际程序是不同的）

[...]

arb_mat_t RMa,RMb,RMI,RMP,RMV,RRes;

arb_mat_init(RMa,Nmax,Nmax); //3 matrices
arb_mat_init(RMb,Nmax,Nmax);
arb_mat_init(RMI,Nmax,Nmax);
arb_mat_init(RMV,Nmax,1);  // 3 vectors
arb_mat_init(RMP,Nmax,1);
arb_mat_init(RRes,Nmax,1);

[...]

//Result= V + ABV +(AB)^2 V + (AB)^3 V + ...
//in my actual program A and B would be j- and k-dependent and would
//include more matrices and vectors

#pragma omp parallel for collapse(1) private(j)
for(j=0; j<jmax; j++){
        arb_mat_one(RMI); //sets the matrix RMI to 1
        for(k=0; k<j; k++){ 
            Qmmd(RMI,RMI,RMa,Nmax,prec); //RMI=RMI*RMa
            Qmmd(RMI,RMI,RMb,Nmax,prec); //RMI=RMI*RMb
            cout << "j=" << j << ",   k=" << k << "\r" << flush;
        }

        Qmvd(RMP,RMI,RMV,Nmax,prec);     //RMP=RMI*RMV
        arb_mat_add(RRes,RRes,RMP,prec); //Result=Result + RMP        
}

[...]

当使用超过 1 个线程时，这当然会崩溃，因为我没有指定减少 RRes。这里Qmmd()和Qmvd()是自写的矩阵-矩阵和矩阵-向量乘积函数，RMa、RMb、RMV是随机矩阵和向量，resp.

现在的想法是减少 RRes，这样每个线程都可以计算 RRes 的私有版本，包括最终结果的一小部分，然后再使用 arb_mat_add 将它们全部相加。我可以写一个函数 matrixadd(A,B) 来计算 A=A+B

void matrixadd(arb_mat_t A, arb_mat_t B) {
    arb_mat_add(A,A,B,2000); 
//A=A+B, the last value is the precision in bits used for that operation
}

然后最终

#pragma omp declare reduction \
(myadd : void : matrixadd(omp_out, omp_in))
#pragma omp parallel for private(j) reduction(myadd:RRes) 
for(j=0; j<jmax; j++){
        arb_mat_one(RMI);
        for(k=0; k<j; k++){
            Qmmd(RMI,RMI,RMa,Nmax,prec);
            Qmmd(RMI,RMI,RMb,Nmax,prec);
            cout << "j=" << j << ",   k=" << k << "\r" << flush;
        }

        Qmvd(RMP,RMI,RMV,Nmax,prec); 
        matrixadd(RRes,RMP);        
}

Gcc 对此不满意：

main.cpp: In function ‘int main()’:
main.cpp:503:46: error: invalid use of void expression
     (myadd : void : matrixadd(omp_out, omp_in))
                                          ^
main.cpp:504:114: error: ‘RRes’ has invalid type for ‘reduction’

Openmp 能否理解我的无效减少并且可以与 ARB 和 GMP 一起工作？如果是这样，如何？谢谢！

（此外，我的程序目前在 j-for 循环中包含一个带中断条件的收敛检查。如果您也知道如何轻松实现这样的事情，我将非常感激，因为我目前openmp 测试我刚刚删除了中断并设置了一个常量 jmax。）

我的问题与this one非常相似。

编辑：抱歉，这是我尝试的一个最小的、完整的和可验证的例子。其他必需的软件包是 arb, flint, gmp, mpfr (available through packetmanagers) and gmpfrxx.

#include <iostream>
#include <omp.h>
#include <cmath>
#include <ctime>

#include <cmath>
#include <gmp.h>
#include "gmpfrxx/gmpfrxx.h"

#include "arb.h"
#include "acb.h"
#include "arb_mat.h"

using namespace std;

void generate_matrixARBdeterministic(arb_mat_t Mat, int N, double w2) //generates some matrix
{
    int i,j;
    double what;
    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
        {
            what=(i*j+30/w2)/((1+0.1*w2)*j+20-w2);
            arb_set_d(arb_mat_entry(Mat,i,j),what);
        }
    }
}

void generate_vecARBdeterministic(arb_mat_t Mat, int N) //generates some vector
{
    int i;
    double what;
    for(i=0;i<N;i++)
    {
        what=(4*i*i+40)/200;
        arb_set_d(arb_mat_entry(Mat,i,0),what);
    }
}

void Qmmd(arb_mat_t res, arb_mat_t MA, arb_mat_t MB, int NM, slong prec)
{   ///res=M*M=Matrix * Matrix

    arb_t Qh1;
    arb_mat_t QMh;

    arb_init(Qh1);

    arb_mat_init(QMh,NM,NM);

    for (int i=0; i<NM; i++){
            for(int j=0; j<NM; j++){
                     for (int k=0; k<NM; k++ ) {
                        arb_mul(Qh1,arb_mat_entry(MA, i, k),arb_mat_entry(MB, k, j),prec);
                        arb_add(arb_mat_entry(QMh, i, j),arb_mat_entry(QMh, i, j),Qh1,prec);
                    }
             }
    }

    arb_mat_set(res,QMh);

    arb_mat_clear(QMh);
    arb_clear(Qh1);
  }

void Qmvd(arb_mat_t res, arb_mat_t M, arb_mat_t V, int NM, slong prec)  //res=M*V=Matrix * Vector
{ ///res=M*V
    arb_t Qh,Qs;
    arb_mat_t QMh;

    arb_init(Qh);
    arb_init(Qs);
    arb_mat_init(QMh,NM,1);

    arb_set_ui(Qh,0.0);
    arb_set_ui(Qs,0.0);
    arb_mat_zero(QMh);

    for (int i=0; i<NM; i++){
            arb_set_ui(Qs,0.0);
            for(int j=0; j<NM; j++){
                    arb_mul(Qh,arb_mat_entry(M, i, j),arb_mat_entry(V, j, 0),prec);
                    arb_add(Qs,Qs,Qh,prec);
             }
            arb_set(arb_mat_entry(QMh, i, 0),Qs);
    }
    arb_mat_set(res,QMh);

    arb_mat_clear(QMh);
    arb_clear(Qh);
    arb_clear(Qs);
}

void QPrintV(arb_mat_t A, int N){ //Prints Vector
    for(int i=0;i<N;i++){
                cout <<  arb_get_str(arb_mat_entry(A, i, 0),5,0) << endl; //ARB_STR_NO_RADIUS
    }
}

void matrixadd(arb_mat_t A, arb_mat_t B) {
    arb_mat_add(A,A,B,2000);
}

int main() {
    int Nmax=10,jmax=300; //matrix dimension and max of j-loop
    ulong prec=2000; //precision for arb

    //initializations

    arb_mat_t RMa,RMb,RMI,RMP,RMV,RRes;

    arb_mat_init(RMa,Nmax,Nmax);
    arb_mat_init(RMb,Nmax,Nmax);
    arb_mat_init(RMI,Nmax,Nmax);
    arb_mat_init(RMV,Nmax,1);
    arb_mat_init(RMP,Nmax,1);
    arb_mat_init(RRes,Nmax,1);

    omp_set_num_threads(1);
    cout << "Maximal threads is " << omp_get_max_threads() << endl;

    generate_matrixARBdeterministic(RMa,Nmax,1.0); //generates some Matrix for RMa
    arb_mat_set(RMb,RMa); // sets RMb=RMa

    generate_vecARBdeterministic(RMV,Nmax); //generates some vector

    double st=omp_get_wtime();

    Qmmd(RMI,RMa,RMb,Nmax,prec);
    int j,k=0;

    #pragma omp declare reduction \
    (myadd : void : matrixadd(omp_out, omp_in))
    #pragma omp parallel for private(j) reduction(myadd:RRes)
    for(j=0; j<jmax; j++){
            arb_mat_one(RMI);
            for(k=0; k<j; k++){
                Qmmd(RMI,RMI,RMa,Nmax,prec);
                Qmmd(RMI,RMI,RMb,Nmax,prec);
                cout << "j=" << j << ",   k=" << k << "\r" << flush;
            }

            Qmvd(RMP,RMI,RMV,Nmax,prec);  
            matrixadd(RRes,RMP);      
    }

    QPrintV(RRes,Nmax);

    double en=omp_get_wtime();
    printf("\n Time it took was %lfs\n",en-st);


    arb_mat_clear(RMa);
    arb_mat_clear(RMb);
    arb_mat_clear(RMV);
    arb_mat_clear(RMP);
    arb_mat_clear(RMI);
    arb_mat_clear(RRes);

    return 0;
}

和

g++ test.cpp -g -fexceptions -O3 -ltbb -fopenmp -lmpfr -lflint -lgmp -lgmpxx -larb -I../../PersonalLib -std=c++14 -lm -o b.out

Answer 1

您可以为每个线程创建一个矩阵数组来求和。您只需将 matrixadd(RRes, RMP) 替换为 matrixadd(RRes[get_omp_thread_num()], RMP)，然后在最后对所有 RRes 求和，其中 RRes 现在将是 std::vector<arb_mat_t>。

或者你可以尝试为包装器定义一个加法运算符class，当然你应该小心避免复制整个矩阵。这感觉更麻烦，因为您必须小心内存管理（因为您使用的是库 - 除非您花时间完成所有操作，否则您不知道它到底做了什么）。

Answer 2

你可以像这样手工缩小

#pragma omp parallel
{
    arb_mat_t RMI, RMP;
    arb_mat_init(RMI,Nmax,Nmax);  //allocate memory
    arb_mat_init(RMP,Nmax,1);     //allocate memory
    #pragma omp for
    for(int j=0; j<jmax; j++){
        arb_mat_one(RMI);
        for(int k=0; k<j; k++){
            Qmmd(RMI,RMI,RMa,Nmax,prec);
            Qmmd(RMI,RMI,RMb,Nmax,prec);
        }
    }
    Qmvd(RMP,RMI,RMV,Nmax,prec);
    #pragma omp critical
    arb_mat_add(RRes,RRes,RMP,prec);
    arb_mat_clear(RMI);  //deallocate memory
    arb_mat_clear(RMP);  //deallocate memory
}

如果您想使用 declare reduction，您需要为 arb_mat_t 创建一个 C++ 包装器。使用 declare reduction 让 OpenMP 决定如何进行缩减。但我非常怀疑您会发现这样的情况会比手动情况提供更好的性能。

使用 GMP/ARB 矩阵减少 OpenMP

Reduction in OpenMP with GMP/ARB matrices

c++

openmp

gmp

reduction