用GMP/ARB矩阵约简OpenMP

Question

我想并行一个程序，我已经写了，它计算了一系列涉及矩阵和向量积，其结果是一个向量。由于争论变得非常小和大，我使用ARB (基于GMP，MPFR和打火机)，以防止失去的意义。此外，由于序列元素是独立的，矩阵维度并不大，但是需要计算高达50k的元素，所以让多个线程计算序列中的几个元素是有意义的，也就是说，每个线程可以并行计算10k元素，然后将结果向量相加。

现在的问题是，将向量和矩阵相加的ARB函数不是一个标准操作，可以很容易地与openmp还原一起使用。当我天真地尝试编写自定义约简时，g++抱怨void类型，因为ARB中的操作没有返回值：

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

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

例如:类似问题的随机循环如下(当然，我的实际程序是不同的)。

[...]

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        
}

[...]

当然，在使用多个线程时会出现故障，因为我没有在RRes上指定裁减。这里，Qmvd()和Qmvd()是自写矩阵-矩阵和矩阵向量乘积函数，RMa、RMb和RMV是随机矩阵和向量。

现在的想法是减少RRes，这样每个线程都可以计算RRes的私有版本，包括最终结果的一小部分，然后再使用arb_mat_add将它们全部添加起来。我可以编写一个函数矩阵(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。)

我的问题非常类似于这一个。

编辑:对不起，这是我尝试的一个最小的，完整的和可验证的例子。其他必需的包是arb、弗林特、gmp、mpfr (可通过包管理器获得)和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

你可以像这样用手做减价

#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决定如何进行缩减。但我非常怀疑你会发现这样的情况比手动情况下的性能更好。

Answer 2

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

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