Параллельные вычислительные методы, страница 2

while(konelem[j][i]!=JG[d+l]) l++;

d+=l;

GG[d]+=loc[k][j];

}

}

}             

}             

}

// create portret global matrix

void portret() {

int k,j,p=0;

IG[0]=0;

for(k=displs[rank];k<displs[rank]+dim[rank];k++) {

IG[k+1-displs[rank]]=IG[k-displs[rank]];

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

if (poisk(k,j)) {

IG[k+1-displs[rank]]++;

JG[p]=j;

p++;

}

}

}

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

DI[j]=0.0;

F[j]=0.0;

}

}

// finding non null elements

int poisk(int fk,int fj) {

int l,q,f1,f2,flag=0;

for(l=0;l<kelem && flag==0;l++) {

f1=0;f2=0;

for(q=0;q<4;q++) {

if (konelem[q][l]==fk) f1=1;

if (konelem[q][l]==fj) f2=1;

}

if (f1+f2==2) flag=1;

}

return flag;

}

// solving slay

void msg() {

int i,j,k,ind;

double r[N],x[N],z[N],temp[N],a,b,norm_r=1.0,p,v,sum,ZZ[N],PP,VV,norm_RR=1.0;

FILE *out;

for(i=0;i<dim[rank];i++) x[i]=0;

for(i=0;i<dim[rank];i++) {

r[i]=F[i];

z[i]=r[i]/DI[i];

}

for(i=0;i<MAX && norm_RR>eps;i++) {

p=v=0;

for(k=0;k<dim[rank];k++)

p+=r[k]*r[k]/DI[k];

ind=0;

MPI_Allgatherv(z,dim[rank],MPI_DOUBLE,ZZ,dim,displs,MPI_DOUBLE,MPI_COMM_WORLD);

for(k=0;k<dim[rank];k++) {

sum=0;

for(j=0;j<IG[k+1]-IG[k];j++) {

sum+=GG[ind]*ZZ[JG[ind]];

ind++;

}

temp[k]=sum;

v+=temp[k]*z[k];    

}

PP=VV=0;

MPI_Allreduce(&p,&PP,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

MPI_Allreduce(&v,&VV,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

a=PP/VV;

v=0;

for(k=0;k<dim[rank];k++) {

x[k]+=a*z[k];

r[k]-=a*temp[k];

v+=r[k]*r[k]/DI[k];

}

VV=0.0;

MPI_Allreduce(&v,&VV,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

b=VV/PP;

norm_r=0;

for(k=0;k<dim[rank];k++) {

z[k]=r[k]/DI[k]+b*z[k];

norm_r+=r[k]*r[k];

}

norm_RR=0.0;

MPI_Allreduce(&norm_r,&norm_RR,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

}

out=fopen("f.txt","w");

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

fprintf(out,"%lf %lf %lf\n",uzel[0][k],uzel[1][k],x[k]);

fclose(out);

}

// kraevie usloviya

void ku() {

int i,ind;

for(i=0;i<dim[rank];i++) {

ind=IG[i];

if (uzel[2][i+displs[rank]]>eps) {

F[i]=CH*u(uzel[0][i+displs[rank]],uzel[1][i+displs[rank]]);

DI[i]=CH;

while(JG[ind]!=i+displs[rank]) ind++;

GG[ind]=CH;

}

}

}

// proverka na prinadlegnost inorodnomu objectu

int prover(int c1,int c2,int c3,int c4) {

int f=0;

if (uzel[1][c1]>=0.2 && uzel[1][c1]<=0.4 && uzel[0][c1]>=0.0 && uzel[0][c1]<=0.5)

f+=1;

if (uzel[1][c2]>=0.2 && uzel[1][c2]<=0.4 && uzel[0][c2]>=0.0 && uzel[0][c2]<=0.5)

f+=1;

if (uzel[1][c3]>=0.2 && uzel[1][c3]<=0.4 && uzel[0][c3]>=0.0 && uzel[0][c3]<=0.5)

f+=1;

if (uzel[1][c4]>=0.2 && uzel[1][c4]<=0.4 && uzel[0][c4]>=0.0 && uzel[0][c4]<=0.5)

f+=1;

if (f==4) return 1;

else return 0;

}

// generation setki

void gen() {

int x_kol,y_kol,i,j;

double x_left,x_right,y_left,y_right,x_shag,y_shag,x,y;

FILE *f1=fopen("uzel4.txt","w");

FILE *f2=fopen("konelem4.txt","w");

x_left=0.0;

x_right=1.0;

y_left=0.0;

y_right=1.0;

x_kol=11;

y_kol=11;

x_shag=(x_right-x_left)/(x_kol-1);

y_shag=(y_right-y_left)/(y_kol-1);

// ãåíåðàöèÿ ñåòêè - óçëîâ

fprintf(f1,"%d\n",x_kol*y_kol);

for(i=0;i<y_kol;i++) {

for(j=0;j<x_kol;j++) {

x=x_left+x_shag*j;

y=y_left+y_shag*i;

if (x==x_left || x==x_right) //|| y==y_left || y==y_right )

fprintf(f1,"%lf %lf %lf\n",x,y,1.0);

else

fprintf(f1,"%lf %lf %lf\n",x,y,0.0);

}

fprintf(f1,"\n");

}

fprintf(f2,"%d\n",(x_kol-1)*(y_kol-1));

for(i=0;i<y_kol-1;i++) {

for(j=0;j<x_kol-1;j++) {

fprintf(f2,"%d %d %d %d\n",i*x_kol+j,(i)*x_kol+(j+1),(i+1)*x_kol+j+1,(i+1)*x_kol+j);

}

fprintf(f2,"\n");

}

fclose(f1);

fclose(f2);

}

int main(int argc, char *argv[]) {

int i;

MPI_Init(&argc,&argv);

MPI_Comm_size(MPI_COMM_WORLD, &size);

MPI_Comm_rank(MPI_COMM_WORLD, &rank);

printf ("Rank=%d, size=%d.\n",rank,size);

gen();      rar();

ost=n%size;

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

dim[i]=n/size;

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

dim[i]++;

displs[0]=0;

for(i=1;i<size;i++) displs[i]=displs[i-1]+dim[i-1];

portret(); glob();     ku();        msg();      MPI_Finalize();     return(0); }