Решение систем линейных алгебраических уравнений методом Гаусса с помощью различных алгоритмов, страница 2

//Получение матрицы от последнего компьютера

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

MPI_Recv(A[j],(n+1),MPI_DOUBLE,size-1,tag,ring,&st);

if(rank==0)

starttime=time(NULL);

}

V=(double*)malloc((n+1)*sizeof(double));

//Приведение матрицы к верхнему треугольному виду

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

{

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

{   

MPI_Recv(V,n+1,MPI_DOUBLE,k,tag,ring,&st);

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

{

if (abs(V[k*m+l])<=eps) kill_all(ring);

d=A[i][k*m+l]/V[k*m+l];

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

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

} 

}

}

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

{

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

MPI_Send(A[l],n+1,MPI_DOUBLE,k,tag,ring);

for(i=l+1;i<m;i++)

{

if (abs(A[l][rank*m+l])<=eps) kill_all(ring);

d=A[i][rank*m+l]/A[l][rank*m+l];

for(j=rank*m+l;j<n+1;j++)

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

}

}

// Проверка...

/*  for(i=0;i<m;i++)

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

printf("\nRank=%d,  A[%2d,%2d]=%lf",rank,i,j,A[i][j]);

*/   

//----------------------------------------------------------//Обратный ход 

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

MPI_Recv(&V[k*m],m,MPI_DOUBLE,k,tag,ring,&st);

for(k=m-1;k>-1;k--)

{

s=0;

j=rank*m+k;

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

s+=A[k][i]*V[i];

V[j]=(A[k][n]-s)/A[k][j];

}

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

MPI_Send(&V[m*rank],m,MPI_DOUBLE,k,tag,ring);

//  if (rank!=0) MPI_Send(&flag,1,MPI_INT,0,tag,ring);

//Вывод результата

if(rank==0)

{

finishtime=time(NULL);

/*  for(k=1;k<size;k++)

{

MPI_Recv(&recv_flag,1,MPI_INT,k,tag,ring,&st);

flag= recv_flag || flag;

}*/

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

/*  if (flag==1) fprintf(fp,"Решение не найдено.");

else

{*/

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

fprintf(fp,"%lf ",V[i]);

fprintf(fp,"\nВремя(в секундах)=%f",difftime(finishtime,starttime));

//  }

}

//Освободить память перед завершением...

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

free(A[i]);

free(A);

free(V);

MPI_Finalize();

return(0);

}

Второй алгоритм

#include <mpi.h>

#include <stdio.h>

#include <time.h>

#define ND 1

//#define N 4  //Number of rows in

//#define M 16 //Number of columns in

double **A, *V;

int n,m;

time_t starttime,finishtime;

main (int argc, char** argv)

{

int tag=11,rank,size;

MPI_Status st; MPI_Comm ring;

int dims[ND],per[ND],reord=0,prev,next;

double d,s;

int i,j,k,l,count0,count1;

FILE *fp;

MPI_Init(&argc,&argv);

MPI_Comm_rank(MPI_COMM_WORLD,&rank); //computer number

MPI_Comm_size(MPI_COMM_WORLD,&size); //total number of computers

per[0]=1; //ring

dims[0]=0; //for auto creation

MPI_Dims_create(size,ND,dims);

MPI_Cart_create(MPI_COMM_WORLD,ND,dims,per,reord,&ring); //создание топологии

//Read data from file (last comp. do it)

if (rank == (size-1))

{ //Read matrix dimension and sent it to others

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

fscanf(fp,"%d",&n);

for (i=0;i<size-1;i++)

MPI_Send(&n,1,MPI_INT,i,tag,ring);

m=n/size; //Calculate the width

//Выделение памяти под массив (A U b) & V

A=(double**)malloc(m*sizeof(double*));

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

A[i]=(double*)malloc((n+1)*sizeof(double));

V=(double*)malloc((n+1)*sizeof(double));

//Считывание матрицы А и вектора В из файла

count0=0;

count1=0;

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

{

for(k=0;k<n+1;k++)

fscanf(fp,"%lf",&V[k]);

if (count0 == (size-1))

{

for (i=0;i<n+1;i++)

A[count1][i]=V[i];

count1++;

count0=0;

}

else

{

MPI_Send(V,(n+1),MPI_DOUBLE,count0,tag,ring);

count0++;

}

}           

fclose(fp);

}

else

{

//Получить размерность от последнего

MPI_Recv(&n,1,MPI_INT,size-1,tag,ring,&st);

//Выделить память под А и В

m=n/size;

A=(double**)malloc(m*sizeof(double));

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

A[i]=(double*)malloc((n+1)*sizeof(double));

//Получение матрицы от последнего компьютера

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

MPI_Recv(A[j],(n+1),MPI_DOUBLE,size-1,tag,ring,&st);

if(rank==0)

starttime=time(NULL);

V=(double*)malloc((n+1)*sizeof(double));

}

//Приведение матрицы к верхнему треугольному виду

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

{

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

{

MPI_Recv(V,n+1,MPI_DOUBLE,k,k,ring,&st);

for (i=l;i<m;i++) //цикл зануления полосы

{ d=A[i][k+l*size]/V[k+l*size];

for(j=k+l*size;j<n+1;j++)

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

}

}

for(i=l+1;i<m;i++)

{d=A[i][rank+l*size]/A[l][rank+l*size];

for(j=rank+l*size;j<n+1;j++)

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

}

for (k=0;k<size-1;k++)

MPI_Send(A[l],n+1,MPI_DOUBLE,(rank+k+1)%size,rank,ring);

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

{ 

MPI_Recv(V,n+1,MPI_DOUBLE,k,k,ring,&st);

for(i=l+1;i<m;i++)

{ d=A[i][k+l*size]/V[k+l*size];

for(j=k+l*size;j<n+1;j++)

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

} 

}

}     

// Проверка...

//  for(i=0;i<m;i++)

//   for(j=0;j<n+1;j++)

//    printf("\nRank=%d,  A[%2d,%2d]=%lf",rank,i,j,A[i][j]);

//----------------------------------------------------------//Обратный ход 

if(rank!=size-1)

{

MPI_Recv(&V[(m-1)*size+rank+1],size-rank-1,MPI_DOUBLE,rank+1,tag,ring,&st);

}   

s=0;

j=(m-1)*size+rank;

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

s+=A[m-1][i]*V[i];

V[j]=(A[m-1][n]-s)/A[m-1][j];

MPI_Send(&V[(m-1)*size+rank], size-rank, MPI_DOUBLE, (rank+size-1)%size, tag, ring);

for(l=m-2;l>=0;l--)

{

MPI_Recv(&V[l*size+rank+1],n-(l*size+rank+1), MPI_DOUBLE, (rank+1)%size, tag, ring, &st);

s=0;

j=l*size+rank;

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

s+=A[l][i]*V[i];

V[j]=(A[l][n]-s)/A[l][j];  

if((rank!=0)||(l!=0))

MPI_Send(&V[l*size+rank],n-(l*size+rank), MPI_DOUBLE, (rank+size-1)%size, tag, ring);

}

//Вывод результата

if(rank==0)

{

finishtime=time(NULL);

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

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

fprintf(fp,"%lf ",V[i]);

fprintf(fp,"\nВремя(в секундах)=%f",difftime(finishtime,starttime));

}

//Освободить память перед завершением...

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

free(A[i]);

free(A);

free(V);

MPI_Finalize();

return(0);

}