[Boost-users] Matrix Vector Multiplication in MPI (Coordinate Storage Format)

Subject: [Boost-users] Matrix Vector Multiplication in MPI (Coordinate Storage Format)
From: ozhan fenerci (fenerci.ozhan_at_[hidden])
Date: 2013-12-29 06:45:29

I wrote a serial code for a matrix-vector multiplication in a Coordinate
Storage Format(COO). COO method has already made a big speed-up. My matrix
is in a tridiagonal form. I want to parallelize the code and see how the
parallelizm improves the code. My matrix could grow up to
2.000.000*2.000.000. Any efficient parallelization approach doesn't come to
my mind (I have written some basic mpi programs but not much more. So my
experience is limited.) Could you give me some headstart?

Thanks.

Özhan Fenerci

#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;