Программа решения СЛАУ методом Рунге-Кутта 4го порядка.

Министерство образования и науки РФ

Новосибирский государственный технический университет

Лабораторная работа №1

дисциплина:

Параллельные вычислительные технологии

Факультет: ПМИ

Группа: ПМ-23

Студенты: Воронов А. С.

                    Липатников Р. Е.

Преподаватель:

                    Куликов

Новосибирск

2006

Постановка задачи

Решить СЛАУ методом Рунге-Кутта 4го порядка.

СЛАУ вида

Текст программы

#include <stdio.h>

#include <math.h>

#include <mpi.h>

const int N=57;

const int S=4;

double tau=0.001;

int size,rank,blocksize,c,n,n1;

MPI_Status st;

int GenMatrix (double A[][N])

{

int i,j,begin;

if (rank<c) begin=rank*n;

else begin=c*n+(rank-c)*n1;

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

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

{

if (begin+i==j) A[i][j]=1;

else A[i][j]=1.0/N;}

return 1;  

}

int GenVect(double *y, int n)

{ int i;

for (i=0;i<n;i++) y[i]=1;

return 1;

}

int PrintVector (double *Vect, int n)

{ int i;

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

printf ("Vect%d[%d] = %e;\n",rank,i,Vect[i]);

return 1;

}

int Multiply (double A[][N], double *Vect, double *Res, int n_b)

{ int i,j,k,m,dn;

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

Res[i]=0;

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

{                  if ((k+rank)%size<=c)

dn=(k+rank)%size*n;

else

dn=(c*n+n1*(k+rank-c))%N;

if ((k+rank)%size<c)

m=n;

else

m=n1;

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

{

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

{             Res[i]+=A[i][dn+j]*Vect[j]; }

}

MPI_Send(Vect,m,MPI_INT,(rank-1+size)%size,16,MPI_COMM_WORLD);

if ((k+(rank+1)%size)%size<c)

m=n;

else

m=n1;

MPI_Recv(Vect,m,MPI_INT,(rank+1)%size,16,MPI_COMM_WORLD,&st);

}                 

return 1;

}

int MultiplyScal (double *Vect, double *Res, double scal,int n)

{

int i;

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

{  Res[i]=scal*Vect[i]; }

return 1;

}

int Sum (double *Vect1, double *Vect2, double scal,double *Res, int n)

{ int i;

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

Res[i]=Vect1[i]+scal*Vect2[i];

return 1;

}

int main(int argc,char **argv)

{ double A[N/S+1][N],k1[N/S+1],k2[N/S+1],k3[N/S+1],k4[N/S+1],y[N/S+1],tmp[N],tmp2[N];

double t;

int i,k; 

MPI_Init(&argc,&argv);

MPI_Comm_rank (MPI_COMM_WORLD, &rank);

MPI_Comm_size (MPI_COMM_WORLD, &size);

c=N%size;

n1=N/size;

n=n1+1;

if (rank<c) blocksize=n;

else blocksize=n1;

GenMatrix(A);

GenVect(y,N);

k=0;

for (t = tau; t < 1; t =t+tau)

{

if (rank==0) printf ("Iter = %d\n",k);

k++;

// k1

Multiply (A,y,tmp,N);

MultiplyScal (tmp,k1,tau,blocksize);

// k2

Sum (y,k1,0.5,tmp,blocksize);

//PrintVector (tmp,2);

Multiply (A,tmp,tmp2,N);

// PrintVector (tmp2,blocksize);

MultiplyScal (tmp2,k2,tau,blocksize);

// k3

Sum (y,k2,0.5,tmp,blocksize);

Multiply (A,tmp,tmp2,N);

MultiplyScal (tmp2,k3,tau,blocksize);

// k4

Sum (y,k3,1,tmp,blocksize);

Multiply (A,tmp,tmp2,N);

MultiplyScal (tmp2,k4,tau,blocksize);

//y+ (k1+2*k2+2*k3+k4) /6

Sum (k1,k2,2,tmp,blocksize);

Sum (tmp,k3,2,tmp,blocksize);

Sum (tmp,k4,1,tmp,blocksize);

Sum (y,tmp,1.0/6.0,y,blocksize); //*/

}

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

printf ("Res%d[%d] = %8.5e, Otvet = %8.5e, Eps = %e;\n",rank,i,y[i], exp(t*(2.0*N-1.0)/N),y[i]-exp(t*(2.0*N-1.0)/N));

MPI_Finalize();

return 1;

}