magma sparse solver has different answer compare then MKL

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

magma sparse solver has different answer compare then MKL

Post by alarm » Sun Mar 05, 2017 9:24 pm

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.
please help me, what's wrong in my magma sample code ?

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

Post by hartwig anzt » Mon Mar 06, 2017 7:56 pm

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

Post by alarm » Mon Mar 06, 2017 9:44 pm

Dear Hartwig Anzt:

Thanks your answer.
you are right.
I using other tools to verify the matrix.
get the same result. ^^
maybe I missing some parameter in MKL pardiso.

Have a good day.
AlarmChang.

Post Reply