#include <stdio.h>
#include "math.h"
#include "mkl_blas.h"
#include "mkl_spblas.h"
#include "mkl.h"
#include "mkl_rci.h"

#define size 128

int main(void) {

	int *ia, *ja;
	double *A, *rhs, *Sol, *b;

	int i,N, NNZ; 
        N = 3314;
	NNZ = 20820;

	ia = (int *)calloc(N+1, sizeof(int));
	ja = (int *)calloc(NNZ, sizeof(int));
	A  = (double *)calloc(NNZ, sizeof(double));
	rhs = (double *)calloc(N, sizeof(double));
	Sol =(double *)calloc(N,sizeof(double));
	b = (double *)calloc(N, sizeof(double));

        ia[0]=1;

	// Read values from the document
	FILE *fp1, *fp2;

	fp1 = fopen("ia.txt","r");
	fp2 = fopen("ja_A.txt","r");
        if((!fp1 )|| (!fp2)) {
	   printf("Can't open txt files!");
	   return -1;
	}
        for(i =0; i<N; i++) 
		fscanf(fp1, "%d %lf\n", &ia[i+1], &b[i]);
	for(i=0; i<NNZ; i++)
		fscanf(fp2, "%d %lf\n", &ja[i], &A[i]);

	fclose(fp1);
	fclose(fp2);

	for(i=N-40; i<N; i++)
		printf("%d %f %d %f\n ", ia[i+1], b[i],ja[i], A[i]);


        int ivar = N;
	char cvar ='N';
	for (i=0; i<N; i++) Sol[i] = 1.0;
        printf("dcsrgemv\n");
	mkl_dcsrgemv(&cvar,&ivar, A, ia, ja, Sol, rhs); 
        
	// call the linear solver 
        printf("Linear Solver begins\n");
	LinearSolver(&N,&NNZ, ia, ja, A, rhs );
        
	double nrm = 0.0;

	for(i=0;i<N; i++) nrm += (rhs[i]-Sol[i]) *(rhs[i]-Sol[i]);
	
	printf("The error is %d\n", sqrt(nrm));

	free(ia);
	free(ja);
	free(A);
	free(rhs);
	free(Sol);
        free(b);
	return 0;
}