case 6:
return - 2*h* (mu*Kx+Ky) *simpsonFx (0.0,a,8, i,j) *simpsonFy (0.0,a,8, i,j);
case 7:
return - 2*h* (Kx+mu*Ky) *simpsonFx (0.0,a,9, i,j) *simpsonFy (0.0,a,9, i,j);
case 8:
return - 2*h* (mu*Kx+Ky) *simpsonFx (0.0,a,8, i,j) *simpsonFy (0.0,a,8, i,j);
case 9:
return 2*h* (Kx*Kx+2*mu*Kx*Ky+Ky*Ky) *simpsonFx (0.0,a,10, i,j) *simpsonFy (0.0,a,10, i,j) +pow (h,3) /12
* (2*simpsonFx (0.0,a,11, i,j) *simpsonFy (0.0,a,10, i,j) +2*mu*simpsonFx (0.0,a,12, i,j)
*simpsonFy (0.0,a,12, i,j) +2*mu*simpsonFx (0.0,a,12, i,j) *simpsonFy (0.0,a,12, i,j) +2*simpsonFx (0.0,a,10, i,j)
*simpsonFy (0.0,a,11, i,j) +8*mu1*simpsonFx (0.0,a,13, i,j) *simpsonFy (0.0,a,13, i,j));
case 10:
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) {
if ( (row >= maxsize) || (column >= maxsize)
|| (row<0) || (column<0)) return false;
data [row * maxsize + column] = newvalue;
return true;
};
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] ;