## magma sparse solver has different answer compare then MKL

Open discussion for MAGMA library (Matrix Algebra on GPU and Multicore Architectures)
alarm
Posts: 4
Joined: Wed Mar 01, 2017 9:12 pm

### magma sparse solver has different answer compare then MKL

Hi All:

I am new in magma,
I copy a CSR matrix from mkl example, try to using magma sparse solver to solve the matrix X
fallowing is my code.

Code: Select all

``````// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

// includes, project
#include "magma_v2.h"
#include "magmasparse.h"
#include "testings.h"

/* ////////////////////////////////////////////////////////////////////////////
-- testing any solver
*/
int main(int argc, char** argv)
{

int n = 8;
int ia[9] = { 0, 4, 7, 9, 11, 14, 16, 17, 18 };//zero base
int ja[18] =	{ 0, 2, 5, 6,	1, 2, 4,	2, 7,	3, 6,	4, 5, 6,	5, 7,	6,	7	};//zero base
double a[18] =	{ 7.0, 1.0, 2.0, 7.0,	-4.0, 8.0, 2.0,	1.0, 5.0,	7.0, 9.0,	5.0, 1.0, 5.0,	-1.0, 5.0,	11.0,	5.0	};

double h_b[8] = { 1, 1, 1, 1, 1, 1, 1, 1 }, h_x[8] = {0, 0, 0, 0, 0, 0, 0, 0};

magma_int_t info = 0;
TESTING_CHECK(magma_init());
magma_print_environment();

magma_dopts zopts;
magma_queue_t queue;
magma_queue_create(0, &queue);

magma_d_matrix A = { Magma_CSR };
magma_d_matrix d_A = { Magma_CSR };
magma_d_matrix x = { Magma_CSR }, b = { Magma_CSR };
magma_d_matrix d_x = { Magma_CSR }, d_b = { Magma_CSR };

//zopts.solver_par.solver = Magma_GMRES;
//zopts.solver_par.solver = Magma_CG;
//zopts.solver_par.solver = Magma_PIDRMERGE;
zopts.solver_par.solver = Magma_BICGSTAB;
zopts.solver_par.restart = 8;
zopts.solver_par.maxiter = 1000;
zopts.solver_par.rtol = 1e-10;
zopts.solver_par.maxiter = 1000;
zopts.precond_par.solver = Magma_ILU;
//zopts.precond_par.solver = Magma_JACOBI;
zopts.precond_par.levels = 0;
zopts.precond_par.trisolver = Magma_CUSOLVE;

TESTING_CHECK(magma_dsolverinfo_init(&zopts.solver_par, &zopts.precond_par, queue));
magma_dcsrset(n, n, &ia[0], &ja[0], a, &A, queue);
magma_dvset(n, 1, h_b, &b, queue);
magma_dvset(n, 1, h_x, &x, queue);

// preconditioner

if (zopts.solver_par.solver != Magma_ITERREF) {
TESTING_CHECK(magma_d_precondsetup(A, b, &zopts.solver_par, &zopts.precond_par, queue));
}

//TESTING_CHECK(magma_dmconvert(A, &B, Magma_CSR, zopts.output_format, queue));

printf("\n%% matrix info: %lld-by-%lld with %lld nonzeros\n\n",
(long long)A.num_rows, (long long)A.num_cols, (long long)A.nnz);

printf("matrixinfo = [\n");
printf("%%   size   (m x n)     ||   nonzeros (nnz)   ||   nnz/m   ||   stored nnz\n");
printf("%%============================================================================%%\n");
printf("  %8lld  %8lld      %10lld             %4lld        %10lld\n",
(long long)A.num_rows, (long long)A.num_cols, (long long)A.true_nnz,
(long long)(A.true_nnz / A.num_rows), (long long)A.nnz);
printf("%%============================================================================%%\n");
printf("];\n");
/**/

// Copy the system to the device(optional, only necessary if using the GPU)
magma_dmtransfer(A, &d_A, Magma_CPU, Magma_DEV, queue);
magma_dmtransfer(b, &d_b, Magma_CPU, Magma_DEV, queue);
magma_dmtransfer(x, &d_x, Magma_CPU, Magma_DEV, queue);

// vectors and initial guess
//TESTING_CHECK(magma_dvinit(&d_A, Magma_DEV, A.num_rows, 1, one, queue));
//magma_dvinit( &x, Magma_DEV, A.num_cols, 1, one, queue );
//magma_d_spmv( one, dB, x, zero, b, queue );                 //  b = A x
//magma_dmfree(&x, queue );
//TESTING_CHECK(magma_dvinit(&x, Magma_DEV, A.num_cols, 1, zero, queue));

info = magma_d_solver(d_A, d_b, &d_x, &zopts, queue);
if (info != 0) {
printf("%%error: solver returned: %s (%lld).\n",
magma_strerror(info), (long long)info);
}
printf("convergence = [\n");
magma_dsolverinfo(&zopts.solver_par, &zopts.precond_par, queue);
printf("];\n\n");

zopts.solver_par.verbose = 0;
printf("solverinfo = [\n");
magma_dsolverinfo(&zopts.solver_par, &zopts.precond_par, queue);
printf("];\n\n");

printf("precondinfo = [\n");
printf("%%   setup  runtime\n");
printf("  %.6f  %.6f\n",
zopts.precond_par.setuptime, zopts.precond_par.runtime);
printf("];\n\n");

/*
for (int i = 0; i < 8; i++)
{
printf("b %d = %f\n", i, b.val[i]);
//printf("x %d = %f\n", i, x.val[i]);
}
*/

magma_dmtransfer(d_A, &A, Magma_DEV, Magma_CPU, queue);
magma_dmtransfer(d_b, &b, Magma_DEV, Magma_CPU, queue);
magma_dmtransfer(d_x, &x, Magma_DEV, Magma_CPU, queue);
//int m = 8;
//double *sol = (double*)calloc(m, sizeof(double));
//magma_dvget(x, &m, &n, &sol, queue);
//magma_dmconvert(x, h_x, Magma_CSR, zopts.output_format, queue);

for (int i = 0; i < 8; i++)
{
//printf("x %d = %f\n", i, sol[i]);
printf("x [%d] = %f\n", i, x.val[i]);
}

//magma_dmfree(&dB, queue);
//magma_dmfree(&B, queue);
magma_dmfree(&A, queue);
magma_dmfree(&d_A, queue);
magma_dmfree(&x, queue);
magma_dmfree(&d_x, queue);
magma_dmfree(&b, queue);
magma_dmfree(&d_b, queue);

magma_queue_destroy(queue);
TESTING_CHECK(magma_finalize());
return info;
}

``````
the output is:

Code: Select all

``````#magma AX=B, matrix x answer.
x [0] = 0.051948
x [1] = -0.195455
x [2] = 0.000000
x [3] = 0.025974
x [4] = 0.109091
x [5] = 0.000000
x [6] = 0.090909
x [7] = 0.200000
``````
I am using MKL pardiso to solve the same question, I got fallow answer

Code: Select all

``````Solve completed ...
The solution of the system is:
x [0] = -0.041860
x [1] = -0.003413
x [2] =  0.117250
x [3] = -0.112640
x [4] =  0.024172
x [5] = -0.107633
x [6] =  0.198720
x [7] =  0.190383
``````
the answer from magma seems very different compare with mkl pardiso.

Have a good day.
AlarmChang.

hartwig anzt
Posts: 90
Joined: Tue Sep 02, 2014 5:44 pm

### Re: magma sparse solver has different answer compare then MK

Dear AlarmChang,

maybe you use MKL in the wrong fashion? From the last row, you can see right away that the MKL computation has quite some rounding.... 1/5 = 0.2 the last time I checked...

Also, if I run the system through matlab, I get very much the same result, see below.

Let me know if you have further questions... Hartwig

M =

7 0 1 0 0 2 7 0
0 -4 8 0 2 0 0 0
0 0 1 0 0 0 0 5
0 0 0 7 0 0 9 0
0 0 0 0 5 1 5 0
0 0 0 0 0 -1 0 5
0 0 0 0 0 0 11 0
0 0 0 0 0 0 0 5

>> b = ones(8,1);
>> x = M\b

x =

0.0519
-0.1955
-0.0000
0.0260
0.1091
0.0000
0.0909
0.2000

alarm
Posts: 4
Joined: Wed Mar 01, 2017 9:12 pm

### Re: magma sparse solver has different answer compare then MK

Dear Hartwig Anzt: