Я пытаюсь вызвать трехзначный решающий код Чжана из openacc.Я делаю следующее:
Я помещаю его код в отдельный файл и компилирую с именем pcr.cu
#include <cstdio>
#include <cuda_runtime.h>
__device__ void Solve_Kernel_PCR(float * alist, float * blist, float * clist, float * dlist, float * xlist, int iter_max, int DMax)
{
int idx_row = blockIdx.x*blockDim.x + threadIdx.x;
int row_max = DMax - 1;
// printf("idx_row = %d iter_max= %d\n",idx_row, iter_max);
int stride = 1;
int next_stride = stride;
float a1, b1, c1, d1;
float k01, k21, c01, a21, d01, d21;
bool next_or_ot = true;
int accum;
for (int iter = 0; iter < iter_max; iter++)
{
if ( next_or_ot ) {
next_stride = stride<<1;
// 1 for updating 'a'
if ((idx_row - stride)<0) {
// 1.1 if it is the 'first' line
a1 = 0.0f;
k01 = 0.0f;
c01 = 0.0f;
d01 = 0.0f;
} else if ((idx_row - next_stride)<0) {
// 1.2 if no place for 'a'
a1 = 0.0f;
k01 = alist[idx_row]/blist[idx_row - stride];
c01 = clist[idx_row - stride]*k01;
d01 = dlist[idx_row - stride]*k01;
} else {
// 1.3 for rest general rows
k01 = alist[idx_row]/blist[idx_row - stride];
a1 = -alist[idx_row - stride]*k01;
c01 = clist[idx_row - stride]*k01;
d01 = dlist[idx_row - stride]*k01;
}
// 2 for updating 'c'
if ((idx_row + stride)>row_max) {
// 2.1 if it is the 'last' line
c1 = 0.0f;
k21 = 0.0f;
a21 = 0.0f;
d21 = 0.0f;
} else if ((idx_row + next_stride)>row_max) {
c1 = 0.0f;
k21 = clist[idx_row]/blist[idx_row + stride];
a21 = alist[idx_row + stride]*k21;
d21 = dlist[idx_row + stride]*k21;
} else {
k21 = clist[idx_row]/blist[idx_row + stride];
c1 = -clist[idx_row + stride]*k21;
a21 = alist[idx_row + stride]*k21;
d21 = dlist[idx_row + stride]*k21;
}
// 3 for updating 'b'
b1 = blist[idx_row] - c01 - a21;
// 4 for updating 'd'
d1 = dlist[idx_row] - d01 - d21;
stride = next_stride;
int pos = idx_row-2*stride;
accum = 0;
for ( size_t iter = 0; iter<5; iter++ ) {
if (pos >=0 && pos < DMax) accum++;
pos+=stride;
}
if (accum < 3) {
next_or_ot = false;//Turn of for ever
}
}
__syncthreads();
alist[idx_row] = a1;
blist[idx_row] = b1;
clist[idx_row] = c1;
dlist[idx_row] = d1;
}
if ( accum==1 ) {
xlist[idx_row] = dlist[idx_row] / blist[idx_row];
} else if ( (idx_row-stride)<0 ) {
int i = idx_row; int k = idx_row+stride;
float f = clist[i]/blist[k];
xlist[i] = (dlist[i]-dlist[k]*f)/(blist[i]-alist[k]*f);
} else {
int i = idx_row - stride; int k = idx_row;
float f = alist[k]/blist[i];
xlist[k] = (dlist[k]-dlist[i]*f)/(blist[k]-clist[i]*f);
}
}
и файл main.cpp:
#pragma acc routine gang bind("_Z16Solve_Kernel_PCRPfS_S_S_S_ii")
__device__ void Solve_Kernel_PCR(float * alist, float * blist, float * clist, float * dlist, float * xlist, int iter_max, int DMax);
int main(int argc, char *argsv[]) {
size_t diagonal_size ;
diagonal_size=atoi(argsv[1]);
float *alist = (float *)malloc(sizeof(float) * diagonal_size);
float *blist = (float *)malloc(sizeof(float) * diagonal_size);
float *clist = (float *)malloc(sizeof(float) * diagonal_size);
float *dlist = (float *)malloc(sizeof(float) * diagonal_size);
float *xlist = (float *)malloc(sizeof(float) * diagonal_size);
float delx=1./(diagonal_size-1);
for (int i = 0; i < diagonal_size; i++) {
alist[i] = 1.0f;
blist[i] = -2.0f;
clist[i] = 1.0f;
// dlist[i] = -1.; // rand() % 100 + 1;
xlist[i] = 0.0f;
}
float pi = atan(1.0) * 4.0;
for (int i = 0; i < diagonal_size; i++) {
dlist[i] = -pi * pi * sin(OMEGA * i * delx * pi) * delx * delx;
}
alist[0] = 0.0;
clist[diagonal_size - 1] = 0.0;
int DMax = diagonal_size;
int iter_max=count_iter(DMax);
printf("iter_max= %d\n", count_iter(DMax) );
int ngang=1;
int N=diagonal_size;
#pragma acc data copy(alist[0:diagonal_size],blist[0:diagonal_size],clist[0:diagonal_size],dlist[0:diagonal_size] ,xlist[0:diagonal_size] )
{
#pragma acc parallel num_gangs(ngang)
Solve_Kernel_PCR(alist, blist, clist, dlist, xlist, iter_max, DMax);
}
for (size_t it = 0; it < diagonal_size; it++) {
// std::cout << alist[it] << " " << blist[it] << " " << clist[it] << " " <<
// xlist[it] << " " << dlist[it] << std::endl;
//printf("%f \n", xlist[it]);
printf("%f \n", dlist[it]);
}
float err0 = 0.0;
float err1 = 0.0;
for (int i = 1; i < diagonal_size-1; i++) {
// printf("r %lf %lf %lf %lf\n ",dl[i], d[i],du[i], r[i]);
err1 = fabs(dlist[i] - sin(OMEGA * i * delx * pi));
if (err0 < err1) {
err0 = err1;
}
}
printf(" l infinity of Error = %lf \n",err0 );
}
Это простой контрольный пример, так как я хотел бы вызвать эту функцию из подпрограммы openacc gang.
, компилируя код cuda с помощью "nvcc -rdc true -arch = sm_60" и связывая с pgc ++, он связывается нормально.
Однако результаты неверны для входного значения> 32
Я предполагаю, что это из-за длины вектора по умолчанию 32. Мне нужно иметь возможность обрабатывать сетки до 1023. Есть ли какое-либо решение этой проблемы?