2
$\begingroup$

I know that is not the first time someone asks this question but I'm really confused.I'm new to MPI, and I tried to implement the Jacobi solver for a linear system $Ax=b$. I want to compare the time between parallel and serial problem, but for large $N$ my MPI code is much slower. I don't know if the problem lies in my code or on my computer. By the way, I am using 2 processes on a dual-core CPU. Any help will be useful!

That's my parallel code:`

#include <iostream>
#include "mpi.h"
#include <cmath>
#include <cstdlib>
using namespace std;

#define n 10000

void CreateMatrix(double *A){

     double value;
     double sum=0;
     int j=0;
     for (int i=0;i<n*n;i++){
             j++;
            //cout<<"Enter matrix value row i="<<i<<endl;
            //cin>>value;
            value=(rand()%1000)+1;

            A[i]=value;
     }

     for(int i=0;i<n;i++){
         A[i*n+i]=(1000*n+128.5)+rand()%8984;//dummy valus in diagonal...
     }                                       //in order to be diagonally dominant   
};

void CreateVector(double *b){


    double value;

        for(int i=0;i<n;i++){
            //cout<<"Enter vector b value row i="<<i<<endl;
            value=(rand()%974)-i*(rand()%20)+(rand()%1000);//random values
            b[i]=value;
        }
};

int JacobiMethod(double *A,double *b,double *x_out,int size,int rank,double tol,int max_iter){

    int n_local=n/size;
    double *xold;
    double *xnew;
    double *temp1=new double [n];
    double *temp2=new double [n];
    double *swap;
    int iterations=0;
    int k;
    double sum,sum1;
    double error=1;
    int flag;


    MPI_Allgather(b,n_local,MPI_DOUBLE,temp1,n_local,MPI_DOUBLE,MPI_COMM_WORLD);//gather all b matrices from each process to temp1
    xold=temp2;
    xnew=temp1;

    while(abs(error)>tol && iterations<max_iter){
        swap=xnew;
        xnew=xold;
        xold=swap;
        iterations+=1;

        for(int i=0;i<n_local;i++){

            k=i+rank*n_local;//tells us the location of diagonal element


            sum=0;
            for(int j=0;j<n;j++){
                if (j!=k){
                    sum+=A[j+n*i]*xold[j];//A is a continuus array so n*i gives us the next row

                }
            }

            x_out[i]=(b[i]-sum)/A[n*i+k];

        }
        MPI_Allgather(x_out,n_local,MPI_DOUBLE,xnew,n_local,MPI_DOUBLE,MPI_COMM_WORLD);
        sum1=0;
        for(int i=0;i<n;i++){
            sum1=sum1+pow((xnew[i]-xold[i]),2);
        }
        error=sqrt(sum1);
    }


    if (iterations>=max_iter){
        flag= -1;
    }
    else{
        flag= 0;

    }
    delete []A;
    delete []b;
    delete []temp1;
    delete  [] temp2;
    return flag;
};
void printResults( double *x_out,int size,int rank){
    int n_local=n/size;
    double *answ=new double [n];
    MPI_Gather(x_out, n_local, MPI_DOUBLE, answ, n_local, MPI_DOUBLE,0,MPI_COMM_WORLD);//gather all data to answ-->in one process(process 0)
    if (rank==0){
        cout<<"The algorith converges"<<endl;
        cout<<"The results are: "<<endl;
        for(int i=0;i<n;i++){
            //cout<<answ[i]<<endl;
        }
    }
    delete[] answ;
};

int main(){


double tol;//tolerance
int max_iter;

int converged;
//Mpi initialization
MPI_Init(NULL,NULL);
double *b_local;
double *A_local;
int size;
MPI_Comm_size(MPI_COMM_WORLD, &size);   
int rank;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);    
if (rank==0){
    b_local=new double[n];
    A_local=new double[n*n];
    cout<<"Enter talerance,number of iterations"<<endl;
    cin>>tol;
    cin>>max_iter;
    //Create A and scatter it to all process
    CreateMatrix(A_local);
    //Create b and scatter it to all process    
    CreateVector(b_local);

}
//data init
double *A=new double[n*n/size];
double *b=new double[n/size];
double *x_out=new double[n/size];
//brocast tol,max_iter to all processes
MPI_Bcast(&tol,1,MPI_DOUBLE,0,MPI_COMM_WORLD);//send to all processes
MPI_Bcast(&max_iter,1,MPI_INT,0,MPI_COMM_WORLD);    
//scatter vector b.Each process takes n/size
MPI_Scatter(b_local,n/size,MPI_DOUBLE,b,n/size,MPI_DOUBLE,0,MPI_COMM_WORLD);//here n_local cause we have only one column
//scatter it to all processes
MPI_Scatter(A_local,(n/size)*n,MPI_DOUBLE,A,(n/size)*n,MPI_DOUBLE,0,MPI_COMM_WORLD);//n_local*n--->number of elements in n/size rows
if (rank==0){
 delete [] b_local;
 delete [] A_local;
}
double time0,time1;
time0=MPI_Wtime();
converged=JacobiMethod(A,b,x_out,size,rank,tol,max_iter);
time1=MPI_Wtime();

if (converged==0){
    cout<<"Time needed for the jacobi algorith to be executed is :"<<time1-time0<<endl;
    printResults(x_out,size,rank);
}
else{
    cout<<"Jacobi doesnt converge"<<endl;
}
MPI_Finalize(); 
return 0;   

}
$\endgroup$
6
$\begingroup$

The first thing you need to ask yourself: is your problem big enough that the overhead of MPI messaging is less than the work that you save. Your problem size is 10k which is small, but on the other hand you seem to have a dense matrix, so there is a good amount of work.

Next: parallelizing a numerical method may change the mathematics. It looks like you're calculating the location of the diagonal element correctly, but it wouldn't hurt to print out residual norms: they should be the same between sequential and parallel. This is not true in general.

The only thing I see that is really wrong is that you use Allgather. This goes against the philosophy of MPI, which is that everything should be distributed. You should never have completely replicated data.

Sorry, these are just some thoughts. Maybe the issue is something else entirely. Start by checking that sequential and parallel run in the same number of iterations.

$\endgroup$
  • $\begingroup$ You are right about Allgather.I used it only due to Jacobi's cause xn+i needs xn (data from other processes) in order to be calculated.I couldn't find another so if you have to suggest something I will appreciate it.The difference between 2 programs are very big using mpirun -np 2 it takes 7 seconds for the programm to be executed and 0.8 for the serial one $\endgroup$ – Spyros Arvanitis Jun 27 at 19:37
  • $\begingroup$ By the way using static arrays instead of dynamics makes the code efficient.I mean the parallel program is faster than the serial one $\endgroup$ – Spyros Arvanitis Jun 27 at 22:55
  • $\begingroup$ @SpyrosArvanitis Are you timing outside the mpirun call or are you using a timer inside your program? Mpirun creates an ssh tunnel, so it takes a lot of time, relatively speaking. $\endgroup$ – Victor Eijkhout Jul 1 at 17:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.