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

Subject: Re: [Boost-users] Matrix Vector Multiplication in MPI (Coordinate Storage Format)
From: oswin krause (oswin.krause_at_[hidden])
Date: 2013-12-30 12:29:15

Hi,

1. using a sparse format for a tri-diagonal matrix is suboptimal as you 
know in advance which elements are the nonzero entries of the matrix. 
Use a dense-tridiagonal format instead. For example: take the rows of A 
and store the 1-3 nonzero elements in memory.

i.e. matrix

a11 a12 0
a21 a22 a23
0     a32 a33

is stored as

a11 a12 |a21 a22 a23 |a32 a33

than you only need the row index to compute the start of the row in the 
array.

2. using this format, parallelizing b=Ax is a breeze as you can compute 
every result value of b independently using a row of A.

Best,
Oswin

On 29.12.2013 12:45, ozhan fenerci wrote:
>
> 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;|
>
>
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users