//	matrix A
//
//	2	4	3  
//	4	2	1
//	3	1	6
//
//	matrix B
//
//	2	1	3
//	4	5	6
//	7	8	9
//	
//  ZERO-BASED INDEXING
//	
//	a = {2 4 3 4 2 1 3 1 6}
//	columns= {0 1 2 0 1 2 0 1 2}
//  idexRow = {0 3 6 9}
//
//	b = {2 1 3 4 5 6 7 8 9} (row order array)
//
//	We print the array in row-major order
//
//	ONE-BASED INDEXING
//	
//	a = {2 4 3 4 2 1 3 1 6}
//	columns={1 2 3 1 2 3 1 2 3}
//  indexRow = {0 3 6 9}
//
//	b = {2 4 7 1 5 8 3 6 9} (column order array)
//  
//	We print the array in column-major order (because the resoult is in column major order, ie transposed)
//
//
//

#include <stdio.h>
#include "mkl_types.h"
#include "mkl_spblas.h"

int main() 
{

#define M 3        
#define NNZ 9 
#define N 3 

		MKL_INT		m = M, nnz = NNZ, n=N;
        float		a[NNZ]	  = {2.0,4.0,3.0,4.0,2.0,1.0,3.0,1.0,6.0};
		MKL_INT		columns[NNZ]  = {0,1,2,0,1,2,0,1,2};
		MKL_INT		rowIndex[M+1] = {0,3,6,9};

		float		b[M][N]	= {2.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
		float		c[M][N]	= {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

        float		alpha = 1.0, beta = 0.0;
        MKL_INT		i, j;
		char		transa;
		char		matdescra[6];

		float		a1[NNZ]	  = {2.0,4.0,3.0,4.0,2.0,1.0,3.0,1.0,6.0};
		MKL_INT		columns1[NNZ]  = {1,2,3,1,2,3,1,2,3};
		MKL_INT		rowIndex1[M+1] = {1,4,7,10};
		float		b1[M][N]	= {2.0, 4.0, 7.0, 1.0, 5.0, 8.0, 3.0, 6.0, 9.0};
		float		c1[M][N]	= {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

//********************************
		//ZERO-BASED INDEXING
//********************************
        
		transa = 'n';
		matdescra[0] = 's';
		matdescra[1] = 'l';
		matdescra[2] = 'n';
		matdescra[3] = 'c';

		mkl_scsrmm(&transa, &m, &n, &m, &alpha, matdescra, a, columns, rowIndex, &(rowIndex[1]), &(b[0][0]), &n,  &beta, &(c[0][0]), &n);

		printf("                             \n");
        printf("   Right Solution: ZERO-BASED: C    \n");
		for (i = 0; i < m; i++)
		{
			for (j = 0; j < n; j++) {
				printf("%7.1f", c[i][j]);
			};
			printf("\n");
		};

		printf("                             \n");
        printf("   ZERO-BASED: C'    \n");
		for (i = 0; i < m; i++) 
		{
			for (j = 0; j < n; j++) 
			{
				printf("%7.1f", c[j][i]);
			};
			printf("\n");
		};

//********************************
		//ONE-BASED INDEXING
//********************************
        
		matdescra[3] = 'f';

		mkl_scsrmm(&transa, &m, &n, &m, &alpha, matdescra, a1, columns1, rowIndex1, &(rowIndex1[1]), &(b[0][0]), &n,  &beta, &(c1[0][0]), &n);

		printf("                             \n");
        printf("   ONE-BASED: C    \n");
		for (i = 0; i < m; i++) 
		{
			for (j = 0; j < n; j++) 
			{
				printf("%7.1f", c1[i][j]);
			};
			printf("\n");
		};

		printf("                             \n");
        printf("   ONE-BASED: C'    \n");
		for (i = 0; i < m; i++) 
		{
			for (j = 0; j < n; j++) 
			{
				printf("%7.1f", c1[j][i]);
			};
			printf("\n");
		};
		return 0;
}