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@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users