printf ("\nValue = %G (total point = %G)\n", _s, totpt);
// Интеграл и кол-во точек
printf ("Current error = %G\n", fabsl(s - _s));
// Абсолютная ошибка
totpt *= mult;
} while (fabsl(s - _s) > eps && totpt <= TOTPT_MAX);
stop = clock(); // Останавливаем таймер
TotTimePrint(start, stop, NULL, 'd');
// Печатаем время выполнения цикла
return _s;
}
///////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral1DMK_par (CURVLINE f, LIMIT left,
LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult)
/* Вычисление однократного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции */
// П Р О Т Е С Т И Р О В А Н А
// НА КЛАСТЕРЕ ! ! ! – 26.05.04 г.
{
// Проверка корректности введенных значений
if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1);
int totproc, self, root = 0;
MPI_Comm comw = MPI_COMM_WORLD;
MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов
MPI_Comm_rank(comw, &self); // Получение процессом своего номера
FUNC_TYPE _s = 0.0, _s_root = 0.0, s;
ARG_TYPE h;
//srand(self);
//MPI_Barrier(comw);
//double mpi_start = MPI_Wtime();
do {
s = _s;
_s = 0.0;
// Инициализация генератора случайных чисел номером процесса self
srand((UINT)(time(NULL) >> self));
do { // Обеспечиваем малый шаг
h = (right - left) / totpt;
if (h > 1) totpt *= mult;
} while (h > 1);
if (!h) exit(1);
/* Распараллеливание цикла накапливания суммы _s значений подынтегральной функции f в случайных точках, лежащих в пределах [left, left + i*h] */
ULONG i = self + 1;
while ((h > 0 && left + i*h <= right) ||
(h < 0 && left + i*h >= right)) {
// Накапливаем сумму
_s += f(random(left + (i-1)*h, left + i*h));
i += totproc;
}
// Синхронизация всех процессов, т.к процессы выполняют
// разное кол-во итераций в цикле
MPI_Barrier(comw);
// Суммирование _s во всех процессах и сохранение полученного
// значения в корневом процессе в переменной _s_root
MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw);
if (!self) { // В корневом процессе:
// Рассчет интеграла по формуле (3.1.10)
_s = (_s_root/totpt)*(right-left);
// Вывод на экран корневого процесса промежуточных данных
printf("\nValue = %G (total point = %G)\n",
_s, totpt);
printf ("Current error = %G\n", fabsl(s - _s));
}
// Синхронизация всех процессов
MPI_Barrier(comw);
/* Широковещание полученного значения интеграла из корнево го процесса всем остальным, включая самого себя */
MPI_Bcast(&_s, 1, MPI_LONG_DOUBLE, root, comw);
/* Если ошибка вычислений больше допустимой, увеличение
числа разыгрываемых точек в mult раз */
totpt *= mult;
} while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX);
//double mpi_stop = MPI_Wtime();
//printf("Total time = %f\n", mpi_stop-mpi_start);
return _s;
}
///////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral1DMK_par1 (CURVLINE f, LIMIT left,
LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult)
/* Вычисление однократного интеграла методом Монте-Карло на кластере MPI Вариант со средним значением функции. Каждый процесс считает интеграл на уменьшенном интервале*/
// П Р О Т Е С Т И Р О В А Н А
// НА КЛАСТЕРЕ ! ! ! – 26.05.04 г.
{
// Проверка корректности введенных значений
if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1);
int totproc, self, root = 0;
MPI_Comm comw = MPI_COMM_WORLD;
MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов
MPI_Comm_rank(comw, &self); // Получение процессом своего номера
FUNC_TYPE _s = 0.0, _s_root = 0.0, s;
ARG_TYPE h, H, left_lim;
// Для каждого из процессов определяем свой интервал интегрирования
H = (ARG_TYPE)(right - left) / totproc;
left += self*H;
right = left + H;
printf("%f\t%f\n", left, right);
//srand(self);
//MPI_Barrier(comw);
//double mpi_start = MPI_Wtime();
do {
s = _s;
// Запоминаем предыдущее значение интеграла на промежутке [left, right]
_s = 0.0; // Обнуляем сумму
// Инициализируем генератор случайных чисел временем
srand((UINT)(time(NULL) >> self));
do { // Обеспечиваем малый шаг
h = (right - left) / totpt;
if (h > 1) totpt *= mult;
} while (h > 1);
if (!h) exit(1);
left_lim = left;
// Накапливаем сумму
while ((h > 0 && left_lim + h <= right) ||
(h < 0 && left_lim + h >= right)) {
_s += f(random(left_lim, left_lim+h));
left_lim += h;
}
// Вычисляем интеграл "через" среднее арифметическое функции
_s = (_s/totpt)*(right-left);
totpt *= mult;
MPI_Barrier(comw); // Cинхронизируем процессы
// Суммируем значения интеграла, полученные на малых промежутках
MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw);
if (!self) {
_s = _s_root;
// Вывод на экран промежуточных данных
printf("\nValue = %G (total point = %G)\n",
_s, totpt); // Интеграл и кол-во точек
printf ("Current error = %G\n", fabsl(s - _s));
// Абсолютная ошибка
}
MPI_Barrier(comw); // Cинхронизируем процессы
MPI_Bcast(&_s, 1, MPI_LONG_DOUBLE, root, comw);
} while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX);
//double mpi_stop = MPI_Wtime();
//printf("Total time = %f\n", mpi_stop-mpi_start);
return _s;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral1DMK_par2 (CURVLINE f, LIMIT left,
LIMIT right, ERR_T eps, PTCOUNT totpt, ULONG mult)
/* Вычисление однократного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции. Каждый процесс вычисляет интеграл на уменьшенном интервале */
// П Р О Т Е С Т И Р О В А Н А
// НА КЛАСТЕРЕ ! ! ! - 26.05.04 г.
{
// Проверка корректности введенных значений
if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1);
int totproc, self, root = 0;
MPI_Comm comw = MPI_COMM_WORLD;
MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов
MPI_Comm_rank(comw, &self); // Получение процессом своего номера
FUNC_TYPE _s = 0.0, _s_root = 0.0, s;
ARG_TYPE h, H, left_lim;
// Для каждого из процессов определяем свой интервал интегрирования
H = (ARG_TYPE)(right - left)/totproc;
left += self*H;
right = left + H;
//srand(self);
//MPI_Barrier(comw);
//double mpi_start = MPI_Wtime();
do {
s = _s;
// Запоминаем предыдущее значение интеграла на промежутке [left, right]
_s = 0.0; // Обнуляем сумму
// Инициализируем генератор случайных чисел временем
srand((UINT)(time(NULL) >> self));
do { // Обеспечиваем малый шаг
h = (right - left) / totpt;
if (h > 1) totpt *= mult;
} while (h > 1);
if (!h) exit(1);
left_lim = left;
// Накапливаем сумму
while ((h > 0 && left_lim + h <= right) ||
(h < 0 && left_lim + h >= right)) {
_s += f(random(left_lim, left_lim+h));
left_lim += h;
}
// Вычисляем интеграл по формуле (3.1.10)
_s = (_s/totpt)*(right-left);
totpt *= mult;
} while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX);
MPI_Barrier(comw);
MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw);
//double mpi_stop = MPI_Wtime();
//printf("Total time = %f\n", mpi_stop-mpi_start);
return _s_root;
}
/////////////////////////////////////////////////////////////////////////////////////////////
UFUNC_TYPE Square
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps)
// Вычисление площади под кривой f(x) по формуле Симпсона
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 24.05.04 г.
{
if (eps <= 0) exit (EPS_ERR);
ARG_TYPE _s, s = 0.0;
ARG_TYPE t0 = fabsl(f(left)) + fabsl(f(right));
ULONG n=1024;
do {
_s = s;
ARG_TYPE h = (right - left) / (2*n);
ARG_TYPE t1 = 0, t2 = 0;
// t1 – сумма нечетных, t2 – сумма четных
for (ULONG j = 1; j <= 2*n-1; j++)
j % 2 == 0 ? t2 += fabsl(f(left+j*h)) : t1 += fabsl(f(left+j*h));
s = (t0 + 4*t1 + 2*t2) * (fabsl(h)/3);
n *= 2; // mult == 2
} while (fabsl(s-_s) > eps && n <= ULONG_MAX);
return s;
}
///////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral2DMK
(SURF surf, LIMIT left, LIMIT right,
CURVLINE bottom, CURVLINE top, unsigned char var,
ERR_T eps, PTCOUNT totpt, ULONG mult)
// Рассчет двойного интеграла для у-правильной области
{
// Проверка корректности введенных значений
if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1);
FUNC_TYPE V = 0.0, _s = 0.0, s;
clock_t start = clock(), stop; // Запускаем таймер
do {
s = _s;
_s = 0.0;
// Инициализируем генератор случайных чисел временем
srand((UINT)(time(NULL)));
ARG_TYPE x, y;
for (ULONG i = 1; i <= totpt; i++) {
x = random(left, right);
y = random(bottom(x), top(x));
_s += surf(x, y);
}
_s = _s/totpt; // Среднее значение функции
// Вывод на экран промежуточных данных
printf("\nValue = %G (total point = %G)\n",
_s, totpt);
printf ("Current error = %G\n", fabsl(s - _s));
totpt *= mult;
} while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX);
// Вычисляем площадь области интегрирования
FUNC_TYPE field = RegionY(bottom, top, left, right, var, eps);
printf("region = %f\n", field);
// Вычисляем двойной интеграл
V = _s*field;
stop = clock();
TotTimePrint(start, stop, NULL, 'd');
return V;
}
///////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral2DMK_par
(SURF surf, LIMIT left, LIMIT right,
CURVLINE bottom, CURVLINE top, unsigned char var,
ERR_T eps, PTCOUNT totpt, ULONG mult)
// Рассчет двойного интеграла для у-правильной области
{
// Проверка корректности введенных значений
if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1);
FUNC_TYPE V = 0.0, _s = 0.0, _s_root = 0.0, s;
int totproc, self, root = 0;
MPI_Comm comw = MPI_COMM_WORLD;
MPI_Comm_size(comw, &totproc); // Получение кол-ва процессов
MPI_Comm_rank(comw, &self); // Получение процессом своего номера
do {
s = _s;
_s = 0.0;
// Инициализируем генератор случайных чисел временем
srand((UINT)(time(NULL) >> self));
ARG_TYPE x, y;
for (ULONG i = 1; i <= totpt; I += totproc) {
x = random(left, right);
y = random(bottom(x), top(x));
_s += surf(x, y);
}
MPI_Barrier(comw);
MPI_Reduce(&_s, &_s_root, 1, MPI_LONG_DOUBLE, MPI_SUM, root, comw);
if (!self) {
_s = _s_root/totpt; // Среднее значение функции
// Вывод на экран промежуточных данных
printf("\nValue = %G (total point = %G)\n", _s, totpt);
printf ("Current error = %G\n", fabsl(s - _s));
}
totpt *= mult;
MPI_Barrier(comw);
MPI_Bcast (&_s, 1, MPI_LONG_DOUBLE, root, comw);
} while (fabsl(s – _s) > eps && totpt <= TOTPT_MAX);
// Вычисляем площадь области интегрирования
FUNC_TYPE field = RegionY(bottom, top, left, right, var, eps);
printf("region = %f\n", field);
// Вычисляем двойной интеграл
V = _s*field; return V; }
[1] Операция понимается в самом широком смысле этого слова: операция интегрирования, операция решения СЛАУ и т.д.