#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <vector>
#include <iostream> 
#include "mpi.h"
#include <stdio.h>


int main(int argc, char *argv[] ) {

using namespace boost::numeric::ublas;
using std::cout;
using std::endl;

int numprocs, rank, chunk_size, i,j;
int max, mymax,rem;
int totalNonZeros; //totalNonZeros of TriDiagonal Matrix
MPI_Status status;
//b=A*x
banded_matrix<double> A (10, 10, 1, 1);
//vector<double> x(10); 
//vector<double> b(10); 
int* x = new int[10];
int* b = new int[10];    


/* Initialize MPI */
MPI_Init( &argc,&argv);
MPI_Comm_rank( MPI_COMM_WORLD, &rank);
MPI_Comm_size( MPI_COMM_WORLD, &numprocs);  
//MPI_Status status; 

if (rank == 0)  {
for (signed i = 0; i < signed (A.size1 ()); ++ i) {
    for (signed j = std::max (i - 1, 0); j < std::min (i + 2, signed (A.size2 ())); ++j)
        A (i, j) = i + j;
    x[i]= 1;
}       
}   


if(rank == 0) {
cout << A <<endl<<endl;


totalNonZeros = 0;
for (signed i = 0; i < signed (A.size1 ()); ++ i) {
    for (signed j = std::max (i - 1, 0); j < std::min (i + 2, signed (A.size2 ())); ++ j)
        if( A(i,j) !=0) 
            totalNonZeros++;
}
}


cout<<"Number ot NonZeros "<<totalNonZeros<<endl;


if(rank == 0) {

     int * row_idx = new int[totalNonZeros]; // row index array
     int * col_idx = new int[totalNonZeros]; // col index array
     int * val = new int[totalNonZeros]; //val array    
     int index=0;

for (signed i = 0; i < signed (A.size1 ()); ++ i) {
    for (signed j = std::max (i - 1, 0); j < std::min (i + 2, signed (A.size2 ())); ++ j)
        if( A(i,j) !=0)   {
            row_idx[index] = i; 
            col_idx[index]=j;
            val[index] = A(i,j);  
            index++;              
        }
}
//Clear out result
for ( int k=0; k<A.size1(); ++k)
{
    b[k]= 0;
}
//Matrix-vector Multiply
for (int index=0 ; index <totalNonZeros;++index){
     b[row_idx[index]] = b[row_idx[index]]+ val[index]*x[col_idx[index]]; cout<<b[3]<<endl;
} 

//delete row_idx;
// delete col_idx;
//delete val;
}                       

MPI_Finalize();
return 0;