magma sparse solver has different answer compare then MKL

Open discussion for MAGMA library (Matrix Algebra on GPU and Multicore Architectures)

magma sparse solver has different answer compare then MKL

Postby 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.
alarm
 
Posts: 4
Joined: Wed Mar 01, 2017 9:12 pm

Re: magma sparse solver has different answer compare then MK

Postby 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
hartwig anzt
 
Posts: 75
Joined: Tue Sep 02, 2014 5:44 pm

Re: magma sparse solver has different answer compare then MK

Postby 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.
alarm
 
Posts: 4
Joined: Wed Mar 01, 2017 9:12 pm


Return to User discussion

Who is online

Users browsing this forum: Bing [Bot] and 2 guests

cron