Алгоритмы решения СЛАУ методом исключений Гаусса.

Условие задачи

Построить два параллельных алгоритма решения СЛАУ методом исключений Гаусса.  В первом случае матрица разрезается на горизонтальные ленты; во втором случае матрица распределяется по строкам – по очереди для всех компьютеров.

Исходные данные заданы в файле input.txt; формат: N –количество строк, далее коэффициенты строки и коэффициент правой части. Результат – в файл output.txt.

Ход работы

В процессе выполнения задания были написаны две программы, реализующие указанные алгоритмы с различным способом распределения исходной матрицы.

Программа 1 (Алгоритм 1)

#include <mpi.h>

#include <stdio.h>

#include <sys/time.h>

#include <stdlib.h>

float **A, *W;

int M,N,C;

main(int argc, char ** argv)

{

int tag,rank,size,ND,dims[1],per[1],soure,dest;

FILE *F;

MPI_Status st;

MPI_Request reg;

MPI_Comm ring;

struct timeval t1,t2;

int l,i,j,k;

int d,q;

ND=1;

per[0]=1;

MPI_Init(&argc,&argv);

MPI_Comm_rank(MPI_COMM_WORLD,&rank);

MPI_Comm_size(MPI_COMM_WORLD,&size);

dims[0]=size;

MPI_Dims_create(size,ND,dims);

MPI_Cart_create(MPI_COMM_WORLD,ND,dims,per,0,&ring);

//printf("Hallo! My rank = %d\n",rank);

//reading of matrix

if (rank==0)

{

F=fopen("input.txt","rt");

fscanf(F,"%d",&M);

MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);

N = M / size;

W = (float*)malloc((M+1)*sizeof(float));

A = (float**)malloc(N*sizeof(float*));

for (i=0; i<N; i++)

A[i]=(float*)malloc((M+1)*sizeof(float));

for (i=0; i<N; i++)

for (j=0; j<M+1; j++)

fscanf(F,"%f",&A[i][j]);

for (i=N; i<M; i++)

{

for (j=0; j<M+1; j++)

fscanf(F,"%f",&W[j]);

MPI_Send(W,M+1,MPI_FLOAT,(int)(i/N),10,MPI_COMM_WORLD);

//printf("0 sent line  %d to %d\n",i,i/N);

}

fclose(F);

//printf("matrix sent succ.\n");

//getchar();

}

else

{

MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);

N = M / size;

W = (float*)malloc((M+1)*sizeof(float));

A = (float**)malloc(N*sizeof(float*));

for (i=0; i<N; i++)

A[i]=(float*)malloc((M+1)*sizeof(float));

for (i=0; i<N; i++)

{

MPI_Recv(W,M+1,MPI_FLOAT,0,10,MPI_COMM_WORLD,&st);

//printf("%d got  line %d from 0\n",rank,i);

for (j=0; j<M+1; j++)

A[i][j]=W[j];

}

//printf("matrix gotten succ.\n");

}

// Direct Gauss moving

//Reading part

if (rank==0) gettimeofday(&t1,(struct timezone *)0);

d=0;

for (l=0; l<rank; l++)

{

for (i=0; i<N; i++)

{

//printf("Wait for %d line from %d\n",10*rank+i,l);

MPI_Recv(W,M+1,MPI_FLOAT,l,10*rank+i,MPI_COMM_WORLD,&st);

//for(k=0; k<M; k++)

// printf("%f ",W[k]);

//printf("\n");

//getchar();

for (j=0; j<N; j++)

{

for (k=d+1; k<M+1; k++)

A[j][k]=A[j][k]-W[k]*A[j][d];

A[j][d]=0;

}

d++;

}

}

// Part of self-transformation

for (i=0; i<N; i++)

{

if (A[i][d+i]==0)

{

printf("Error!\n");

break;

}

for(j=d+i+1; j<M+1; j++)

A[i][j]=A[i][j]/A[i][d+i];

A[i][d+i]=1;

for(j=i+1; j<N; j++)

{

for(k=d+i+1; k<M+1; k++)

A[j][k]=A[j][k]-A[j][d+i]*A[i][k];

A[j][d+i]=0;

}

}

// Sending Part

for(i=0; i<N; i++)

{

for(j=rank+1; j<size; j++)

{

MPI_Send(A[i],M+1,MPI_FLOAT,j,10*j+i,MPI_COMM_WORLD);

//printf("Sent %d line to %d\n",20+i,j);

}

}

//Indirect Gauss moving

if (rank!=size-1)

MPI_Recv(W,M,MPI_FLOAT,rank+1,5,MPI_COMM_WORLD,&st);

for (i=N-1; i>=0; i--)

{

d=rank*N+i;

W[d]=A[i][M];

for(j=M-1; j>d; j--)

{

W[d]=W[d]-W[j]*A[i][j];

}

W[d]=W[d]/A[i][d];

}

if (rank!=0)

MPI_Send(W,M,MPI_FLOAT,rank-1,5,MPI_COMM_WORLD);

if (rank==0) gettimeofday(&t2,(struct timezone *)0);

// OUTPUT

if (rank==0)

{

F=fopen("output.txt","wt");

for(i=0; i<M; i++)

fprintf(F,"%f\n",W[i]);

fprintf(F,"Time = %ld msec",(t2.tv_sec-t1.tv_sec)*1000000+(t2.tv_usec-t1.tv_usec));

fclose(F);

}

MPI_Finalize();

return(0);

}

Программа 2 (Алгоритм 2)

#include <mpi.h>

#include <stdio.h>

#include <sys/time.h>

#include <stdlib.h>

float **A, *W;

int M,N,C;

main(int argc, char ** argv)

{

int tag,rank,size,ND,dims[1],per[1],soure,dest;

FILE *F;

MPI_Status st;

MPI_Request reg;

MPI_Comm ring;

struct timeval t1,t2;

int l,i,j,k;

int d,q,z;

ND=1;

per[0]=1;

MPI_Init(&argc,&argv);

MPI_Comm_rank(MPI_COMM_WORLD,&rank);

MPI_Comm_size(MPI_COMM_WORLD,&size);

dims[0]=size;

MPI_Dims_create(size,ND,dims);

MPI_Cart_create(MPI_COMM_WORLD,ND,dims,per,0,&ring);

//printf("Hallo! My rank = %d\n",rank);