Алгоритмы параллельных процессов при исследовании устойчивости подкрепленных пологих оболочек

Одно из главных правил оптимизации заключается в упрощении действий, выполняемых многократно в циклах. Логично было бы применить этот же принцип и для распараллеливания.Т. е. "рутинная" многократная работа будет разделена между определенным количеством процессоров. Так как нет смысла распараллеливать всю программу, можно выделить "головную" ЭВМ, которая будет выполнять основные шаги вычисления. Остальные же машины будут обрабатывать пакет задач, принимаемые от ведущей машины.Т. е. получаем схему коммуникации процессоров, похожую на схему "клиент-сервер", за исключением того, что клиенты также могут обмениваться сообщениями, но только при согласовании своих действий с сервером (рис.14).

87

Рис.14. Схема коммуникационной среды (ЦП - "головная" ЭВМ, П - "ведомая" ЭВМ)

Основные проблемы, которые могут возникнуть в процессе распределенных вычислений:

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

возникновение тупиковых ситуаций при приеме/передаче сообщений между процессорами;

неправильное разделение заданий между процессорами.

Ведущая ЭВМ на этапах распараллеливания будет освобождена от выполнения "рутинных" вычислений. Взамен она будет контролировать распределение заданий между ведомыми ЭВМ, фиксировать транзакции, отслеживать загруженность процессоров, принимать решение в случае возникновения ошибок, перенаправлять задания / повторно высылать задания в случае возникновения сбоя в среде коммуникации, собирать и проверять результаты вычислений.

Данная схема рассчитана на большое количество ЭВМ в коммуникационной среде.

На данном этапе для отработки схем распараллеливания и расчета выигрыша по времени используются две ЭВМ со схемой коммуникации типа "точка-точка", где обе машины являются равноправными.

Целью этой работы является создание программного интерфейса "PSS" (parallel solving for Shell) на основе MPI, позволяющего, не обращаясь к низкоуровневым командам MPI, производить "прозрачное" распараллеливание вычислений. Таким образом, общая блок-схема программы не изменится, за исключением того, что отдельные блоки в ней выполняются параллельно.

Ниже показана принципиальная блок-схема выполнения программы (рис.15). Блоки 2 и 6 можно производить в коммуникационной среде посредством нескольких процессоров, используя интерфейс PSS. Блок 3 можно распараллелить при использовании метода Гаусса для решения систем алгебраических уравнений.

Выигрыш по времени определяется следующими основными параметрами:

количество ЭВМ (процессоров) в коммуникационной среде;

отношение части кода программы, пригодной для распараллеливания, ко всему объему кода программы;

степень сложности логики проверки и контроля за вычислениями со стороны программного интерфейса PSS;

особенности аппаратной части (пропускная способность канала, скорость чтения/записи, частота процессоров, …).

87

Рис.15. Блок-схема выполнения программы

Если мы будем использовать N ЭВМ с одинаковой производительностью, то выигрыш по времени составит приблизительно (T/N) *1.3, где T - время выполнения программы на одном процессоре, а 1.3 - 30% -ая поправка, включающая в себя проверку и контроль со стороны PSS и особенности аппаратуры.

Все вышеуказанные положения особенно важны, принимая во внимание исследовательский характер комплекса "Оболочка", что подразумевает проведение большого количества вычислений с использованием различных параметров оболочек, нагружений, закреплений, ребер и т.д. Поэтому даже небольшое ускорение процесса вычисления одного эксперимента будет ощутимой экономией времени при многократном повторе этого эксперимента.

1. Антонов А.С. Параллельное программирование с использованием технологии MPI - М.: Изд-во МГУ, 2004. - 71 с.

2. Бондаренко В.М., Бондаренко С.В. Инженерные методы нелинейной теории железобетона. - М.: Стройиздат, 1982. - 288 с.

3. Васильев А.Н. Самоучитель C++ с задачами и примерами - М.: Наука и Техника, 2010. - 480 с.

4. Воеводин В.В. Вычислительная математика и структура алгоритмов - М.: Изд-во МГУ, 2006. - 112 с.

5. Воеводин В.В. Математические модели и методы в параллельных процессах - М.: Наука. Гл. ред. физ. - мат. лит., 1986. - 296 с.

6. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления - СПб.: БХВ-Петербург, 2002. - 608 с.

7. Воеводин Вл.В., Жуматий С.А. Вычислительное дело и кластерные системы - М.: Изд-во МГУ, 2007. - 150 с.

8. Жгутов В.М. Математические модели и алгоритмы исследования устойчивости пологих ребристых оболочек при учете различных свойств материала // Известия Орловского гос. техн. ун-та. Серия "Строительство, транспорт". - 2007. №4. - С. 20-23.

9. Ильюшин А.А. Пластичность. М.: Гостехиздат. 1948. - 376 с.

10. Карпов В.В. Математическое моделирование, алгоритмы исследования модели, вычислительный эксперимент в теории оболочек: учеб. пособие - СПб.: СПбГАСУ, 2006. - 330 с.

11. Карпов В.В., Баранова Д.А., Беркалиев Р.Т. Программный комплекс исследования устойчивости оболочек - СПб.: СПбГАСУ, 2009. - 102 с.

12. Карпов В.В., Сальников А.Ю. Вариационные методы и вариационные принципы механики при расчете строительных конструкций: учеб. пособие - СПб.: СПбГАСУ, 2009. - 75 с.

13. Керниган Б, Ритчи Д. Язык программирования C - М.: Вильямс, 2009. - 304 с.

14. Климанов В.И., Тимашев С.А. Нелинейные задачи подкрепленных оболочек. Свердловск: УНЦ АН СССР, 1985. - 291 с.

15. Петров В.В., Овчинников И.Г., Ярославский В.И. Расчет пластинок и оболочек из нелинейно-упругого материала. - Саратов: Изд-во Сарат. ун-та, 1976. - 136 с.

16. Хьюз К., Хьюз Т. Параллельное и распределенное программирование с использованием C++ - М.: Вильямс, 2004. - 672 с.

Приложение 1. Код программы

Main. cpp:

/*

Main. cpp

Сизов А.С. ПМ5

2010 г.

*/

// Подключение необходимых заголовочных файлов

#include "mpi. h" // библиотека mpi

#include <stdio. h>

/*

stdio. h (от англ. standard input/output header -

стандартный заголовочный файл ввода/вывода) заголовочный файл

стандартной библиотеки языка Си, содержащий определения макросов,

константы и объявления функций и типов, используемых для различных

операций стандартного ввода и вывода. Функциональность унаследована

от "портативного пакета ввода/вывода" ("portable I/O package"),

написанного Майком Леском из Bell Labs в начале 1970-х.

Функции, объявленные в stdio. h, являются весьма популярными благодаря тому,

что являясь частью Стандартной библиотеки языка Си, они гарантируют работу

на любой платформе, поддерживающей Си. Приложения на отдельных платформах

могут, тем не менее, иметь причины для использования функций ввода/вывода

самой платформы вместо функций stdio. h.

*/

#include <stdlib. h>

/*

stdlib. h - заголовочный файл стандартной библиотеки

общего назначения языка Си, который содержит в себе функции,

занимающиеся выделением памяти, контроль процесса выполнения программы,

преобразования типов и другие.

*/

#include <math. h>

/*

math. h - заголовочный файл стандартной библиотеки языка программирования С,

разработанный для выполнения простых математических операций. Большинство

функций привлекают использование чисел с плавающей точкой. Все эти функции

принимают double, если не определено иначе. Для работы с типами float и

long double используются функции с постфиксами f и l соответственно.

Все функции, принимающие или возвращающие угол, работают с радианами.

*/

#include "matrix. h"

/*

matrix. h - заголовочный файл библиотеки, содержащей определения

класса матрицы matrix, использованного в дальнейшем в программе.

Содержит объявление класса, методы для работы с ним, такие как:

произведение матрицы на матрицу;

сложение матрицы с матрицей;

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

выгрузка данных;

загрузка данных;

копирование матрицы;

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

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

*/

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

double *X1, *X2, *X3, *Y1, *Y2, *Y3; // Коэффициенты аппроксимирующих // функций

double *C1, *C2; // Коэффициенты системы

int N, n, m, E;

double q, h, mu, r, Pi, mu1, Kx, Ky, R1, R2, W;

double a,b;

/*

N - количество членов в приближении по методу Ритца;

n - количество ребер по x;

m - количество ребер по y;

E - модуль упругости материала;

q - поперечная равномерно-распределенная нагрузка;

h - толщина оболочки;

mu - коэффициент Пуассона;

r - ширина ребер;

Pi - число Пи;

R1 - радиус кривизны по x;

R2 - радиус кривизны по y;

a - длина оболочки по x;

b - длина оболочки по y;

*/

// Процедуры вычисления производных функций

double sin_0 (double k, double x)

{

return sin (k * x);

}

double sin_1 (double k, double x) // первая производная

{

return k * cos (k * x);

}

double sin_2 (double k, double x) // вторая производная

{

return - k * k * sin (k * x);

}

double Fx (int k, int i, int j, double arg) // выборка произведений

{ // аппроксимирующих

switch (k) // функций для х

{

case 1:

return sin_1 (X1 [i], arg) * sin_1 (X1 [j], arg);

case 2:

return sin_0 (X1 [i], arg) * sin_0 (X1 [j], arg);

case 3:

return sin_1 (X1 [i], arg) * sin_0 (X2 [j], arg);

case 4:

return sin_0 (X1 [i], arg) * sin_1 (X2 [j], arg);

case 5:

return sin_1 (X1 [i], arg) * sin_0 (X3 [j], arg);

case 6:

return sin_0 (X2 [i], arg) * sin_0 (X2 [j], arg);

case 7:

return sin_1 (X2 [i], arg) * sin_1 (X2 [j], arg);

case 8:

return sin_0 (X2 [i], arg) * sin_0 (X3 [j], arg);

case 9:

return sin_0 (X3 [i], arg) * sin_1 (X1 [j], arg);

case 10:

return sin_0 (X3 [i], arg) * sin_0 (X3 [j], arg);

case 11:

return sin_2 (X3 [i], arg) * sin_2 (X3 [j], arg);

case 12:

return sin_2 (X3 [i], arg) * sin_0 (X3 [j], arg);

case 13:

return sin_1 (X3 [j], arg) * sin_1 (X3 [j], arg);

case 14:

return sin_0 (X3 [j], arg);

default:

return 0;

}

}

double simpsonFx (double a, double b, int k, int i, int j)

{ // метод Симпсона для х

int n = 1000;

double locI, h, xi, s1 = 0.0, s2 = 0.0, y [10000] ;

int f;

h = (b - a) / n;

xi = a;

for (f = 0; f <= 2 * n; f++)

{

y [f] = Fx (k, i,j,xi);

xi = xi + h/2;

}

locI = y [0] + y [2*n] ;

for (f = 2; f < 2*n; f = f + 2)

s1 = s1 + y [f] ;

for (f = 1; f < 2*n; f = f + 2)

s2 = s2 + y [f] ;

locI = (locI + 2.0 * s1 + 4.0 * s2) * (h / 6);

return (locI);

}

double Fy (int k, int i, int j, double arg) // выборка произведений

{ // аппроксимирующих

switch (k) // функций для y

{

case 1:

return sin_0 (Y1 [i], arg) * sin_0 (Y1 [j], arg);

case 2:

return sin_1 (Y1 [i], arg) * sin_1 (Y1 [j], arg);

case 3:

return sin_0 (Y1 [i], arg) * sin_1 (Y2 [j], arg);

case 4:

return sin_1 (Y1 [i], arg) * sin_0 (Y2 [j], arg);

case 5:

return sin_0 (Y1 [i], arg) * sin_0 (Y3 [j], arg);

case 6:

return sin_1 (Y2 [i], arg) * sin_1 (Y2 [j], arg);

case 7:

return sin_0 (Y2 [i], arg) * sin_0 (Y2 [j], arg);

case 8:

return sin_1 (Y2 [i], arg) * sin_0 (Y3 [j], arg);

case 9:

return sin_0 (Y3 [i], arg) * sin_0 (Y1 [j], arg);

case 10:

return sin_0 (Y3 [i], arg) * sin_0 (Y3 [j], arg);

case 11:

return sin_2 (Y3 [i], arg) * sin_2 (Y3 [j], arg);

case 12:

return sin_2 (Y3 [i], arg) * sin_0 (Y3 [j], arg);

case 13:

return sin_1 (Y3 [i], arg) * sin_1 (Y3 [j], arg);

case 14:

return sin_0 (Y3 [j], arg);

default:

return 0;

}

}

double simpsonFy (double a, double b, int k, int i, int j)

{ // метод Симпсона для y

int n = 1000;

double locI, h, xi, s1 = 0.0, s2 = 0.0, y [10000] ;

int f;

h = (b - a) / n;

xi = a;

for (f = 0; f <= 2 * n; f++)

{

y [f] = Fy (k, i,j,xi);

xi = xi + h/2;

}

locI = y [0] + y [2*n] ;

for (f = 2; f < 2*n; f = f + 2)

s1 = s1 + y [f] ;

for (f = 1; f < 2*n; f = f + 2)

s2 = s2 + y [f] ;

locI = (locI + 2.0 * s1 + 4.0 * s2) * (h / 6);

return (locI);

}

double C (int k, int i, int j) // коэффициэнты в ФПЭД

{

switch (k)

{

case 1:

return 2*h*simpsonFx (0.0,a,1, i,j) *simpsonFy (0.0,a,1, i,j) +2*mu1*h*simpsonFx (0.0,a,2, i,j) *simpsonFy (0.0,a,2, i,j);

case 2:

return 2*h*mu*simpsonFx (0.0,a,3, i,j) *simpsonFy (0.0,a,3, i,j) +2*mu1*h*simpsonFx (0.0,a,4, i,j) *simpsonFy (0.0,a,4, i,j);

case 3:

return - 2*h* (Kx+mu*Ky) *simpsonFx (0.0,a,5, i,j) *simpsonFy (0.0,a,5, i,j);

case 4:

return 2*h*mu*simpsonFx (0.0,a,3, i,j) *simpsonFy (0.0,a,3, i,j) +2*mu1*h*simpsonFx (0.0,a,4, i,j) *simpsonFy (0.0,a,4, i,j);

case 5:

return 2*h*simpsonFx (0.0,a,6, i,j) *simpsonFy (0.0,a,6, i,j) +2*mu1*h*simpsonFx (0.0,a,7, i,j) *simpsonFy (0.0,a,7, i,j);

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) (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