«Н.Н ОЛЕНЁВ, Р.В ПЕЧЁНКИН, А.М. ЧЕРНЕЦОВ ПАРАЛЛЕЛЬНОЕ ПРОГРАММИРОВАНИЕ В MATLAB И ЕГО ПРИЛОЖЕНИЯ ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР ИМ. А.А ДОРОДНИЦЫНА РАН МОСКВА 2007 УДК 512.643 + 519.86 Ответственный редактор академик РАН А.А. ...»

РОССИЙСКАЯ АКАДЕМИЯ НАУК

ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР ИМ. А.А. ДОРОДНИЦЫНА

Н.Н ОЛЕНЁВ, Р.В ПЕЧЁНКИН, А.М. ЧЕРНЕЦОВ

ПАРАЛЛЕЛЬНОЕ ПРОГРАММИРОВАНИЕ В

MATLAB И ЕГО ПРИЛОЖЕНИЯ

ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР ИМ. А.А ДОРОДНИЦЫНА РАН

МОСКВА 2007

УДК 512.643 + 519.86

Ответственный редактор

академик РАН А.А. Петров Работа посвящена вопросам использования распределенных и параллельных вычислений в среде MATLAB. Изложена технология настройки кластера для использования MATLAB. Технология апробирована на стандартных задачах линейной алгебры. Предложен промышленный подход идентификации математических моделей сложных систем на основе параллельных средств MATLAB.

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

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (коды проектов 07-01-00563-а, 07-01-12032-офи), Российского гуманитарного научного фонда (код проекта 06-02-91821-а/G), программы Intel Education Project 2006-2007, программы Президиума РАН № 15, программы ОМН РАН № 3, по программе государственной поддержки ведущих научных школ (код проекта НШ-5379.2006.1).

Рецензенты Г.М. Михайлов, В.С. Молоствов Научное издание c Вычислительный центр им. А.А. Дородницына Российской академии наук, Оглавление Введение I. Администрирование и конфигурирование 1.1 Используемые термины и обозначения..... 1.2 Требования к аппаратному и программному обеспечению.................... 1.3 Инсталляция и запуск MDCE.......... 1.3.1 Работа в Windows............. 1.3.2 Работа в Unix................ 1.4 Запуск планировщика............... 1.4.1 Запуск в Windows............. 1.4.2 Запуск в Unix............... 1.5 Запуск рабочих процессов............ 1.6 Объекты-ссылки на системные процессы.... 1.7 Эквиваленты между функциями DCT и MPI. 1.8 Использование сторонних планировщиков... 1.9 Заключение..................... II. Программирование параллельных задач 2.1 Режим pmode, первая параллельная задача.. 2.2 Вычисление определенного интеграла...... 2.3 Объект ”параллельная задача”.......... 2.3.1 Объект ”параллельное задание”..... 2.3.2 Матричное умножение.......... 2.3.3 Решение СЛАУ методом Гаусса..... III. Параллельные вычисления в моделировании экономики 3.1 Проблема идентификации внешних параметров модели экономики.......... 3.2 Численная реализация задачи идентификации. 3.2.1 Листинги программ............ 3.2.2 Результаты идентификации модели, их графическое представление........ 3.3 Сценарные расчеты с моделью.......... 3.3.1 Базовый, он же пессимистический сценарий.................. 3.3.2 Оптимистический сценарий....... 3.3.3 Сравнение, выводы............ Заключение Введение Данная работа представляет собой первую попытку в виде научной монографии, доступной начинающему пользователю MATLAB, описать параллельную технологию вычислений, реализованную компанией MathWorks с помощью двух взаимосвязанных пакетов расширений (приложений): MATLAB Distributed Computing toolbox [1] и MATLAB Distributed Computing Engine [2].

В настоящее время не только компания MathWorks обладает технологиями, позволяющими реализовывать параллельные вычисления. Следует упомянуть об аналогичных средствах, имеющихся в математических пакетах других фирм. Так, Maple, начиная с версии 9.5, имеет в своем составе toolbox HPC-Grid [3], который предлагает пользователю этой среды аналогичные возможности. Для пакета Mathematica имеется Parallel Computing Toolkit [4]. Также имеется независимый toolkit Pooch MPI [5], который работает только под операционной системой MacOS, что резко снижает его область применения. Часто применяемый пакет Mathcad на момент написания работы не имеет каких-либо средств, позволяющих выполнять распределенные вычисления.

Несмотря на наличие подобного рода инструментов в пакетах компьютерной алгебры Maple и Mathematica, среда технических расчетов MATLAB является, по мнению авторов, более привычной и удобной для программирования, отладки и реализации параллельных алгоритмов.

За последние годы различными научными коллективами были созданы многочисленные пакеты расширений (Toolbox) MATLAB, реализующих какую-либо функциональность распределенных вычислений. С их обзором можно ознакомиться в [6]. Эти пакеты в своей реализации использовали различные способы (парадигмы):

• Использование модели передачи сообщений.

• Работа с общей памятью.

• Использование специальных процедур, с помощью которых производится обмен данных между последовательными сессиями MATLAB.

• Перекомпиляция кода MATLAB на С и создание параллельного кода уже для программ С.

В зависимости от реализации при использовании этих пакетов расширений требовалось иметь разное количество лицензий на запуск клиентских сессий MATLAB.

Пакеты расширений для параллельного программирования сторонними производителями начали разрабатываться в то время, когда необходимость в параллельном программировании для MATLAB не являлась критической. С тех пор в связи с необходимостью расчета более сложных задач эта необходимость резко возросла.

Некоторые toolbox работали только под определенную операционную систему, чаще всего – под Unix. Поэтому использование единого подхода для всех операционных систем стало необходимым. Большая часть разработанных пакетов уже не поддерживается. Особо следует отметить пакет MPITB [7], написанный в 2000 г. и получивший достаточно широкое распространение для кластерных систем.

В нем реализовано применение функций Message Passing Interface – интерфейса передачи сообщений MPI – в программах MATLAB. Последняя на момент написания работы версия - для MATLAB 7.0.1 R14SP1 и библиотеки LAM MPI v. 7.1.1 [8] - датируется 2005 г. Основными его недостатками являются:

• Необходимость получения лицензии на каждое выполняемое ядро.

• Невозможность применения пакета для неоднородных вычислительных архитектур.

• Отсутствие официальной поддержки от компании MathWorks, постоянно развивающей систему MATLAB.

• Исключение работы пакета с операционными системами, отличающимися от ОС Unix.

Данный пакет был импортирован в среду Windows и реализацию MPICH2 [9]. Его дистрибутивы можно скачать с [10].

Последняя версия от 02.2006 работает с MATLAB v. 7.1.0. (R14) Service Pack 3, т.е. также с устаревшей версией MATLAB.

Учитывая все перечисленное выше, реализация подобных пакетов расширений уже самой фирмой MathWorks выглядит вполне логичным шагом.

Авторы ставили перед собой цель написать такую работу, которая могла бы значительно ускорить процесс освоения технологии параллельного программирования читателем, который вообще не знаком ни с навыками параллельного программирования в системе MATLAB, ни в целом с системой.

В первом разделе описаны основные процедуры установки, настройки и конфигурирования кластерной части MATLAB Distributed Computing Engine. Этот раздел будет интересен в первую очередь системным администраторам и администраторам кластерных систем. В данном разделе авторы попытались представить информацию, содержащуюся "между строчек" в стандартной технической документации. Как показал опыт сотрудничества с другими исследовательскими группами, информации, представленной в этом разделе, достаточно для настройки высокопроизводительной вычислительной среды.

Во втором разделе описаны первоначальные сведения для создания параллельных программ с использованием средств библиотеки MPI и ее реализации на основе платформы MATLAB. Рассмотренные простые задачи позволяют пользователю убедиться в приросте производительности вычислений при увеличении числа процессоров и понять базовые принципы написания параллельных программ в MATLAB.

В третьем разделе авторы в качестве примера приложения параллельных вычислений в MATLAB для научных исследований рассматривают их применение в математическом моделировании экономических систем. Это приложение ни в коей мере не ограничивает применение параллельных вычислений в математическом моделировании сложных систем только данной областью науки. Во-первых, рассмотренные примеры легко могут быть обобщены на аналогичные задачи в других областях знаний. Во-вторых, экономические приложения к современной жизни интересны для всех любознательных людей.

Рассмотрение приложений параллельных вычислений очень важно, поскольку в связи с развитием многопроцессорных и многоядерных [11] архитектур вычислительных систем в настоящее время параллельное программирование планируют применять повсеместно во всех областях науки, техники и бизнеса. В настоящее время фактическим стандартом параллельного программирования служит довольно сложный для первоначального изучения интерфейс передачи сообщений MPI – библиотека передачи сообщений, собрание функций, облегчающих коммуникацию (обмен данными и синхронизацию задач) между процессами параллельной программы с распределенной памятью [12, 13]. Качественные реализации MPI обеспечивают асинхронную коммуникацию, эффективное управление буфером сообщения, эффективные группы и богатые функциональные возможности. MPI включает большой набор коллективных операций коммуникации, виртуальных топологий и различных способов коммуникации и, кроме того, MPI поддерживает неоднородные сети.

Приложения, разработанные компанией MathWorks для создания параллельных и распределенных программ с использованием средств библиотеки MPI, и их реализация на платформе MATLAB упрощают практическое применение параллельных вычислений на многоядерных компьютерах, кластерах и GRID-системах.

Авторы выражают искреннюю благодарность департаменту MathWorks компании Softline за предоставленные дистрибутивы MATLAB и возможность ознакомиться и протестировать технологические возможности, заложенные в приложении MATLAB Distributing Computing Engine.

I. Администрирование и Изучение настоящего раздела начинающим пользователем MATLAB даст возможность ощутить на собственном опыте преимущества решения практических задач с помощью распределенных вычислений на многопроцессорных, многоядерных или даже неоднородных архитектурах в сравнении с обычным (последовательным) вариантом решения этих задач при использовании одного компьютера.

1.1 Используемые термины и обозначения Цель данного раздела описать принципиальную схему работы двух Toolbox MATLAB, которые реализуют технологию распределенных и параллельных вычислений (в дальнейшем, если не будет оговорено отдельно, под термином распределенные вычисления будут пониматься как параллельные, так и распределенные вычисления).

В 2005 г. для того, чтобы занять определенный сегмент рынка инженерного программного обеспечения, предназначенного для программирования распределенных задач, компанией MathWorks были разработаны и выпущены на рынок два новых продукта под названием Distributing Computing Toolbox (DCT) и MATLAB Distributing Computing Engine (MDCE). Эти два пакета расширений нельзя рассматривать отдельно друг от друга, и применяются они только в связке. Оба эти toolbox изначально были предназначены для решения одной задачи – увеличения производительности (сокращения времени счета) при использовании нескольких компьютеров, объединенных в сеть, вместо одного.

Для дальнейшего изложения потребуется дать несколько определений основным понятиям, которые встретятся нам на пути освоения параллельной технологии.

Client – сессия MATLAB, в которой определяются и из которой отправляются задачи. Это сессия MATLAB, в которой обычно работает разработчик– программист. Также распространено название MATLAB-клиент.

Client computer – компьютер, на котором запущен MATLABклиент.

Cluster – набор компьютеров, которые соединены в сеть с целью решать общую вычислительную задачу..

Head node – обычно это узел кластера, предназначенный для запуска job manager и license manager. Часто бывает полезным запускать все сервисные процессы для распределенных вычислений, не связанные с рабочими процессами, на одном компьютере.

Coarse-grained application – приложения, для которых время исполнения (счета) значительно больше, чем время, требуемое для создания, отправки и получения данных расчетной задачи.

Distributed application – приложения, аналогичные coarse-grained application, запускаются независимо на различных узлах, возможно с различными входными параметрами. Не совершают обмен сообщениями, не имеют общих данных или точек синхронизации между узлами. Распределенные вычисления могут быть как coarse-grained, так и fine-grained (разветвленные задачи).

Distributed computing – вычисления, производимые с помощью распределенных приложений, запускаемых на нескольких узлах одновременно.

Job – полностью описанная вычислительная операция большой размерности для выполнения в MATLAB, состоящая из набора подзадач.

License Manager – компонент сервера лицензий, который отвечает за проверку и поддержку имеющихся лицензий на ядро/различные toolbox.

Task – подзадача, некоторый сегмент основной задачи.

Именно подзадача отправляется на счет worker.

Mdce – служба, которая должна быть запущена на всех машинах, перед тем, как там будут запущены job manager или worker-процессы. Это основополагающая среда для запуска всех остальных процессов.

Job Manager – фирменный планировщик The MathWorks, это процесс, который образует систему очередей из задач (job) и распределяет подзадачи для workers. Под этим термином понимаются также планировщики, разрабатываемые сторонними фирмами производителями.

Worker – системный процесс MATLAB, который выполняет вычисление подзадач. Также распространенное название – MATLAB worker или worker-процесс.

Scheduler – планировщик. Основной функционал заключается в процессе планировки и управления очередями задач, В качестве планировщика может выступать MathWorks job manager или продукты сторонних фирм разработчиков.

По мере необходимости мы будем вводить дополнительные определения, пока введенных терминов должно быть достаточно, чтобы рассмотреть простейшую схему (см. 1.1), положенную в основу Distributing Computing Toolbox. При использовании промышленного кластера в качестве MATLAB client может выступать любой компьютер, подключенный к сети и имеющий доступ к кластеру, для scheduler обычно выделяют один из узлов кластера (в принципе, на этом же узле можно запустить и процесс worker, это не должно сильно сказаться на производительности), на всех остальные узлах кластера можно запустить по одному процессу. Запускать на одном узле несколько процессов имеет смысл только в случаях, когда узел многопроцессорный либо многоядерный.

Детали процесса запуска процессов job manager и worker, а также небольшие ухищрения (описанные в официальной документации [2]) мы рассмотрим, когда перейдем к описанию соответствующего этапа.

Distributing Computing Job Manager Рис. 1.1. Простейшая конфигурация, реализующая механизм распределенных вычислений MATLAB Для первоначального знакомства с распределенной технологией с объектами, их свойствами и методами, службами и сервисами необязательно иметь в своем распоряжении промышленный кластер. Более того, для первоначального знакомства достаточно будет персонального компьютера, рабочей станции или ноутбука. Для того чтобы почувствовать прирост производительности, разумеется, необходим многопроцессорный комплекс, состоящий минимум из двух компьютеров. В его роли может выступать как обычный многоядерный процессор, так и многопроцессорная система независимо от способа организации - кластер или SMP.

Для инсталляции и работы DCT компьютер должен удовлетворять нескольким условиям. Во-первых, на нем возможно запустить MATLAB (разд. Требования к аппаратному и программному обеспечению). Во-вторых, накладываются определенные требования к сети (разд. Требования к сети).

1.2 Требования к аппаратному и программному Аппаратная платформа (общие требования):

• CD(DVD) дисковод (для инсталляции пакета).

• Для некоторых типов лицензий FlexLm Manager версии 10.8.0.1 (поставляется вместе с MATLAB).

• В случае использования Flexlm требуется применение протокола TCP/IP на всех платформах, использующих сервер лицензий.

• Для лицензий, использующих аппаратный ключ, требуется наличие USB-порта.

Минимальные требования к объему оперативной памяти и дисковому пространству изложены в таблице ниже.

Дисковое пространство 460 Мб Минимальное дисковое пространство 460 Мб требуется для инсталляции ядра и справочной системы.

Для 32-битной архитектуры возможна работа на следующих микропроцессорах фирм Intel и AMD.

Фирма Микропроцессор Intel Pentium III, Pentium 4, Pentium Xeon, Pentium M AMD Athlon, Athlon MP, XPAMD, Opteron Для 64-битной архитектуры можно использовать микропроцессоры: у Intel: EM64T, у AMD: AMD64.

Операционная система:

• Для Windows - Windows 2000 SP3/XP SP1/2003/Vista.

• Для Linux: начиная с ядра версии 2.4.x и GLIBC 2.3.2.

• Для Solaris : Solaris 8,9,10 для SPARC/UltraSPARC.

• Для Macintosh: powerPC G4 - Mac OS X 10.4.7, intel Mac OS X 10.4.8.

В данной работе применение Solaris и Macintosh не рассматривается по нескольким причинам:

• эти системы значительно менее распространены, чем Windows и Linux, • для достижения поставленных целей описания технологии распределенного программирования от компании MathWorks вполне достаточно рассмотреть две системы, • операции для Solaris аналогичны операциям для Общие требования к аппаратной платформе и операционной системе определены. Это позволяет запускать ядро MATLAB. Для запуска Distributing Computing Toolbox накладываются также опеределенные требования к сети.

• Распределенные вычислительные процессы должны быть способны определять друг друга по именам хостов, что требует наличия в сети службы разрешения имен (DNS, например). Кроме того, имя хоста, под которым компьютер виден остальным компьютерам сети, должно совпадать с именем хоста, под которым компьютер определяет сам себя.

• Имя хоста для каждого узла кластера должно ограничиваться IP-адресом, которому соответствует адрес одной из сетевых карт. Также поддерживаются машины с несколькими сетевыми картами.

• Рекомендуется минимум 5 ГБ дискового пространства для хранения рабочих файлов при использовании встроенного планировщика JobManager.

• Распределенные вычислительные процессы используют несколько TCP портов. Для узла, на котором запущено всего n JobManager и рабочих процессов, служба MDCE резервирует для собственного использования 5+n последовательных портов.

• Дополнительные TCP порты открываются для межпроцессного обмена между рабочими процессами во время выполнения параллельного задания.

• Чтобы запускать параллельные приложения, которые используют межпроцессный обмен, все узлы кластера, выполняющие приложение, должны иметь одинаковый размер слов и архитектуру команд процессора. Гетерогенные кластерные системы, на которых возможен запуск - Solaris и Macintosh. Системы Windows и Linux в версии MDCE 3.0 не поддерживаются.

MATLAB 2006b поставляется с реализацией MPI-2 [13] mpich2 v. 1.0.3. Вместо стандартной библиотеки можно подключить собственную реализацию. Для этого она должна удовлетворять следующим условиям:

• Быть собранной как динамическая библиотека.

• Поддерживать все функции MPI-1.

• Поддерживать пустые аргументы в MPI_Init, в соответствии с разд. 4.2 стандарта MPI-2.

• Иметь заголовочный файл mpi.h, полностью совместимый с mpich2.

• Для работы встроенного JobManager от MathWorks дополнительно требуется поддержка функций MPI- MPI_Open_port, MPI_Comm_accept, MPI_Comm_connect.

Цель, которую мы ставим в данном разделе - это запустить на нашей рабочей станции службу, которая, в свою очередь, позволит запустить процесс job manager. После запуска с его помощью рабочего процесса станет возможным на одном компьютере ознакомиться с основными свойствами и методами объектов.

Во-первых, необходимо удостовериться, что на вашем компьютере установлены соответствующие toolbox. Для этого, запустив MATLAB, в командном окне введите команду >>ver. В списке проинсталированных toolbox обязательно должны быть следующие два toolbox:

Distributed Computing Toolbox Version 3.0(R2006b), MATLAB Distributed Computing Engine Version 3.0(R2006b).

Остальные пакеты расширений в данный момент нас не интересуют.

В зависимости от того, как у вас установлен дистрибутив MATLAB, С:\MATLAB\R2006B или С:\MATLABR2006B или еще как-то, будем называть первую часть пути к директории MATLAB просто путь – path.

В каталоге path\toolbox\distcomp\bin\ содержатся *.bat файлы, которые позволяют запускать служебные сервисы и процессы. Для запуска службы mdce пользователь должен либо в окне сеансе cmd, либо непосредственно в MATLAB запустить файл mdce.bat install.

Например, в приведенном примере, когда путь имел вид С:\MATLAB\R2006B или С:\MATLABR2006B запуск (инсталляция) службы mdce непосредственно из командной строки MATLAB выглядела бы следующим образом:

!C:\MATLAB\R2006B\toolbox\distcomp\bin\mdce install или >>!C:\MATLAB\R2006B\toolbox\distcomp\bin\mdce install знак ! указывает MATLAB, что выполняемая команда системная.

В обоих случаях результат должен быть одним и тем же, в списке служб Windows должна появится новая служба – MATLAB Distributing Computing Engine.

Аналогичный результат может быть получен с использованием функции eval(coommand). В качестве параметра command необходимо передать строку, составив ее так, как если бы мы вводили команду в сеансе командной строки.

Корневой каталог установки MATLAB можно определить, например, следующим образом:

>>command1=matlabroot command1 = C:\MATLAB Затем создадим переменную command2:

>> command2=’\toolbox\distcomp\bin\’ Теперь склеим наши переменные:

>> command=[’!’ command1 command2 ’mdce install’];

И передадим функции eval значение command >> eval(command) Аналогично будем поступать и в других случаях, когда потребуется запустить другие *.bat файлы.

В том случае, если вам требуется получить быструю справку по какой-либо команде, достаточно ввести >> help function_name в этом случае будет выведен листинг help информации, расположенной в заголовке функции command, или ввести >> doc function_name в этом случае откроется страница помощи, на которой будет отображена более чем подробная информация по данной команде.

В том случае, если вы забыли, как точно называется требуемая функция, можно, набрав только первые буквы, нажать кнопку "Tab", и перед вами откроется ниспадающий список с функциями, начинающимися с введенных ранее букв. Теперь можно быстро, Рис. 1.3. Ниспадающая, быстрая подсказка перемещаясь по списку, найти нужную команду (см. рис 1.3).

Инсталляция MDCE в Unix происходит практически аналогично с Windows. Настройки MCDE в отличие от Windows находятся в файле mdce_def.sh.

Методы лицензирования в средах Windows и Unix отличаются.

Если в Windows применяется пара лицензионный ключ + аппаратный ключ, то в Unix используется механизм лицензионных файлов и широко известный Flexlm Manager. В лицензионном файле, в частности, обязательно должны быть прописаны строки, относящиеся к Distributing Computing toolbox и MDCE, например:

INCREMENT Distrib_Computing_Toolbox MLM 16 10-feb-2007 0 \ CDE6E0ABF6779416F1F4 HOSTID=DEMO SN=DEMO # LicenseNo: DEMO HostID:

INCREMENT TMW_Archive MLM 16 01-jan-0000 0 \ CD86C100118CC4209A0F \ VENDOR_STRING=908ff75f02aae4679ac61e7ccff HOSTID=DEMO SN= INCREMENT MATLAB_Distrib_Comp_Engine MLM 16 15-feb-2007 16 \ 9DE651D0F3651D1D608B HOSTID=ID=0 SN=DEMO Покажем последовательность действий для инсталляции mdce. Пусть MATLABROOT - корневой каталог, в который установлен MATLAB. Отметим, что при работе в конфигурации по умолчанию установка и запуск mdce требует прав суперпользователя root.

[root@broody bin]# cd $MATLABROOT/toolbox/distcomp/bin [root@broody bin]#./mdce start Creating LOGBASE directory (/var/log/mdce).

Creating CHECKPOINTBASE directory (/var/lib/mdce).

Starting the MATLAB Distributed Computing Engine in the background.

Признаком нормальной работы службы mdce является следующее сообщение при запросе статуса:

[root@broody bin]#./mdce status The MATLAB Distributed Computing Engine is running with PID 13592.

Use nodestatus to obtain more information.

Обратимся теперь ко второму этапу конфигурирования, а именно к процедуре запуска процесса Job Manager. C помощью файла startjobmanager.bat, который расположен в директории MATLAB\toolbox\distcomp\bin\, производится запуск системного процесса, который будет по сути нашим планировщиком. Запуск можно произвести только в том случае, если служба mdce уже запущена.

Сначала рассмотрим использование встроенного планировщика от MathWorks - JobManager. Об использовании других планировщиков можно прочесть в разд. 1.7.

Команде запуска startjobmanager можно передавать различные параметры с использованием ключей. Приведем описание только некоторых ключей. Полное описание пользователь может найти в стандартной документации.

-name – имя Job Manager, идентификатор. Если ключ name не задан, то процессу присвоится значение переменной DEFAULT_JOB_MANAGER_NAME=default_jobmanager, которое определено в файле MATLAB\toolbox\distcomp\ bin\mdce_def. Пользователь может сам поменять значения переменных в этом файле на свое усмотрение.

-remotehost – имя удаленной рабочей станции, узла, на котором будет произведен запуск этого процесса. Обязательным требованием является запуск службы mdce на этой удаленной сессии. В качестве значения этого ключа можно передавать IP адрес компьютера (узла). Если ключ не используется, процесс будет запущен на локальной станции (с которой производится запуск).

-v – режим verbose, подробного листинга, во время запуска отображаются системные сообщения.

При запуске в среде Unix синтаксис команды startjobmanager совпадает с запуском в Windows. В приведенном ниже примере производится запуск JobManager с именем MyJobManager на локальной машине. Включен режим выдачи отладочной информации.

[root@broody bin]#./startjobmanager -name MyJobManager -v Contacting the mdce service on the host broody.ccas.ru to start the job manager lookup process.

Started the job manager lookup process on the host broody.ccas.ru.

Contacting the mdce service on the host broody.ccas.ru to start a job manager process.

Started the job manager MyJobManager on the host broody.ccas.ru.

Именно системные процессы, называемые workers или "удаленные сессии MATLAB", обеспечивают решение подзадач (task).

Запуск worker аналогичен запуску процесса JobManager. С помощью файла startworker.bat, который так же, как и в предыдущем случае, запускается с использованием ключей -name – имя MATLAB worker. Если ключ не задан, то имя по умолчанию это параметр DEFAULT_WORKER_NAME заданный в файле mdce_def file.

-jobmanager – значением этого ключа должно быть имя jobmanager, с которым будет ассоциирован (от которого будет получать подзадачи) запускаемый процесс worker. По умолчанию присваивается значение параметра DEFAULT_JOB_MANAGER_NAME, которое задано в файле mdce_def file.

-jobmanagerhost – IP адрес или имя компьютера на котором запущен процесс jobmanager.

-remotehost – имя удаленного компьютера (IP адрес) на котором должен быть запущен данный процесс worker.

-v – режим verbose, подробного листинга, во время запуска отображаются системные сообщения.

В нашем примере необходимо выполнить команду >>!C:\MATLAB\R2006B\toolbox\distcomp\bin \startworker -name worker1 -remotehost mylaptop -jobmanager MyJobManager -jobmanagerhost mylaptop, после чего будет запущен процесс с именем worker1.

После того как в сети (на кластере) запущен хотя бы один процесс jobmanager, можно использовать команду nodestatus.

Данная команда позволяет получать информацию с различной детализацией о запущенных на конкретном узле (node):

• процессах, • Job Manager, • службах MDCE.

Синтаксис этой команды очень прост. Имеется два ключа для передачи параметров этой процедуре:

-remotehost – имя удаленного компьютера (IP адрес), на котором должен быть запущен данный процесс worker, -infolevel – уровень информативности, 1 - минимальная информация, 3 - подробная.

Выполнив команду >>!C:\MATLAB\R2006B\toolbox\distcomp\bin\nodestatus -remotehost mylaptop -infolevel при одном запущенном jobmanager и одном процессе worker, ассоциированном с этим планировщиком, вы увидите следующее сообщение Job manager lookup process:

Job manager:

Worker:

Summary:

The mdce service on mylaptop manages processes:

Поменяв значение infolevel на 3, пользователь может увидеть более детальную информацию о процессах, которые запущены на узле. Выводимой информации достаточно, чтобы выявить некорректные моменты в работе кластера, например выявить запущенные рабочие процессы, не ассоциированные с JobManager.

Для того чтобы завершить создание нашего минимального кластера, нам необходимо запустить процесс worker на втором компьютере, который соединен через сеть с первым (служба mdce должна быть запущена и на этом компьютере тоже). Пусть имя второго worker будет worker2. Для запуска worker2 на втором компьютере, нужно в команде startworker в качестве ключа remotehost указать имя второго компьютера (пусть имя компьютера будет mydesktop).

Теперь можно считать, что минимальный кластер создан.

1.6 Объекты-ссылки на системные процессы Итак, как уже было сказано, наши процессы jobmanager и worker живут отдельно от клиентской сессии MATLAB. При остановленном клиентском приложении MATLAB эти процессы тем не менее продолжают существовать в системе.

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

В этом случае в рабочей области MATLAB мы получим перемен ную – объект распределенных вычислений – и с помощью ее свойств и методов сможем управлять этими объектами.

Для поиска ссылки на процесс jobmanager необходимо воспользоваться функцией ndResource (данная функции принадлежит MATLAB Distributing Computing toolbox) в качестве параметров передать тип искомого процесса и имя узла.

jm = findResource(’scheduler’,’type’,’jobmanager’,...

’Name’,’MyJobManager’,’LookupURL’,’mynote’) Результатом выполнения данной команды будет переменная jm в текущей области MATLAB.

distcomp.jobmanager Посмотреть свойства объекта можно, как и в обычном случае, воспользовавшись графическим режимом, дважды щелкнув курсором по объекту jm или воспользовавшись командой get(jm), которая является предопределенной для всех классов системы MATLAB.

>>get(jm) NumberOfBusyWorkers: NumberOfIdleWorkers: Вообще говоря, методы, которыми обладает какой-либо объект, можно узнать с помощью команды method. Так, например, передав в качестве параметра объект jm, мы можем узнать, какими методами он обладает.

>> methods(jm) Methods for class distcomp.jobmanager:

createParallelJob findJob promote Ссылка на рабочий процесс worker Аналогично c помощью функции findResource мы можем найти ссылку на системный процесс исполнителя - worker.

>> worker = findResource(’worker’,’LookupURL’,’mynote’) distcomp.worker >> get(worker) HostAddress: {0x1 cell} JobManager: [1x1 distcomp.jobmanager] Вообще говоря, получить объект worker можно и с помощью объекта jobmanager. Так, у объекта (переменной) jm есть свойство IdleWorkers (число свободных рабочих процессов), значение которого в нашем случае [1x1 distcomp.worker]. Для того чтобы получить объект worker, обратимся через точку к соответствующему свойству >> w=jm.IdleWorkers >> get(w) JobManager: [1x1 distcomp.jobmanager] В данном случае переменная w будет равносильна переменной worker – обе эти переменные представляют собой ссылки на один и тот же системный процесс worker.

Следует заметить, что эти объекты – ссылки на системные процессы – обладают занимательным свойством, а именно:

размер этих переменных в байтах равен нулю.

>> whos jm Grand total is 1 element using 0 bytes >> whos worker Grand total is 1 element using 0 bytes Несмотря на то, что сам системный процесс jobmanager использует 512 Мб оперативной памяти, JOB_MANAGER_MAXIMUM_ MEMORY = 512 m - как это определено в файле mdce_def.bat, процессу worker доступно всего 60 Мб. Каким образом можно изменить этот параметр - размер процессов, можно узнать, воспользовавшись стандартной документацией [2]. В нашем сквозном примере эти параметры изменятся не будут.

В том случае, если требуется программно запустить несколько рабочих процессов на нескольких узлах, можно написать цикл, в котором будет происходить запуск соответствующего процесса worker на соответствующем узле. Причем текст команды будет формироваться (склеиваться) программно.

mroot=matlabroot;

% путь установки MATLAB catalog=’\toolbox\distcomp\bin\win32\’;

% каталог, в котором установлен toolbox mdcepath=[mroot catalog] % соединяем две переменные (склеиваем строку) nodes=[{’ip1’},{’ip2’},{’ip3’},{’ip4’},{’ip5’}];

% массив ячеек - каждая из которых IP адрес или имя узла MyJobManager = ’MyJobManager’;

% имя уже запущенного jobmanager MyJMhost = ’mynote’;

% имя узла, на котором запущен jobmanager for i=1:length(ips) % цикл для запуска worker на соответствующем узле name=[’worker’ num2str(i)];

% уникальное имя worker command=[’!’ mdcepath ’startworker -name ’ name...

% склеили команду, которая запустит i-й процесс eval(command) % функция eval выполняет команду, % которую ей передали в виде строки - входного параметра, end Результатом работы этого участка кода являются запущенные процессы на соответствующих узлах. Аналогичный код можно написать, например, и для команд stopworker, nodestatus.

Служба MDCE, а также планировщик и рабочие процессы worker, генерируют различные рабочие файлы в процессе своей работы. Рабочие файлы бывают двух типов: файлы журнала и файлы контрольных точек (checkpoints).

Windows "TEMP\MDCE\log" TEMP\MDCE\Checkpoint Unix /var/log/mdce /var/lib/mdce Еще раз отметим, что согласно данным производителя, для хранения файлов на жестком диске необходимо иметь не менее 5 Гбайт свободного места в рабочих каталогах. Также следует иметь в виду, что после сбоя работа рабочих процессов и jobmanager начнется с момента последней контрольной точки.

Чтобы этого избежать, необходимо использовать при запуске ключ -clean.

В процессе работы размер рабочего каталога (/var/lib/mdce или TEMP\MDCE\Checkpoint) стремительно растет. Чтобы избежать переполнения диска, мы рекомендуем выполнять после каждого запуска уничтожение созданного задания через функцию destroy(job).

1.7 Эквиваленты между функциями DCT и MPI В приложении Distributed Computing toolbox (DCT) системы MATLAB версии 2006b имеется встроенный набор функций для передачи сообщений между процессами, основанный на интерфейсе передачи сообщений Message Passing Interface (MPI). Эквиваленты между функциями DCT системы MATLAB и библиотеки MPI [12] для языка программирования С/С++ приведены в табл. 1.1.

Таблица 1.1. Сравнение основных функций передачи сообщений в DCT MATLAB с соответствующими функциями в MPI Функция/ переменная MATLAB Функция/ переменная MPI в С/С++ shared_data=... MPI_Bcast(void* buffer, int count, LabBroadcast(root,buffer) MPI_Datatype datatype, int root, LabSend(buf,dest) MPI_Send(void *buf, int count, LabSend(data,dest,tag) MPI_Datatype datatype, int dest, data=LabReceive(source) MPI_Recv(void *buf,int count, data=LabReceive(’any’,tag) MPI_Datatype datatype,int source, data=LabReceive(source,tag) int tag, MPI_COMM_WORLD, [data,source,tag]=LabResceive(.) MPI_Status *status) is_data_available=... MPI_Probe(int source,int tag, LabProbe(source,tag) MPI_COMM_WORLD,MPI_Status *status) received=labSendReceive(... MPI_Sendrecv(void *sendbuf, labTo,labFrom,data) int sendcount, MPI_Datatype sendtype, received = labSendReceive(... int dest, int sendtag, void *recvbuf, labTo,labFrom,data,tag) int recvcount, MPI_Datatype recvtype, Из табл. 1.1 можно заметить, что в отличие от MPI в приложении DCT системы MATLAB пока не предусмотрена возможность работы с различными коммуникаторами, разделяющими процессы на группы. Здесь имеется только один коммуникатор, который включает все процессы. Соответствующий коммуникатор в MPI называется MPI_COMM_WORLD.

Механизм тегирования – использование меток (тегов) для разбиения сообщений по типам – в MATLAB не является обязательным, однако он представлен целочисленными тегами в диапазоне от нуля до 32767. Если теги явно не указывать, то MATLAB самостоятельно присваивает всем рабочим процессам тег, равный нулю.

Функция mpigateway системы MATLAB дает возможность входа во внутреннюю функциональность MPI. Функция mpiLibConf позволяет локализовать реализацию MPI, которая будет использоваться. Функция mpiSettings(option, varargin) может использоваться для смены стандартных опций MATLAB для работы с коммуникациями MPI.

Отметим, что часто используемая в параллельном программировании операция Sendrecv была реализована только в MDCE версии 3.0. Также реализованы такие функции, как MPI_Reduce с операцией MPI_SUM и подобные ей.

За более подробной информацией можно обратиться к стандартной документации по Distributed Computing Toolbox [1].

1.8 Использование сторонних планировщиков MATLAB Distributing Computing Engine поддерживает следующие планировщики: встроенный от MathWorks (JobManager), Windows CCS - планировщик от Microsoft для Windows Server 2003 Cluster Edition, LSF, mpiexec либо планировщик общего вида. Под планировщиком общего вида может выступать, например, планировщик системы PBS [14] или Sun Grid Engine [15]. Среди преимуществ встроенного планировщика от MathWorks можно выделить:

• поддержку гетерогенности, • поддержку контрольных точек, • обеспечение автоматической доставки данных в процессы (свойство FileDependencies).

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

Вообще говоря, перечисленные процедуры (запуск службы mdce, запуск процессов jobmanager и worker) должен выполнять администратор кластера, а не тот пользователь-программист, который будет заниматься программированием самих вычислительных задач. Для программиста не имеет значения архитектура кластера, ему не нужно знать имена рабочих процессов, имена узлов, на которых они запущены и т.д.

Единственный объект, который нужен программисту-исследователю – это ссылка на объект планировщика – JobManager. Этого будет достаточно, чтобы создать объекты задач, подзадач, чтобы отправить задачу на счет.

Тем не менее материал, изложенный в данном разделе, сам по себе является достаточно интересным, чтобы его рассматривать отдельно.

Следует заметить также, что при использовании распределенных вычислений необходимо обеспечить идентичное наличие на каждом исполнительном узле как всех необходимых функций из разных toolbox, так и пользовательских файлов. Что касается toolbox, то лицензия на MDCE автоматически предполагает возможность установки всех возможных toolbox на узел.

Сделать это для пользовательских файлов можно тремя путями:

• либо пользователь самостоятельно рассылает все необходимые файлы в соответствующие места, • либо пользователь располагает файлы в общей для сервера и узлов папке; в этом случае для Windows имеет смысл сделать общую папку, а для Linux использовать файловую систему NFS, • либо пользователь использует механизм встроенного планировщика JobManager и MATLAB самостоятельно занимается распределением файлов (при отправке задания, MATLAB автоматически упаковывает рабочие файлы в zip-архив и распаковывает его в соответствующих сессиях рабочих процессов).

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

• Пользователь имеет права администратора на том компьютере, на котором происходит запуск соответствующего процесса.

• Брандмауэр операционной системы должен быть либо отключен, либо следует произвести соответствующую ему настройку.

• В среде Windows желательно минимизировать количество антивирусных программ, а лучше всего их вообще отключить.

II. Программирование После изучения настоящего раздела начинающий пользователь MATLAB сможет провести требуемые процедуры для создания распределенных и параллельных задач и понять, как описываются сами объекты – задачи.

Минимальное оборудование, необходимое для выполнения примеров из данного раздела – один компьютер с установленным на нем MATLAB Distributing Computing Toolbox, процесс jobmanager и один рабочий процесс. Этого будет достаточно, чтобы освоить технологию создания задачи, ее отправки на обсчет и возврата результатов.

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

В MATLAB в настоящее время существует два подхода для решения параллельных задач. Первый подход основан непосредственно на процедуре отправки задания jobmanager, в инструкциях (m - файле) которой описана последовательность команд, которая будут выполняться рабочими процессами. В этом m-файле помимо основных команд MATLAB могут быть использованы функции MPI для коммуникаций между рабочими процессами.

Второй подход для решения параллельных задач основан на режиме pmode. С помощью этого режима непосредственно из командного окна MATLAB становится возможным обращение к процессам workers, просмотр их локальных переменных, обмен данными между ними. В режиме pmode команды, вводимые в рабочем окне MATLAB, будут исполняться всеми рабочими процессами, ассоциированными с соответствующим jobmanager.

Режим pmode, по мнению авторов, следует использовать исключительно с двумя целями: как удобный пользовательский режим, предназначенный для первоначального знакомства с элементами параллельного программирования, для понимания многопроцессорной архитектуры и парадигмы параллельного программирования и как средство отладки параллельных программ.

В предыдущем разделе уже был рассмотрен один из способов получения ссылки на планировщика. Второй способ поиска ссылки на планировщика заключается в использовании свойства ’configuration’, которое есть описание того планировщика, ссылку на который необходимо найти.

Под описанием, имеются в виду следующие свойства:

минимальное и максимальное количество рабочих процессов, FileDependencies, LookupURL и т.д. Пользователь может создать (описать) несколько свойств непосредственно в m-файле matlabroot/toolbox/distcomp/user/distcompUserConfig.m, дав им соответствующие имена, и в дальнейшем осуществлять поиск ссылки на планировщик, используя только свойство ’configuration’, а не имя самого планировщика и узла, на котором он запущен.

Команда jm = findResource(’scheduler’,’configuration’, ’myconfig’) осуществляет поиск ссылки на системный процесс jobmanager, обладающего свойством ’myconfig’.

Пользователь может включить режим pmode, введя команду:

>> pmode start conf numlabs где conf - имя конкретной конфигурации планировщика, numlabs - количество рабочих процессов, которые должны быть запущены. Другие варианты включения режима pmode описаны в [1]. Итак, команда pmode start myconfig 2 включает режим pmode, ассоциирует его с планировщиком, описанным в свойстве myconfig, при этом используя только два рабочих процесса.

Следует отметить, что при исполнении этой команды ссылка на планировщик в рабочем пространстве MATLAB не появится.

После включения режима pmode обычный указатель на текущую строку MATLAB поменяет представление с привычного ">"на "P>>". Следующие две функции MPI будут рассмотрены в первую очередь:

labindex – уникальный идентификатор, порядковый номер рабочего процесса. labindex принимает значение от 1 до n, где n – общее число рабочих процессов, доступных в данной сессии (ассоциированных с данным планировщиком), numlabs – общее число рабочих процессов, ассоциированных с данным планировщиком.

При вводе в командной строке команды labindex будет получен следующий ответ системы:

P>> labindex В данном случае в рабочем окне отображен ответ системы так, как если бы команда labindex была введена (выполнена) непосредственно в каждом из процессов workers.

Введя команду numlabs, очевидно, мы получим следующий ответ системы:

P>> numlabs Иными словами, находясь в любой из двух сессий рабочих процессов, можно определить, сколько всего процессов доступно в данный момент.

Создадим теперь переменную a, которая будет определена в сессиях (рабочих пространствах) всех процессов. Пусть переменная a будет принимать значение ID процесса, т.е. labindex.

P>> a=labindex P>> whos Командой whos выводится список всех переменных, которые определены в текущей сессии. В обоих процессах в данном примере определена только переменная a.

Следует заметить, что в сессии рабочих процессов доступны только функции тех toolbox, которые установлены на клиентской части (на том компьютере, где производится запуск параллельной задачи ). Вообще говоря, в силу того, что любой рабочий процесс есть не что иное, как системный процесс MATLAB, ассоциированный с программой MATLAB, установленной на соответствующем узле, в его сессии доступны все функции всех пакетов расширений, в том числе и SIMULINK.

2.2 Вычисление определенного интеграла Информации, представленной в предыдущем разделе, достаточно для решения первой "классической" параллельной задачи - вычисления значения числа. Как известно, значение числа равно значению определенного интеграла:

Чтобы решить эту задачу, проинтегрировать функцию Рис. 2.1. Схематичный рисунок к задаче вычисления числа F (x) = 1+x2, следует воспользоваться тем свойством, что каждый из доступных рабочих процессов может интегрировать функцию F (x) на своем участке интервала (схематически это проиллюстрировано на рис. 2.1).

Очевидно, что первым шагом решения задачи должно быть распространение функции F (x) среди всех рабочих процессов.

Вводим команду P>> F = @(x) 4./(1 + x.^2) 1: @(x) 4./(1 + x.^2) 2: @(x) 4./(1 + x.^2) P>> Теперь функция F (x) становится переменной, определенной в каждом рабочем процессе.

Все, что осталось сделать перед тем, как запустить задачу вычисления интеграла, – необходимо определить границы интегрирования для каждого рабочего процесса. Пусть a и b – переменные, которые определяют границы интегрирования для каждого рабочего процесса. Вводим команду P>> a = (labindex - 1)/numlabs;

P>> b = labindex/numlabs;

Эти переменные будут определены в сессиях каждого рабочего процесса, причем значения их различны.

P>> [a, b] 1: ans = 2: ans = Теперь, воспользовавшись функцией quadl(F,a,b), в параметры которой передаются подынтегральная функция и граница интегрирования, вычислим значение определенного интеграла в каждом из рабочих процессов.

P>> myIntegral = quadl(F, a, b) 1: myIntegral = 2: myIntegral = Теперь, после того как в каждой сессии определена переменная myIntegral, необходимо воспользоваться функцией gplus и произвести глобальное суммирование переменной myIntegral по всем процессам.

gplus – gplus(x) возвращает сумму значений переменной x по всем процессам. Результат тиражируется на всех сессиях процессов.

P>> piApprox = gplus(myIntegral) 1: piApprox = 2: piApprox = Если требуется, переменную, хранящуюся в сессии рабочего процесса, можно экспортировать в текущую (пользовательскую) сессию MATLAB. Предварительно нужно переключиться из режима pmode в обычный режим MATLAB с помощью команды suspend.

Команда lab2client позволяет копировать определенную переменную из любого процесса в пользовательскую сессию.

P>> pmode suspend >> pmode lab2client piApprox Итогом рассмотренной задачи является переменная piApprox, определенная в локальной сессии MATLAB. Для сравнения полученного значения со значением константы pi, которое хранится в системе MATLAB, рассмотрим разность piApprox - pi >> piApprox - pi 2.4866e- Если пользователю после работы в локальной сессии необходимо снова переключится в режим pmode, то это можно сделать, введя команду >> pmode resume P>> Суммируя выше сказанное в данном разделе, можно сделать следующие выводы:

• Режим pmode работает без ссылки на jobmanager в клиентской сессии.

• Существует возможность копировать переменные из клиентской сессии в сессии рабочих процессов и наоборот.

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

• Рассматривать режим pmode следует исключительно как средство для первоначального знакомства с элементами параллельного программирования в MATLAB и как отладочное средство.

В данном разделе будут описана технология создания параллельных программ с использованием m-файлов. Данная технология существенно отличается от режима pmode и позволяет пользователю при написании параллельных программ использовать редактор кода MATLAB. При этом подходе становится возможным программировать в клиентской сессии MATLAB параллельные программы практически любой сложности.

Перейдем непосредственно к описанию технологии создания параллельных задач. Одним из методов объекта планировщика (jm, jobmanager) является метод createParallelJob 1. Методом createParallelJob можно воспользоваться двумя способами: передав функции createParallelJob в качестве методы любого объекта можно узнать, введя команду methods(var), где var - переменная MATLAB.

входного параметра переменную-ссылку на объект планировщик pjob=createParallelJob(jm), либо через операцию "точка", обратившись напрямую к этому методу pjob= jm.createParallelJob. Результаты выполнения обеих команд должены быть идентичными.

>> pjob = jm.createParallelJob pjob = distcomp.paralleljob Переменная pjob является ссылочной переменной и ссылается на определенный участок оперативной памяти, выделенной для системного процесса планировщика. По умолчанию в объекте pjob заполнены следующие свойства:

>> get(pjob) MaximumNumberOfWorkers: Inf MinimumNumberOfWorkers: FileDependencies: {0x1 cell} PathDependencies: {0x1 cell} Смысл большинства свойств понятен по названию, с детальным описанием всех свойств можно ознакомится в [1].

Обратимся для начала к свойствам MaximumNumberOfWorkers и MinimumNumberOfWorkers.

В силу того, что задача будет отправляться на вычисление через планировщик jm (у которого в свойствах указано, сколькими занятыми и сколькими свободными рабочими процессами он обладает в данный момент времени), уже в самой задаче (в объекте pjob) должно быть указано, сколько свободных процессов требуется планировщику для начала вычислений, описанных в задаче pjob.

Значение свойств MaximumNumberOfWorkers и MinimumNumberOf Workers есть положительное целое число. Установим его равным двум, т.е. вычисления в задаче pjob начнутся, как только у планировщика jm появятся два свободных рабочих процесса.

Записывать команды, которые устанавливают свойствам определенные значения, можно как в командном окне, так и непосредственно в m-файле, что является более удобным на наш взгляд. Напомним что с версии MATLAB R2006 редактор кода позволяет выделять ячейки с помощью двух символов %. Так, например, участок кода set(pjob,’MinimumNumberOfWorkers’,2) set(pjob,’MaximumNumberOfWorkers’,2) устанавливает нужные нам значения свойств. Для проверки значения свойств можно воспользоваться командой get и вернуть значение свойства, например выполнив get(pjob,’MinimumNumberOfWorkers’).

2.3.1 Объект ”параллельное задание” В данном разделе описывается, пожалуй, один из самых важных объектов в рассматриваемой нами архитектуре, а именно непосредственно параллельное задание. По сути параллельное задание есть не что иное, как m-файл, который будет исполняться одновременно всеми рабочими процессами. Передать m-файл напрямую рабочему процессу нельзя, для этой цели необходимо использовать объект - Task (параллельное задание). Объект Task - является дочерним объектом по отношению к объекту distcomp.paralleljob и может быть создан посредством метода createTask.

Создадим m-файл и опишем в нем функцию (процедуру), которая будет выполняться каждым процессом (по сути эта процедура состоит из команд, которые были введены при описании режима pmode).

При написании этого m-файла следует постоянно помнить о том, что данный m-файл будет выполняться одновременно несколькими процессами. Переменные, которые будут определены в этом m-файле, будут локальными для каждой сессии процесса, к ним не будет открыт доступ от других процессов без использования соответствующих функций MATLAB для передачи сообщений между процессами.

Следующий участок кода представляет собой m-файл (функцию), который состоит из уже описанных команд:

function piapprox=par_pi(F,a,b) % инициализация a=(labindex - 1)/numlabs; b=labindex/numlabs;

% вычисление интеграла myInt=quadl(F,a,b);

% суммирование по всем workers piapprox=gplus(myInt);

После того как написали m-файл, следует указать в свойстве FileDependencies объекта pjob имя файла, который будет выполняться рабочими процессами (вместо имени файла может быть указана и директория, если задача требует создания дополнительных m-файлов) set(pjob, ’FileDependencies’, {’par_pi.m’});

После этого создаем переменную, которая будет является подзадачей. Во входных параметрах требуется указать переменную - параллельную задачу, количество выходных аргументов (в рассматриваемом примере - 1) и входные аргументы {F,a,b} (переменные F,a,b должны уже быть определены в локальной сессии MATLAB):

obj = createTask(pjob,’par_pi’,1,{F,a,b}) Теперь все, что осталось сделать - это выполнить следующий участок кода:

submit(pjob); waitForState(pjob);

Здесь команда submit отправляет задачу pjob планировщику jobmanager. Задача начинает решаться рабочими процессами не сразу, а только тогда, когда свободными окажутся столько процессов, сколько указано в свойстве MinimumNumberOfWorkers.

До этого момента задача находится в очереди. Как только количества свободных процессов оказалось достаточно, задача pjob начинает "считаться"(иными словами, рабочие процессы начинают выполнять один и тот же m-файл). В этот момент свойство State объекта pjob принимает значение running.

Свойство поменяет свое значение на finished только после того, как задача, описанная в m-файле, выполнится. Как только задача поменяла статус на finished, можно воспользоваться методом getAllOutputArguments и получить данные в локальную сессию MATLAB.

results = getAllOutputArguments(pjob) Следует отдельно отметить, что при отправке задачи планировщик предварительно упаковывает в zip-архив файлы, указанные в свойстве FileDependencies, после чего отправляет его всем процессам, где распаковывает его в локальных сессиях, что гарантирует наличие исполняемого m-файла у всех рабочих процессов.

В целом m-файл, в котором описаны основные шаги для создания m-файла (программы на языке MATLAB, которая непосредственно отправляет задачу планировщику и возвращает результаты расчета) выглядит следующим образом.

%% Параллельная программа для вычисления числа pi jm = findResource(’scheduler’,’type’,’jobmanager’,...

’Name’,’MyJobManager’) %% Создание параллельной задачи pjob = jm.createParallelJob;

P=2;

set(pjob,’MinimumNumberOfWorkers’,P);

set(pjob,’MaximumNumberOfWorkers’,P);

set(pjob, ’FileDependencies’, {’par_pi.m’});

F= @(x) 4./(1 + x.^2);

a=0;

b=1;

%% Создание параллельной подзадачи obj = createTask(pjob,’par_pi’,1,{F,a,b});

submit(pjob); waitForState(pjob);

results = getAllOutputArguments(pjob) В дальнейшем мы не будем столь подробно останавливаться на этой последовательности, а обратимся непосредственно к самому важному, на наш взгляд, моменту, а именно к написанию mфайлов с использованием функций передачи сообщений (тех mфайлов (функций на языке MATLAB), которые исполняются непосредственно рабочими процессами).

В данном разделе будет описана широко известная задача параллельного матричного умножения. Фактически будет описано, каким образом, используя несколько процессов, можно Рис. 2.2. Кольцо из numlabs процессоров решить следующую задачу:

Данная задача подробно описана в разд. 6.1. [16]. Здесь мы остановимся на основных моментах, связанных с ней. По сути, технология, предлагаемая при использовании Distributing Computing toolbox и MATLAB Distributing Computing Engine, есть система с распределенной памятью (так как каждый рабочий процесс обладает своим адресным пространством). Схематично эта система (топология) представлена на рис. 2.2.

На начальной стадии перед началом вычислений в каждой сессии (в каждом рабочем пространстве процесса) доступен соответствующий блок матрицы A, вектора x и вектора y.

Вычисления начинаются с того, что каждый рабочий процесс отправляет доступный ему блок вектора x соседнему (по номеру ID) процессу. Следующим шагом является получение с помощью функции labReceive блока вектора x от соседнего рабочего процесса, который к этому моменту уже отправил его. После получения блока x происходит непосредственное вычисление вектора yk = yk + Ap, x, где – номер доступного блока вектора x (на каждом шаге, у каждого процесса рабочего процесса определена своя переменная ). Схематично на примере трех процессов перемещение подвектора xk представлено в таб. 2.1.

Таблица 2.1. Перемещение подвектора xk для трех процессов Ниже приведен листинг программы (m-файла), выполняемой каждым рабочим процессом. В этой программе переменная tau определяет номер текущего доступного блока вектора x.

Функция параллельного матричного умножения function z=par_multiplication(A,x,y) % A - матрица % x,y,z ветор-столбцы % z=Ax+y;

q=mod(n,numlabs);% остаток от деления if q== r=n/numlabs;

x=x((labindex-1)*r+1:labindex*r);

y=y((labindex-1)*r+1:labindex*r);

A=A((labindex-1)*r+1:labindex*r,:);

else % делится с остатком r=(n-q)/numlabs;

if labindex== x=x((labindex-1)*r+1+q:labindex*r+q);

y=y((labindex-1)*r+1+q:labindex*r+q);

A=A((labindex-1)*r+1+q:labindex*r+q,:);

end for t=1:numlabs if labindex==numlabs labSend(x,1); % последний сразу отправил первому labSend(x,labindex+1);

% остальные отправляют своим соседям (справа) if labindex== x=labReceive(numlabs);

% первый получил от последнего x=labReceive(labindex-1);

% остальные получают от своих соседей (слева) tau=labindex-t;

K0 > 0, a (0, 1), b > 1.

Однородная производственная функция с постоянной эластичностью замещения успешно применялась при исследовании экономики СССР и США (см., например, [25]), ее можно применить для исследования экономики современной России. Для учета научно-технического прогресса в [25] применялась модификация производственной функции (3.1), в которую добавлялся сомножитель в виде показательной функции времени exp(t).

Здесь мы рассматриваем короткий период времени и пока абстрагируемся от учета научно-технического прогресса.

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

Но, во-первых, такой подход полностью оправдан только при отсутствии модели, в которой используется данная производственная функция, поскольку все параметры модели должны быть согласованы. Найти параметры производственной функции, внутренне согласованные с поведением других переменных модели, весьма сложно в силу чрезвычайно большого Эластичность замещения труда и капитала для CES-функции (3.1) выражается через параметр b: = 1/(1 + b).

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

Во-вторых, что более существенно, при рассмотрении современной российской экономики нужно иметь ввиду, что значения второго фактора производственной функции, капитала K(t), представляемые статистическими органами, вызывают большое сомнение (см., например, [26]).

Статистические данные по капиталу (основным фондам страны) практически не меняются от года к году, а их стоимость чрезвычайно завышена. Фактически эти данные относятся к основным фондам, созданным во времена СССР, а они в настоящее время в большой степени представляют собой "бесплатные"ресурсы, подобные некоторым природным ресурсам, которые еще можно использовать без оплаты (воздух, а в некоторых случаях вода и земля). На выпуск оказывает влияние только капитал, вовлеченный в процесс производства, имеющий объективную стоимость, некий "эффективный" капитал, который при идентификации модели будет выражен в постоянных ценах 2000 г. и который мы и пытаемся здесь оценить.

Такое представление об основных фондах (совокупном физическом капитале страны) соответствует общепринятому в экономической теории определению физического капитала, как величины, равной запасу произведенных товаров, участвующих в производстве товаров и услуг (см., например, [27, c.322]).

Здесь мы измеряем труд среднегодовым числом занятых в народном хозяйстве, что отличает нашу работу от [26], где труд измерялся общим числом отработанных часов для учета его сезонных колебаний от квартала к кварталу. На основе статистических данных (см. ниже табл. 3.1) будем полагать, что труд L(t) растет с постоянным темпом > 0.

Будем считать, что капитал (эффективная стоимость производственных фондов) K(t) меняется в силу обычно применяемого в макроэкономических моделях уравнения (см., например, [28]) где µ темп4 выбытия капитала, а J(t) скорость прироста нового капитала (инвестиции в основной капитал).

Пусть в каждый момент времени t выполняется основной макроэкономический продуктовый баланс в текущих ценах:

сумма выпуска ВВП pY (t)Y (t) и импорта pI (t)I(t) равна сумме конечного потребления населения, правительства и некоммерческих предприятий с добавлением чистого накопления богатств и прироста материальных запасов, pC C(t), инвестиций в основной капитал pJ (t)J(t) и экспорта pE (t)E(t).

где через pY (t), pI (t), pC (t), pJ (t), pE (t) обозначены дефлятор ВВП и индексы цен на импорт, конечное потребление, инвестиции и экспорт, pY (2000) = pI (2000) = pC (2000) = pJ (2000) = pE (2000) = 1. Поскольку нас интересуют значения величин Обычно этот темп отождествляют с темпом амортизации капитала.

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

выпуска, инвестиций, экспорта и импорта, выраженные в постоянных ценах (ценах 2000 г.), то удобно перейти к продуктовому балансу, выраженному в индексах относительных цен Y (t) + I (t)I(t) = Q(t) + J (t)J(t) + E (t)E(t), где индексы относительных цен импорта, инвестиций и экспорта заданы отношениями а выражение, балансирующее равенство (3.5), дает величину Для решения системы уравнений (3.1) - (3.7) нужно определить объемы инвестиций J(t), экспорта E(t) и импорта I(t) в постоянных ценах 2000 г. Предполагаем здесь, что эти объемы определяются постоянными параметрами. Объем инвестиций в постоянных ценах J(t) определяется долей текущей стоимости инвестиций в сумме текущих стоимостей выпуска и импорта:

Объем экспорта в постоянных ценах E(t) определяется долей экспорта в выпуске (по их текущим стоимостям):

Объем импорта в постоянных ценах I(t) определяется отношением импорта к разнице ВВП и экспорта (по их текущим стоимостям):

Из (3.8)-(3.10) видно, что объемы экспорта E(t), импорта I(t) и инвестиций J(t) в постоянных ценах определяются параметрами,,, относительными ценами E (t), I (t), J (t) и выпуском (валовым внутренним продуктом) Y (t).

Итак, для идентификации модели надо задать изменение внешних интенсивных параметров модели: трех относительных цен E (t), I (t), J (t), - а также определить семь постоянных параметров a, b,, µ,,, и три начальных значения Y 0, K0, и L0 таким образом, чтобы расчетные временные ряды макропоказателей (переменных модели) были близки к статистическим временным рядам соответствующих макропоказателей экономики России.

Используемые статистические данные представлены в табл. 3.1. Источник данных: Федеральная служба государственной статистики РФ http://www.gks.ru. Составляющие ВВП приведены в постоянных ценах 2000 г., в млрд. руб. Индексы относительных цен получены расчетом на основе данных по изменению составляющих ВВП в текущих и постоянных ценах и имплицитного дефлятора ВВП. Данные для величины Q(t), определенной (3.7), рассчитаны из баланса (3.5).

На основе имеющихся в табл. 3.1 статистических данных можно непосредственно определить часть параметров модели, во всяком случае некоторые границы их изменения. По ходу такой идентификации параметров мы уточним и само описание модели.

Определим экспоненциальную линию тренда для числа занятых L(t) согласно статистическим данным по занятости, представленным в табл. 3.1 (как среднеквадратическое отклонение расчетных и статистических значений). Ниже приведен листинг участка кода на языке MATLAB, в котором реализована подгонка методом наименьших квадратов статистических данных для числа занятых в экономике России в 2000-2006 гг. L(t) к экспоненциальной функции.

Подгонка L(t) экспоненциальной функцией %% Подгонка экспоненциальной моделью L=[65.273 65.124 66.358 67.247 67.244 68.719 69.6];

t=[2000 2001 2002 2003 2004 2005 2006];

tm2000=t-2000;

% exp1 - имя модели из библиотеки моделей cflibhelp [Lfit,expgodness] = fit(tm2000’,L’,’exp1’);

% полученная подгонка используется в идентификации Lestim = Lfit(tm2000);

Lfit % характеристики подгонки expgodness % рисунок "Сравнение расчета со статистикой по труду L" figure plot(t,L,’ob’,t,Lestim,’or’,’LineWidth’,2); hold on h = legend(’L stat’,’L estim’,2);title(’L’) plot(t,L,’-b’,t,Lestim,’-r’,’LineWidth’,2);

В результате экспоненциальной подгонки в командном окне получим данные, представленные ниже, что дает следующее за ними соотношение для L(t):

Результат подгонки параметров функции L(t) Lfit = General model Exp1:

Lfit(x) = a*exp(b*x) Coefficients (with 95% confidence bounds):

a = 64.84 (64.14, 65.55) b= 0.01122 (0.008283, 0.01419) expgodness = rsquare: 0. adjrsquare: 0. Для уравнения (3.2) соотношение (3.14) дает значение темпа роста занятого населения параметра = 0.01124 и значение начальной величины числа занятых L 0 = L(2000) = 64.84, при этом последнее несколько меньше статистического значения из табл. 3.1. Для дискретного с единичным шагом варианта уравнения (3.2): Lt+1 = Lt (1 + ), - имеем = exp0.01124 1 = 0.0113, L0 = 64.84. В дальнейшем при совокупной оценке параметров модели в дискретном времени можно, например, искать два данных параметра, исходя из доверительных границ в 95%: L0 [64.14, 65.55], = [0.00832, 0.01420].

Можно также найти линии тренда для индексов относительных цен E (t), I (t) и J (t) согласно данным табл. 3.1 (как средL, correlation=0. Рис. 3.1. График экспоненциальной подгонки труда L(t) на основе неквадратическое отклонение расчетных и статистических значений), например по формулам В этих уравнениях условия нормировки E (2000) = I (2000) = J (2000) = 1 соблюдаются автоматически. Фрагмент кода для вычисления параметров функции E (t) представлен ниже.

Подгонка индекса относительных цен E (t) %% Подгонка piE моделью x(1)+(1-x(1))*exp(-x(2)*(t-2000)) piE=[1 0.84442 0.76610 0.72863 0.68475 0.69651 0.67010];

t=[2000 2001 2002 2003 2004 2005 2006];

tm2000=t-2000;

ydata=piE;lb=[0 -inf]; ub=[inf inf];

[xE resnorm]= lsqcurvefit(@funE,[0 0],tm2000,ydata,lb,ub) piEestim=funE(xE,tm2000);

figure plot(t,piE,’ob’,t,piEestim,’r-’,’Marker’,’s’) title(’piE’) function piX=funE(x,xdata) piX=x(1) + (1-x(1))*exp(-x(2)*xdata);

%%% Результат, отображаемый в командном окне %%% resnorm = 4.2293e- Необходимый код для вычисления параметров функций I (t), J (t) легко написать по аналогии, что начинающему пользователю MATLAB рекомендуется сделать самостоятельно.

Ниже представлен код соответствующих этому коду вспомогательных функций на языке MATLAB.

Рис. 3.2. Подгонка относительной цены E по данным табл. 3. Вспомогательные функции для оценки индексов относительных цен на импорт I и на инвестиции J function piX=funI(x,xdata) piX=1-x(1) *xdata.*xdata.*exp(-x(2)*xdata);% x(1)> function piX=funJ(x,xdata) piX=x(1) + (1-x(1))*(1+xdata).*exp(-x(2)*xdata);% x(1)> Рис. 3.3. Подгонка относительной цены I по данным табл. 3. Результаты подгонки относительных цен представлены на рис. 3.2 - 3.4.

Среднее значение отношения объема инвестиций в основной капитал к сумме ВВП и импорта (3.8) с 2001 г. по 2006 г. остается практически постоянным: = 0.1346 ± 0.0026 (первая цифра среднее значение, вторая - среднеквадратичное отклонение), так Рис. 3.4. Подгонка относительной цены J по данным табл. 3. что [0.1320, 0.1372]. Значение этой величины за 2000 г., = 0.1286, чуть-чуть не укладывается в этот интервал, что можно объяснить тем, что 2000 г. еще следует относить к переходному периоду экономики к новой структуре после дефолта 1998 г.

Для периода с 2001 г. по 2006 г. доля экспорта в ВВП (3.9) = 0.3511 ± 0.0103, а для переходного 2000 г. = 0.4406. Если учитывать только период времени с 2001 г. по 2006 г., интервал для поиска наиболее близкого к статистике будет меньше:

[0.3407, 0.3614].

То же самое касается параметра отношения импорта к остатку от ВВП, после вычета из него объема экспорта (3.10) с 2001 г. по 2006 г. = 0.3532 ± 0.0264, а для переходного для установления импорта 2000 г. = 0.4296. Если учитывать только период времени с 2001 г. по 2006 г., интервал для поиска наиболее близкого к статистике меньше: [0.3268, 0.3796].

Для точного определения пятнадцати параметров модели a, b, K(0), µ, L(0),, aE, bE, aI, b(I), aJ, bJ,,, можно использовать косвенный подход, сравнивая рассчитанные по модели экономики России выходные временные ряды ее переменных с доступными статистическими временными рядами этих переменных для 2000 - 2006 гг. Временные ряды считаются похожими, если они близки как функции времени (другими словами, между значениями временных рядов существует сильная, возможно нелинейная, связь). Поскольку длины статистических временных рядов, которым мы доверяем, составляют здесь небольшую величину в шесть значений, мы не будем здесь использовать индекс похожести, введенный на основе вейвлет-коэффициентов для сравнения нелинейных временных рядов в [29-31]. Вместо него будем использовать коэффициент корреляции Пирсона D(X, Y ), который является мерой силы и направленности линейной связи между сравниваемыми временными рядами X, Y, и чем он ближе к единице, тем более схоже поведение этих рядов. При этом следует учитывать, что инфляционная составляющая может преувеличивать линейную связь рядов, поэтому при использовании коэффициента корреляции нужно сравнивать показатели в реальных величинах. Если длины сравниваемых временных рядов равны n, то имеем следующее выражение для коэффициента корреляции Пирсона:

где T – знак транспонирования.

Индекс Тейла E(X, Y ) измеряет несовпадение временных рядов Xt и Yt и чем ближе он к нулю, тем ближе сравниваемые ряды [32]. Для удобства проведения расчетов со сверткой критериев в дальнейшем вместо индекса Тейла будем использовать коэффициент близости U (X, Y ) = 1 E(X, Y ):

Чем выше он (чем ближе он к единице), тем более близки ряды.

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

Для однозначности выбора оптимального варианта можно использовать ту или иную свертку коэффициентов близости U (X, Y ) и корреляции D(X, Y ); например, если подгонка расчетных и статистических данных для всех макропоказателей имеет примерно равную важность, можно максимизировать среднегеометрическую величину всех коэффициентов.

В формальной записи где множество параметров задано на параллелепипеде а функционал представляет собой среднегеометрическое всех критериев близости и корреляции:

Здесь m - число макропоказателей; j - номер макропоказателя, При этом при переборе следует рассматривать только те варианты значений параметров, при которых коэффициенты близости и корреляции выше некоторых заданных положительных величин; например, в нашем случае можно использовать ограничения Dj > 0.7, Uj > 0.85, где j = 1,..., m.

В первом приближении при определении параметров модели можно для параметров,, взять их средние значения и использовать выражения (3.14), (3.18)-(3.20) для задания внешних величин L(t), E (t), I (t), J (t). Тогда вместо пятнадцати неопределенных параметров останется только четыре, a, b, K(0), µ, которые можно найти простым перебором.

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

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

3.2 Численная реализация задачи идентификации Сложные математические модели экономических систем, как правило, рассчитывают по численной схеме с дискретным временем. Поскольку наша задача показать основные проблемы и пути их решения в общем случае, мы и для расчета нашей простейшей модели будем применять дискретную схему. Более того, для простоты изложения мы шаг по времени примем равным единице. Тогда уравнение (3.3) в дискретном виде выглядят так:

Здесь дискретный момент времени t = 0 соответствует 2000 г., величина K0 равна искомому начальному запасу капитала, а параметр µ, если он положителен, определяет долю амортизируемого в течение года капитала. Заметим, что единичный шаг по времени удобен для сравнения рассчитываемых по модели показателей с их статистическими аналогами, представленными в табл. 3.1.

Для остальных уравнений дискретная и непрерывная записи совпадают (для отличия двух форм записи мы в дискретной записи момент времени указываем с помощью нижнего индекса, а численные значения начальных данных помечаем верхней звездочкой). Простые алгебраические преобразования уравнений (3.4) - (3.20) дают выражение всех макропоказателей простейшей модели в момент времени t экспорта E t, импорта It, инвестиций в основной капитал Jt, балансирующей величины конечного потребления домашних хозяйств и правительственных организаций Qt как величин, пропорциональных валовому внутреннему продукту (выпуску) Yt и деленных на соответствующие индексы относительных цен.

Для положительности величины Q, в которую включено конечное потребление, требуется выполнение условия продуктивности системы.

Относительные индексы цен в (3.27) - (3.30) определяются соотношениями (3.18) - (3.20).

Теперь для упрощения работы с моделью перейдем в выражениях для труда Lt, капитала Kt и выпуска Yt к относительным величинам lt, kt, yt соответственно:

Начальные значения всех этих величин равны единице: l 0 = k0 = y0 = 1. Тогда (3.1) и (3.32) дают Из (3.14), (3.26), (3.32) получим В (3.34) = 0.01124. В (3.35) введены обозначения где aJ = 0.811, bJ = 0.5276.

Наша основная задача подобрать такой временной ряд для капитала, который наилучшим образом приближает временные ряды для макропоказателей, рассчитанных по модели, к их статистическим аналогам, представленным в табл. 3.1, поэтому численную реализацию идентификации (нахождения внешних параметров) модели мы и начинаем с варианта с наименьшим числом параметров. Мы предполагаем, что представленные ниже параметры фиксированы: L0 = 64.84, Y0 = 7305.6, = 0.01124, = 0.3532, = 0.3511, = 0.1346, aJ = 0.811, bJ = 0.5276. Тогда, в соответствии с (3.36) = 0.1569.

Итак, имеем соотношения модели (3.33)-(3.35), которые при каждом заданном наборе параметров a, b, µ, с помощью выражений (3.32)-(3.35) и (3.27)-(3.30) дают искомые временные ряды макропоказателей Yt, Lt, Kt, It, Qt, Jt, Et. Для сравнения близости расчетных временных рядов указанных макропоказателей с их статистическими аналогами надо вычислить критерии корреляции (3.21) и близости (3.22) для выпуска, потребления, инвестиций, экспорта и импорта - Yt, Qt, Jt, Et, It - за период 2001-2006 гг.

и вычислить свертку этих критериев (3.25).

Возможный интервал изменения оцениваемых параметров:

a (0, 1), b (1, 2), µ (0.2, 0.1), (0, 3). Для поиска параметров с помощью параллельных вычислений надо взять сетку по каждому из интервалов, устроить перебор всех возможных сочетаний, распараллелить этот перебор на доступное число процессоров. На каждом из процессоров отбросить варианты, в которых коэффициенты корреляции и близости не превышают 0.4. Среди оставшихся вариантов выбрать вариант с наибольшей сверткой (3.25), совокупным критерием F (a), отправить его номер процессору-мастеру, вычислить самый большой критерий среди полученных наибольших у процессоров-рабочих, а затем для полученного оптимального варианта рассчитать все временные ряды макропоказателей модели по формулам (3.27)-(3.30), (3.32)-(3.35), вывести все графики, сравнивающие расчет со статистикой. Схема алгоритма приведена ниже.

1. В цикле по сеткам для интервалов задания параметров задаются их значения через равный интервал (эти циклы и надо распараллелить для расчета разных групп наборов параметров в параллельных процессах, идущих на разных процессорах многопроцессорной вычислительной системы или разных ядрах многоядерной системы).

2. Для каждого набора параметров рассчитываются временные ряды всех макропоказателей L, Y, I, Q, J, E.

3. Для каждого макропоказателя, участвующего в сравнении расчетных и статистических данных, рассчитываются критерии (3.21)-(3.22) близости U и корреляции D между расчетными и статистическими временными рядами на промежутке времени 2001-2006 гг., для которого, как мы предполагаем, справедливы выдвинутые нами гипотезы о постоянности долей (3.8)-(3.10).

4. Рассчитываем свертку критериев (3.25).

5. В каждом процессе определяем лучший номер набора параметров по критерию (3.23).

6. На процессе-мастере выбираем лучший номер из лучших, полученных от рабочих процессов. Рассчитываем и распечатываем для него временные ряды всех показателей и значения всех критериев.

Описанный алгоритм идентификации динамической системы реализован в виде основного файла сценария, в котором происходит инициализация переменных, первичная статистическая подгонка данных и кластерной части, набора m-файлов, исполняемых на кластере.

Основной m-файл, выполняемый на кластере, состоит из трех логических блоков:

1. Инициализация сетки для каждого рабочего процесса (см.

код на языке MATLAB ниже).

slice=(0.9-0.1)/numlabs; % участок для каждого a1=0.1+(labindex-1)*slice;% начало отрезка a2=0.1+labindex*slice; % и конец отрезка for a = linspace(a1,a2,slice) % верхний цикл index1=index1+1;

2. Непосредственно вложенный цикл, в котором происходит вычисление значения целевой функции F (см. ниже код на языке MATLAB).

for a = linspace(a1,a2,slice) index1=index1+1;

for b = linspace(-0.9,-0.6,n) index2=index2+1;

index3=0;

for alphaa=linspace(0.04,0.9,N) K0=Y0/alphaa;

index3=index3+1;

index4=0;

for muu= linspace(-0.2,0.05,N) index4=index4+1;

[Yestim Kestim]=...

forecastYK(Y0, K0, Y0, K0, L0, a,...

betaa, Lestim, piJestim);

Jestim = betaa * Yestim./ piJestim’;

Eestim = deltaa * Yestim./ piEestim’;

Iestim = rho * (1 - deltaa) *...

Qestim = ((1-sigmaa) * (1 + rho*...

Ycor=pearson(Y’,Yestim);

Yomt=omth(Y’,Yestim);

Jcor=pearson(J’,Jestim);

Jomt=omth(J’,Jestim);

Yres=[Ycor Yomt Jcor Jomt];

if sum(Yres= F(index1,index2,index3,index4)=-1;

F(index1,index2,index3,index4)=...

end Внутри цикла происходит вызов вспомогательных функций (функций для расчета коэффициента близости по индексу Тейла, коэффициента корреляции Пирсона и прогноза выпуска и капитала по уравнениям модели), которые доступны каждому рабочему процессу, в силу того что планировщик автоматически доставляет и распаковывает их в каждой сессии рабочего процесса. Код этих вспомогательных функций на языке MATLAB представлен ниже.

Функция близости рядов по индексу Тейла function omth=omth(x,y) % omth = One minus Theil’s index numerator=sum((x-y).^2);

denominator=sum(x.^2)+sum(y.^2);

omth=1-sqrt(numerator/denominator);

Функция корреляции Пирсона function pearson=pearson(x,y) n=numel(x);

numerator=n*(sum(x.*y))- sum(x)*sum(y);

denominator= sqrt((n*sum(x.^2)-sum(x)^2)*...

pearson=numerator/denominator;

Функция прогноза выпуска и капитала function [Y K]=forecastYK(Yt, Kt, Y0, K0, L0,...

if length(L)==length(PiJ) N=length(L);

warning(’In function forecastYK l=L/L0;

y=zeros(N,1); Y=zeros(N,1);

k=zeros(N,1); K=zeros(N,1);

y(1)=Yt/Y0; k(1)=Kt/K0; % Yt, Kt for forecast for i=2:N k(i)=(1-muu)*k(i-1)+alphaa*betaa*y(i-1)/PiJ(i-1);

% to avoid 1/0 error % Cobb-Douglas production function y(i) = (l(i)^a) * (k(i)^(1-a));

% CES production function geterogenity c y(i)=(a*l(i)^(-b)+(1-a)*k(i)^(-b))^(-c/b);

Y=y*Y0; K=k*K0;

3. Финализация – этот шаг реализован в виде небольшого участка кода, в котором все рабочие процессы отправляют свои расчетные данные, максимальное значение целевой функции F и значения параметров, на которых это значение достигнуто, процессу-мастеру. На последнем происходит выбор максимального из присланных значений, что и является конечным результатом, т.е. выходным значением функции (см. код на языке MATLAB ниже).

(мастер выбирает наибольшее значение) OptimMatrix=...

[FFoptim aoptim boptim alphaaoptim K0 muuoptim];

if labindex ~=1 % отправляем первому worker (мастеру) for otherLab=2:numlabs labSend(OptimMatrix,1);

OptimParams=0;

% переменная OptimParams будет отличаться else % мастер принимает данные for otherLab=2:numlabs OptimMatrixGet=labReceive(otherLab) OptimMatrix=[OptimMatrix;OptimMatrixGet];

%выбираем наибольшее значение по F [FFoptim I]=max(OptimMatrix(:,1));

OptimRow=OptimMatrix(I,:);% вектор содержащий %оптимальные параметры % Инициализируем структуру, которая будет % являться выходным значением % OptimParams отлична от нуля только у мастера OptimParams.FFoptim=OptimRow(1);

OptimParams.aoptim=OptimRow(2);

OptimParams.boptim=OptimRow(3);

OptimParams.alphaaoptim=OptimRow(4);

OptimParams.K0=OptimRow(5);

OptimParams.muuoptim=OptimRow(6);

end Некоторые вспомогательные функции, введенные ранее, используются в цикле перебора вариантов вместо встроенных функций MATLAB, чтобы ускорить расчет. Встроенные функции MATLAB (например, коэффициент корреляции Пирсона corr(x’,y’)), удобно использовать в последовательных частях программы, если нужно подсчитать дополнительные характеристики сравниваемых временных рядов. В частях программы, которые обрабатываются большое количество раз, подготовка быстро исполняемого кода дает существенный выигрыш во времени.

Из кода для основного вложенного цикла видно, что коэффициенты близости и похожести расчетных и статистических временных рядов используются в свертке критериев и сами по себе только для двух макропоказателей, выраженных в постоянных ценах 2000 г.: валового внутреннего продукта Y (t) и инвестиций в основной капитал J(t). Это снова сделано для ускорения расчета параллельной части программы. Здесь учтен тот факт, что близость расчетных и статистических рядов для остальных показателей обеспечивается автоматически в силу соотношений (3.27)-(3.30), если близки расчетные и статистические временные ряды для макропоказателя Y (t).

В коде на языке MATLAB (m-файлах) для инициализации и основного вложенного цикла приведен пример распараллеливания только самого верхнего цикла. Можно распараллелить и большее число циклов, например объединив эти циклы в один цикл и распараллелив последний по приведенной выше схеме.

Начинающий пользователь MATLAB может проделать такое упражнение на распараллеливание объединение нескольких циклов в один и разбиение общего цикла на порции для используемых параллельно процессов самостоятельно, если для расчета задачи идентификации параметров математической модели сложной системы у него есть возможность использовать многопроцессорную технику, например кластерный суперкомпьютер, с числом процессоров, исчисляющихся десятками или сотнями.

3.2.2 Результаты идентификации модели, их Результаты идентификации модели представлены в графическом виде на рис. 3.5 - 3.7, где дано сравнение рассчитанных по модели и статистических временных рядов для макропоказателей экономики России.

Рис. 3.5. Данные для капитала K(t), рассчитанные на основе модели (Kestim), и данные для выпуска Y(t), рассчитанные по модели (Y estim) и статистические (Y stat) В результате идентификации модели экономики России получена оптимальная в смысле свертки (3.25) критериев близости и Рис. 3.6. Исторические данные и данные, рассчитанные на основе модели, для импорта (Istat, Iestim) и экспорта (Estat, Eestim)) корреляции сравниваемых макропоказателей Y (t) и J(t) оценка ее параметров a = 0.84, b = 0.78, µ = 0.175, = 0.41. На основе этих параметров из (3.36) определено начальное значение для капитала K0 = Y0 / = 17819 млрд. руб 2000 г., по уравнению (3.26) восстановлен временной ряд для капитала K(t), а с помощью уравнений (3.1), (3.27)-(3.30) восстановлены временные ряды рассматриваемых в модели макропоказателей для ВВП Y (t), инвестиций J(t), потребления Q(t), импорта I(t) и экспорта E(t).

Представленный на рис. 3.5 капитал является эффективным Рис. 3.7. Исторические данные и данные, рассчитанные на основе модели, для инвестиций в основной капитал (Jstat, Jestim) и капиталом, реально используемым в производстве. Он был получен при идентификации рассмотренной здесь модели экономики. Эффективный капитал отличается от значения, формально рассчитываемого статистическими органами на основе многократных переоценок в условиях огромной инфляции стоимости существующих производственных фондов.

Как ранее уже упоминалось, отрицательное значение параметра µ означает, что капитал прирастает с большей скоростью, чем это обусловлено инвестициями, за счет вовлечения в процесс производства простаивавшего ранее капитала (производственных фондов). Но этот процесс вовлечения не может продолжаться долго, поскольку объем неиспользуемого в процессе производства капитала, доставшегося современной экономике от советских времен, ограничен.

Оценим время T с 2000 г., за которое будет вовлечен весь неиспользуемый до конца этого года капитал. Будем исходить из следующих правдоподобных предположений: (1) объем новых инвестиций совпадает с объемом капитала, выбывающего вследствие физического и морального износа; (2) за весь процесс вовлечения капитал может вырасти в четыре раза в сравнении с его начальным уровнем в 2000 г. Тогда общее время процесса вовлечения старого капитала истечет через T = (1/|µ|)ln(K T /K0 ) = ln(4)/0.175 8 лет после 2000 г., т.е. в 2008 г. исчерпается лимит вовлечения простаивавших производственных фондов.

Другая проблема, связанная с вовлечением старых производственных фондов, состоит в том, что они чрезвычайно изношены.

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

На рис. 3.6 - 3.7 представлены статистические и рассчитанные по модели временные ряды макропоказателей, характеризующих структуру использования реального валового внутреннего продукта в ценах 2000 г. (экспорта E(t), импорта I(t), инвестиций J(t) и потребления Q(t)). В 2000-2006 гг. все рассмотренные макропоказатели растут (см. рис. 3.5 - 3.7). Поведение этих показателей после 2006 г. можно оценить, задав сценарии будущего развития.

В сценарных расчетах с моделью прогноз изменения внешних параметров модели задается сценарием. Будем считать, что индексы относительных цен во всех рассматриваемых сценариях заданы функциями (3.18)-(3.20), так что мы имеем рис. 3.8 - 3.10.

Рис. 3.8. Заданный индекс относительных цен на инвестиции в 3.3.1 Базовый, он же пессимистический сценарий Обычно в качестве базового сценария развития экономики рассматривают такой сценарий, в котором на прогнозный период времени предполагают продолжение тенденций, выявленных за период времени, на котором оценивалась модель. Однако в нашем случае так сделать нельзя, поскольку на прогнозный период Рис. 3.9. Заданный индекс относительных цен на экспорт E (t) времени с 2007 г. по 2020 г. по мнению многих экономистов выполняются следующие условия.

1. Источник вовлечения производственных фондов из наследства, доставшегося от советских времен, вот-вот будет исчерпан.

2. Производственные фонды, используемые в производстве, сильно изношены доля производственных фондов, которые скоро выйдут из строя и которые поэтому надо срочно менять, в большинстве отраслей народного хозяйства превысила 50%.

Рис. 3.10. Заданный индекс относительных цен на импорт I (t) 3. Прирост числа занятых в экономике, наблюдавшийся в период оценки, в прогнозный период практически невозможен, так как начнет сказываться демографический кризис, в который Россия попала в начале 90-х гг. XX в.

Поэтому базовый вариант прогноза оказывается пессимистическим.

Зададим базовый сценарий прогноза с 2007 г. до 2020 г. (он же пессимистический) следующими условиями, накладываемыми на внешние параметры модели:

Рис. 3.11. Данные для капитала K(t) (базовый сценарий) 1. Считаем, что параметры a, b,,,, остаются постоянными, принимающими определенные ранее значения для 2000-2006 гг.: a = 0.84, b = 0.78, = 0.41, = 0.3511, = 0.1346, = 0.3532, так что = 0.1569.

2. Параметр µ резко меняет свое значение и экономический смысл с 2009 г.: µ = 0.175 < 0 до 2008 г., а начиная с 2009 г. он становится положительным µ = = J 0 /K0 = 0.0678 и означает темп выбытия мощностей вследствие износа. Предполагаем, что источник, из которого в производство вовлекались бесплатные производственные фонды, Рис. 3.12. Прогноз по базовому (пессимистическому) сценарию к началу 2009 г. иссякнет. Капитал меняется в силу уравнения (3.26) или в относительных величинах в силу уравнения 3. Считаем, что объем используемого в производстве труда до 2008 г. возрастает в силу соотношения (3.14), достигает максимального значения L8 = 70.94 млн. человек в 2008 г.