tridiagonal solver in oneMKL for dpc++ (linear solver routine getrs())

[Oluwatosin]

Hello,

I have been trying to use the getrs() routine, but I can't seem to get pass, size = 50 when executing the code, I need to be able to measure for size close to 100,000,000. I have to call first the getrf() routine according to the documentation, and I think this is where I miss it. I would appreciate any suggestions to improve the getrf() section of program; thank you. 

I am also following the matrix_multiplication with gemm code sample to generate my code;

The basic idea of what I am working on is:  
1. Create a matrix, A, of size n * n;
2. Initialize the matrix A, to zero;
3. Populate the sub, super and main diagonals with -1.0, -1,0 and 2.0  or random numbers;
4. Create a vector, for the right hand side and initialize it randomly.
4. Call the getrf() routine to factorize the system (device code);
5. Call the getrs() routine to solve the system (device code);
6. Compare the results of getrs() with results of host tridiagonal solver algorithm.

Here's what I have for the first four steps (getrf() )

int main(int argc, char* argv[])

{

size_t size = atoi(argv[1]);

sycl::queue device_queue{sycl::default_selector{}};

sycl::context context  = device_queue.get_context();

sycl::device device = device_queue.get_device();

int_64 i , j;

int64_t m = size, n = size, nrhs = size, lda = size;

// Allocate memory
double **A = sycl::malloc_shared<double*>(m * n, device_queue);

for(i = 0; i < lda; i++){

A[i] = sycl::malloc_shared<double>(lda, device_queue);}

// Initialize array to zero;
for(i = 0; i < m; i++){

for(j = 0; j < n; j++){
A[i][j] = 0.0;}}

// Populate sub,super, main diagonals of A

for(i = 0; i < lda; i++){

A[i][i] = 2.0; // Main diagonal

if(i < lda){ 

A[i][i+1] = -1.0;} // Super diagonal

if(i > 0){

A[i][i-1] = -1.0;} // Sub diagonal

}

// Other arrays

int64_t *ipiv = sycl::malloc_shared<std::int64_t>(ldb, device_queue);

auto b = sycl::malloc_shared<double>(ldb, device_queue);

//Initialize result array, b

for(i = 0; i < ldb; i++){

b[i] = 1.0;}

int64_t info = 0; 

try

{

int64_t scratchpad_size_getrf = oneapi::mkl::lapack::getrf_scratchpad_size<double>(device_queue, m, n, lda);

double* scratchpad_getrf = sycl::malloc_shared<double>(scratchpad_size_getrf, device, context);

auto event1 = oneapi::mkl::lapack::getrf(device_queue, m, n, *A, lda, ipiv, scratchpad_getrf, scratchpad_getrf, scratchpad_size_getrf);

event1.wait_and_throw();

cl::sycl::fre(scartchpad_getrf, context);

} catch(oneapi::mkl::lapack::exception const& e1){

std::cout << ............

if(e1.info() > 0){

info = e1.info();

}

return info;}

}

I would appreciate any suggestions to improve this section of the program.Thank you!

[MRajesh_intel · Accepted solution]

Hi,

>>I would appreciate any suggestions to improve the getrf() section of program

Please refer to the sample(getri.cpp) present in $MKL_ROOT/examples/dpcpp/lapack/source folder.

Regards

Rajesh.

[Oluwatosin]

Hi, thanks. I am working on the DevCloud. Can I get an external link, like a git repository, to download this example? Thank you.

[MRajesh_intel]

Hi,

You can view the MKL samples in the dir:

/opt/intel/oneapi/mkl/latest/examples

You can unzip the examples_dpcpp.tgz and check the getri.cpp example for reference.

Regards

Rajesh.

[Oluwatosin]

Alright! Thanks!

[MRajesh_intel]

Hi,

Can you please let us know if your issue has been resolved?

Regards

Rajesh.

[Oluwatosin]

Hello,
Yes, thank you!

[MRajesh_intel]

Hi,

Thanks for the confirmation!

As this issue has been resolved, we will no longer respond to this thread. If you require any additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.

Have a Good day.

Regards

Rajesh