return 2* (1-pow (mu,2)) *simpsonFx (0.0,a,14, i,j) *simpsonFy (0.0,a,14, i,j);

default:

return 0;

}

}

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

{

// объявление переменных.

int myid, numprocs;

int namelen;

char processor_name [MPI_MAX_PROCESSOR_NAME] ;

double startwtime = 0.0, endwtime;

int rc;

MPI_Status status;

rc = MPI_Init (&argc,&argv);

rc|= MPI_Comm_size (MPI_COMM_WORLD,&numprocs);

rc|= MPI_Comm_rank (MPI_COMM_WORLD,&myid);

if (rc! = 0)

printf ("error initializing MPI and obtaining task ID information\n");

MPI_Get_processor_name (processor_name,&namelen);

fprintf (stdout,"Process%d of%d is on%s\n",

myid, numprocs, processor_name);

fflush (stdout);

// функция начала замера времени вычисления.

if (myid == 0) {

startwtime = MPI_Wtime ();

}

N = 4; // количество членов

a = 1, b = 1;

n = 6, m = 6, E = 21000;

q = 0.0037, h = 0.001, mu = 0.3, r = 0.2;

Pi = 3.14;

int i, j;

double rv;

bool xyz;

R1 = 440*h;

R2 = 440*h;

mu1 = (1-mu) /2;

Kx = 1/ (double) R1;

Ky = 1/ (double) R2;

// выделение памяти под массивы для аппроксимирующих функций

X1 = (double*) malloc (N * sizeof (double));

X2 = (double*) malloc (N * sizeof (double));

X3 = (double*) malloc (N * sizeof (double));

Y1 = (double*) malloc (N * sizeof (double));

Y2 = (double*) malloc (N * sizeof (double));

Y3 = (double*) malloc (N * sizeof (double));

int sqrtN = pow (N, 0.5);

/*

вычисление коэффициентов аппроксимирующих функций на

нескольких процессах.

*/

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

{

for (j = 1; j <= sqrtN; j++)

{

if (myid == 0) {

X1 [sqrtN * (i - 1) + j - 1] = 2 * i * Pi;

printf ("X1 [%d] =%.3f\n",sqrtN * (i - 1) + j - 1,X1 [sqrtN * (i - 1) + j - 1]);

X2 [sqrtN * (i - 1) + j - 1] = (2 * i - 1) * Pi;

printf ("X2 [%d] =%.3f\n",sqrtN * (i - 1) + j - 1,X2 [sqrtN * (i - 1) + j - 1]);

X3 [sqrtN * (i - 1) + j - 1] = (2 * i - 1) * Pi;

printf ("X3 [%d] =%.3f\n",sqrtN * (i - 1) + j - 1,X3 [sqrtN * (i - 1) + j - 1]);

}

if (myid == 1) {

Y1 [sqrtN * (i - 1) + j - 1] = (2 * j - 1) * Pi;

printf ("Y1 [%d] =%.3f\n",sqrtN * (i - 1) + j - 1,Y1 [sqrtN * (i - 1) + j - 1]);

Y2 [sqrtN * (i - 1) + j - 1] = 2 * j * Pi;

printf ("Y2 [%d] =%.3f\n",sqrtN * (i - 1) + j - 1,Y2 [sqrtN * (i - 1) + j - 1]);

Y3 [sqrtN * (i - 1) + j - 1] = (2 * j - 1) * Pi;

printf ("Y3 [%d] =%.3f\n",sqrtN * (i - 1) + j - 1,Y3 [sqrtN * (i - 1) + j - 1]);

}

}

}

/*

пересылка результатов вычислений на "головную" машину

*/

if (myid == 1) {

MPI_Send (Y1, N, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);

MPI_Send (Y2, N, MPI_DOUBLE, 0, 2, MPI_COMM_WORLD);

MPI_Send (Y3, N, MPI_DOUBLE, 0, 3, MPI_COMM_WORLD);

MPI_Recv (X1, N, MPI_DOUBLE, 0, 4, MPI_COMM_WORLD,

&status);

MPI_Recv (X2, N, MPI_DOUBLE, 0, 5, MPI_COMM_WORLD,

&status);

MPI_Recv (X3, N, MPI_DOUBLE, 0, 6, MPI_COMM_WORLD,

&status);

}

if (myid == 0) {

MPI_Recv (Y1, N, MPI_DOUBLE, 1, 1, MPI_COMM_WORLD,

&status);

MPI_Recv (Y2, N, MPI_DOUBLE, 1, 2, MPI_COMM_WORLD,

&status);

MPI_Recv (Y3, N, MPI_DOUBLE, 1, 3, MPI_COMM_WORLD,

&status);

MPI_Send (X1, N, MPI_DOUBLE, 1, 4, MPI_COMM_WORLD);

MPI_Send (X2, N, MPI_DOUBLE, 1, 5, MPI_COMM_WORLD);

MPI_Send (X3, N, MPI_DOUBLE, 1, 6, MPI_COMM_WORLD);

}

// вывод времени вычисления аппрокс. функций

if (myid == 0) {

endwtime = MPI_Wtime ();

printf ("\napp func clock time =%f\n", endwtime-startwtime);

}

printf ("\n - --------------- - BEGIN - -----------------\n");

/*

выделение памяти под массивы коэффициентов ФПЭД

и вычисление их на разных процессах.

*/

C1 = (double*) malloc (3 * N * 3 * N * sizeof (double));

C2 = (double*) malloc (3 * N * 3 * N * sizeof (double));

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

{

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

{

// обнуление всех значений для удобства пересылки

C1 [i*N+j] =0;

C1 [i*N+N+j] =0;

C1 [i*N+2*N+j] =0;

C1 [N+i*N+j] =0;

C1 [N+i*N+N+j] =0;

C1 [N+i*N+2*N+j] =0;

C1 [2*N+i*N+j] =0;

C1 [2*N+i*N+N+j] =0;

C1 [2*N+i*N+2*N+j] =0;

if (myid == 0) {

C1 [N*i+j] =C (1, i,j);

C1 [N*i+N+j] =C (2, i,j);

C1 [N*i+2*N+j] =C (3, i,j);

C1 [N+N*i+j] =C (4, i,j);

C1 [N+N*i+N+j] =C (5, i,j);

}

if (myid == 1) {

C1 [N+N*i+2*N+j] =C (6, i,j);

C1 [2*N+N*i+j] =C (7, i,j);

C1 [2*N+N*i+N+j] =C (8, i,j);

C1 [2*N+N*i+2*N+j] =C (9, i,j);

}

}

}

// пересылка массивов на "головную" машину

if (myid == 1) {

MPI_Send (C1, 3*N*3*N, MPI_DOUBLE, 0, 7, MPI_COMM_WORLD);

}

if (myid == 0) {

MPI_Recv (C2, 3*N*3*N, MPI_DOUBLE, 1, 7, MPI_COMM_WORLD, &status);

printf ("\n\nC2 [1]%.3f\n",C2 [0]);

}

printf ("\n - --------------- - END - -----------------\n");

if (myid == 0) {

matrix <double> M1 (3*N,3*N);

printf ("-------------------- - BEGIN FIRST - -----------------\n");

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

{

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

{

M1. setvalue (i,j,C (1, i,j));

printf ("C1 [%d,%d]: =%.5f ", i,j,C (1, i,j));

M1. setvalue (i,N+j,C (2, i,j));

printf ("C2 [%d,%d]: =%.5f ", i,N+j,C (2, i,j));

M1. setvalue (i,2*N+j,C (3, i,j));

printf ("C3 [%d,%d]: =%.5f\n", i,2*N+j,C (3, i,j));

M1. setvalue (N+i,j,C (4, i,j));

printf ("C4 [%d,%d]: =%.5f ",N+i,j,C (4, i,j));

M1. setvalue (N+i,N+j,C (5, i,j));

printf ("C5 [%d,%d]: =%.5f ",N+i,N+j,C (5, i,j));

M1. setvalue (N+i,2*N+j,C (6, i,j));

printf ("C6 [%d,%d]: =%.5f\n",N+i,2*N+j,C (6, i,j));

M1. setvalue (2*N+i,j,C (7, i,j));

printf ("C7 [%d,%d]: =%.5f ",2*N+i,j,C (7, i,j));

M1. setvalue (2*N+i,N+j,C (8, i,j));

printf ("C8 [%d,%d]: =%.5f ",2*N+i,N+j,C (8, i,j));

M1. setvalue (2*N+i,2*N+j,C (9, i,j));

printf ("C9 [%d,%d]: =%.5f\n",2*N+i,2*N+j,C (9, i,j));

}

}

printf ("-------------------- - END FIRST - -----------------\n");

printf ("-------------------- - BEGIN SECOND - -----------------\n");

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

{

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

{

M1. setvalue (i,j,C1 [N*i+j]);

printf ("C1 [%d,%d]: =%.5f ", i,j,C1 [N*i+j]);

M1. setvalue (i,N+j,C1 [N*i+N+j]);

printf ("C2 [%d,%d]: =%.5f ", i,N+j,C1 [N*i+N+j]);

M1. setvalue (i,2*N+j,C1 [N*i+2*N+j]);

printf ("C3 [%d,%d]: =%.5f\n", i,2*N+j,C1 [N*i+2*N+j]);

M1. setvalue (N+i,j,C1 [N+i+j]);

printf ("C4 [%d,%d]: =%.5f ",N+i,j,C1 [N+N*i+j]);

M1. setvalue (N+i,N+j,C1 [N+N*i+N+j]);

printf ("C5 [%d,%d]: =%.5f ",N+i,N+j,C1 [N+N*i+N+j]);

M1. setvalue (N+i,2*N+j,C2 [N+N*i+2*N+j]);

printf ("C6 [%d,%d]: =%.5f\n",N+i,2*N+j,C2 [N+N*i+2*N+j]);

M1. setvalue (2*N+i,j,C2 [2*N+N*i+N+j]);

printf ("C7 [%d,%d]: =%.5f ",2*N+i,j,C2 [2*N+N*i+N+j]);

M1. setvalue (2*N+i,N+j,C2 [2*N+N*i+N+j]);

printf ("C8 [%d,%d]: =%.5f ",2*N+i,N+j,C2 [2*N+N*i+N+j]);

M1. setvalue (2*N+i,2*N+j,C2 [2*N+N*i+2*N+j]);

printf ("C9 [%d,%d]: =%.5f\n",2*N+i,2*N+j,C2 [2*N+N*i+2*N+j]);

}

}

printf ("-------------------- - END SECOND - -----------------\n");

// заполнение массивов свободных членов

matrix <double> Pow (3*N,3*N);

for (i=0; i < Pow. getactualsize (); i++)

{

for (j=0; j< Pow. getactualsize (); j++)

{

Pow. setvalue (i,j,0);

}

}

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

{

Pow. setvalue (i,0,0);

Pow. setvalue (N+i,0,0);

Pow. setvalue (2*N+i,0,q/E*C (10, i,0));

}

printf ("\n\n\nM1: ");

for (i=0; i < M1. getactualsize (); i++)

{

printf ("\n%d: ", i);

for (j=0; j< M1. getactualsize (); j++)

{

M1. getvalue (i,j,rv,xyz);

printf ("%.10f ",rv);

}

}

// выделение памяти под матрицы

matrix <double> M2 (3*N,3*N);

matrix <double> M3 (3*N,3*N);

// копирование матрицы коэффициентов

M2. copymatrix (M1);

// обращение матрицы

M1. invert ();

// произведение матриц

M3. settoproduct (M1,M2);

// сравнение полученной единичной матрицы с эталоном единичной матрицы

M3.comparetoidentity ();

printf ("\n\n\ninverse: ");

for (i=0; i < M3. getactualsize (); i++)

{

printf ("\n%d: ", i);

for (j=0; j< M3. getactualsize (); j++)

{

M3. getvalue (i,j,rv,xyz);

printf ("%.10f ",rv);

}

}

for (i=0; i < M1. getactualsize (); i++)

{

printf ("\n%d: ", i);

for (j=0; j< M1. getactualsize (); j++)

{

M1. getvalue (i,j,rv,xyz);

printf ("%.10f ",rv);

}

}

// выделение памяти для матрицы результатов

matrix <double> Ans (3*N,3*N);

Ans. settoproduct (M1,Pow);

printf ("\nanswer \n");

for (i=0; i < Ans. getactualsize (); i++)

{

printf ("%d: ", i);

Ans. getvalue (i,0,rv,xyz);

printf ("%.10f ", rv);

printf ("\n");

}

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

W = 0;

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

{

Ans. getvalue (2*N+i,0,rv,xyz);

printf ("!!%.10f%.10f%.10f\n", rv, X3 [i], Y3 [i]);

W += rv * sin_0 (X3 [i],a/2) * sin_0 (Y3 [i],a/2);

}

printf ("W: =%.10f", W);

// вывод времени вычисления программы

endwtime = MPI_Wtime ();

printf ("\nwall clock time =%f\n", endwtime-startwtime);

fflush (stdout);

}

MPI_Finalize ();

return 0;

}

Matrix. h:

#ifndef __mjdmatrix_h

#define __mjdmatrix_h

template <class D> class matrix{

int maxsize;

int actualsize;

D* data;

void allocate () {

delete [] data;

data = new D [maxsize*maxsize] ;

};

matrix () {};

matrix (int newmaxsize) {matrix (newmaxsize,newmaxsize); };

public:

matrix (int newmaxsize, int newactualsize) {

if (newmaxsize <= 0) newmaxsize = 5;

maxsize = newmaxsize;

if ( (newactualsize <= newmaxsize) && (newactualsize>0))

actualsize = newactualsize;

else

actualsize = newmaxsize;

data = 0;

allocate ();

};

~matrix () { delete [] data; };

void settoproduct (matrix& left, matrix& right) {

actualsize = left. getactualsize ();

if (maxsize < left. getactualsize ()) {

maxsize = left. getactualsize ();

allocate ();

}

for (int i = 0; i < actualsize; i++)

for (int j = 0; j < actualsize; j++) {

D sum = 0.0;

D leftvalue, rightvalue;

bool success;

for (int c = 0; c < actualsize; c++) {

left. getvalue (i,c,leftvalue,success);

right. getvalue (c,j,rightvalue,success);

sum += leftvalue * rightvalue;

}

setvalue (i,j,sum);

}

}

void combine (matrix& left, matrix& right) {

actualsize = left. getactualsize ();

if (maxsize < left. getactualsize ()) {

maxsize = left. getactualsize ();

allocate ();

}

for (int i = 0; i < actualsize; i++)

for (int j = 0; j < actualsize; j++) {

D sum = 0.0;

D leftvalue, rightvalue;

bool success;

left. getvalue (i,j,leftvalue,success);

right. getvalue (i,j,rightvalue,success);

sum = leftvalue + rightvalue;

setvalue (i,j,sum);

}

}

void setactualsize (int newactualsize) {

if (newactualsize > maxsize)

{

maxsize = newactualsize;

allocate ();

}

if (newactualsize >= 0) actualsize = newactualsize;

};

int getactualsize () { return actualsize; };

void getvalue (int row, int column, D& returnvalue, bool& success) {

if ( (row>=maxsize) || (column>=maxsize)

|| (row<0) || (column<0))

{ success = false;

return; }

returnvalue = data [row * maxsize + column] ;

success = true;

};

bool setvalue (int row, int column, D newvalue) (row<0) ;

void dumpMatrixValues () {

bool xyz;

double rv;

for (int i=0; i < actualsize; i++)

{

std:: cout << "i=" << i << ": ";

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

{

M. getvalue (i,j,rv,xyz);

std:: cout << rv << " ";

}

std:: cout << std:: endl;

}

};

void comparetoidentity () {

int worstdiagonal = 0;

D maxunitydeviation = 0.0;

D currentunitydeviation;

for (int i = 0; i < actualsize; i++) {

currentunitydeviation = data [i*maxsize+i] - 1.;

if (currentunitydeviation < 0.0) currentunitydeviation *= - 1.;

if (currentunitydeviation > maxunitydeviation) {

maxunitydeviation = currentunitydeviation;

worstdiagonal = i;

}

}

int worstoffdiagonalrow = 0;

int worstoffdiagonalcolumn = 0;

D maxzerodeviation = 0.0;

D currentzerodeviation;

for (int i = 0; i < actualsize; i++) {

for (int j = 0; j < actualsize; j++) {

if (i == j) continue;

currentzerodeviation = data [i*maxsize+j] ;

if (currentzerodeviation < 0.0) currentzerodeviation *= - 1.0;

if (currentzerodeviation > maxzerodeviation) {

maxzerodeviation = currentzerodeviation;

worstoffdiagonalrow = i;

worstoffdiagonalcolumn = j;

}

}

}

printf ("Worst diagonal value deviation from unity:%0.5f at row/column%0.3f\n", maxunitydeviation, worstdiagonal);

printf ("Worst off-diagonal value deviation from zero:%0.5f at row%0.3f, column%0.3f\n", maxzerodeviation, worstoffdiagonalrow,worstoffdiagonalcolumn);

};

void copymatrix (matrix& source) {

actualsize = source. getactualsize ();

if (maxsize < source. getactualsize ()) {

maxsize = source. getactualsize ();

allocate ();

}

for (int i = 0; i < actualsize; i++)

for (int j = 0; j < actualsize; j++) {

D value;

bool success;

source. getvalue (i,j,value,success);

data [i*maxsize+j] = value;

}

};

void invert () {

if (actualsize <= 0) return;

if (actualsize == 1) return;

for (int i=1; i < actualsize; i++) data [i] /= data [0] ;

for (int i=1; i < actualsize; i++) {

for (int j=i; j < actualsize; j++) {

D sum = 0.0;

for (int k = 0; k < i; k++)

sum += data [j*maxsize+k] * data [k*maxsize+i] ;

data [j*maxsize+i] - = sum;

}

if (i == actualsize-1) continue;

for (int j=i+1; j < actualsize; j++) {

D sum = 0.0;

for (int k = 0; k < i; k++)

sum += data [i*maxsize+k] *data [k*maxsize+j] ;

data [i*maxsize+j] =

(data [i*maxsize+j] -sum) / data [i*maxsize+i] ;

}

}

for (int i = 0; i < actualsize; i++) // invert L

for (int j = i; j < actualsize; j++) {

D x = 1.0;

if (i! = j) {

x = 0.0;

for (int k = i; k < j; k++)

x - = data [j*maxsize+k] *data [k*maxsize+i] ;

}

data [j*maxsize+i] = x / data [j*maxsize+j] ;

}

for (int i = 0; i < actualsize; i++)

for (int j = i; j < actualsize; j++) {

if (i == j) continue;

D sum = 0.0;

for (int k = i; k < j; k++)

sum += data [k*maxsize+j] * ( (i==k)? 1.0: data [i*maxsize+k]);

data [i*maxsize+j] = - sum;

}

for (int i = 0; i < actualsize; i++)

for (int j = 0; j < actualsize; j++) {

D sum = 0.0;

for (int k = ( (i>j)? i: j); k < actualsize; k++)

sum += ( (j==k)? 1.0: data [j*maxsize+k]) *data [k*maxsize+i] ;

data [j*maxsize+i] = sum;

}

};

};

#endif

Страницы: 1, 2, 3, 4, 5, 6