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

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

MPI_Dims_create(size,ND,dims);

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

MPI_Cart_shift(ring,0,-1,&prev,&next); //определение соседей

//Read data from file (last computer do it)

if (rank == (size-1))

{ //Read matrix dimension and sent it to others

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

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

tau=1.5/n;

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

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

L=n/size+1; //Calculate the width

/* Начальное количество элементов у вектора Y0 для каждого компьютера

Начальное значение указателя откуда нужно начинать умножение */

start=size-(L*size-n);

//Выделение памяти под матрицу A и вектор b

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

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

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

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

//для rank<start

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

{

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

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

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

if (i!=size-1) 

{

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

MPI_Send(A[j],(n+1),MPI_DOUBLE,i,tag,ring);

}

}

//для rank>start  

for(i=start;i<size;i++)

{

for (j=0;j<L-1;j++)

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

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

if (i!=size-1) 

{

for(j=0;j<L-1;j++)

MPI_Send(A[j],(n+1),MPI_DOUBLE,i,tag,ring);

}

}

fclose(fp);

}

else

{

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

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

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

L=n/size+1;

tau=1.5/n;

start=size-(L*size-n);

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

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

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

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

if(rank<start)

{

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

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

}

else

{

for(j=0;j<L-1;j++)

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

}

if(rank==0)

starttime=time(NULL);

}

Y0=(double*)malloc((L+1)*sizeof(double));

Y1=(double*)malloc((L+1)*sizeof(double)); 

//Инициализация Y0 and Y1

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

{

Y0[i]=0;

Y1[i]=0;

}

//printf("L=%d\n",L);

Continue:

if (rank < start)

{ 

Y0[L]=L;

d=rank*L;

}

else

{ 

Y0[L]=L-1;

d=start*L+(rank-start)*(L-1);

}

iter=iter+1;

flag=0;

if (rank < start) dim=L;

else dim=L-1;

// Умножение матрицы на вектор

for(k=0;k<size;k++) // Цикл по компьтерам

{

if (rank < start) dim=L;

else dim=L-1;

for (j=0;j<dim;j++) // Цикл по строкам

for (m=0,i=d;i<d+Y0[L];i++,m++) // Цикл по столбцам

Y1[j]+=A[j][i]*Y0[m];

d=(d+(int)Y0[L])%n;     

MPI_Sendrecv_replace(Y0,L+1,MPI_DOUBLE,next,tag,prev,tag,ring,&st);

}

// Вычисление решения на текущей итерации

if (rank < start) dim=L;

else dim=L-1;

for (j=0;j<dim;j++) // Цикл по строкaм

Y1[j]=Y0[j]-tau*(Y1[j]-A[j][n]);

//  for(j=0;j<dim;j++)

//  printf("Rank=%d Y1[%d]=%lf iter=%d\n",rank,j,Y1[j],iter);

//Проверка условия выхода

if (rank < start) dim=L;

else dim=L-1;

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

if (fabs((Y1[i]-Y0[i])/Y1[i])>Eps) flag=1;

//Посылка нулевому компьютеру данных, чтобы он принял решение

//об окончании или продолжении счета.

// if (flag==1)

if (rank!=0)

{

MPI_Send(&flag,1,MPI_INT,0,tag,ring); //send data to comp. number zero

MPI_Recv(&flag,1,MPI_INT,0,tag,ring,&st); //receive decision

}

else

{

//Receive datas from all computers

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

{

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

if(recv_flag==1) flag=1;

}

//Send to all computers the decision

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

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

}

if (flag==1)

{

if (rank < start) dim=L;

else dim=L-1;

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

{

Y0[i]=Y1[i];

Y1[i]=0;

}

goto Continue;

}

else

goto Finish;

Finish:

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

if(rank==0)

{

finishtime=time(NULL);

Result=(double*)malloc((size*(L+1))*sizeof(double));

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

Result[i]=Y1[i];

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

{

if (i < start)

MPI_Recv(&Result[i*L],L,MPI_DOUBLE,i,tag,ring,&st);

else

MPI_Recv(&Result[start*L+(i-start)*(L-1)],L-1,MPI_DOUBLE,i,tag,ring,&st);   

}

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

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

fprintf(fp,"%lf\n",Result[j]);

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

fprintf(fp,"\nЧисло итераций=%d",iter);

fclose(fp);

}

else

if (rank < start)

MPI_Send(Y1,L,MPI_DOUBLE,0,tag,ring);

else

MPI_Send(Y1,L-1,MPI_DOUBLE,0,tag,ring);

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

free(A[i]);

free(A);

free(Y0);

free(Y1);

MPI_Finalize();

return 0;

}