В данной строке имя файла с интерфейсной процедурой является вторым параметром mex. После выполнения данной команды в текущем каталоге появится файл mysum.dll, который и содержит требуемую МЕХ-функ- цию. Проверим ее работу, обратившись к ней из командной строки:

>> s=mysum(l,2) s =

3

Пример создания МЕХ-функции из простейшей процедуры Fortran, приведенный выше, демонстрирует основные принципы использования интерфейса приложений Matlab API. Пользователю Matlab, не имеющему достаточного опыта в области разработки интерфейса приложений, но желающему вызывать процедуры Fortran из среды Matlab, достаточно знать небольшой набор функций Matlab API. Требуются только функции, которые позволяют обмениваться основными типами данных. Интерфейсные процедуры, как правило, принципиально не отличаются для различных процедур на Fortran и состоят из трех частей:

•получение значений входных аргументов;

•вызов процедуры на Fortran;

•возвращение значений в среду Matlab.

=========================================================

Пример из документации Matlab

Создание МЕХ-файла, использующего программу на Фортране

Напишем простую программу на Фортране timestwo, которая принимает в качестве входного параметра действительное значение переменной, а возвращает это значение, умноженное на 2:

subroutine timestwo(y_output, x_input) real*8 x_input, y_output

y_output = 2.0 * x_input return

end

Создадим интерфейсную программу для вызова этой подпрограммы на Фортране. Откроем MATLAB Editor и введем следующую информацию

C=======================================================

C=======================================================

7

Добавимзаголовочныйфайл Фортранаfintrf.h, позволяющийпод-

ключить API MATLAB.

#include "fintrf.h"

Сохраним файл в рабочую директорию MATLAB, например, на сменный носитель f:\work, и назовем его timestwo.f.

Длядоступак подпрограмме наФортранеиспользуется интерфейсная процедура mexFunction, которая должна содержать следующий код

C Интерфейсная процедура

subroutine mexFunction(nlhs, plhs, nrhs, prhs)

C Declarations

C Statements

return end

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

implicit none

Явное объявление типа переменных необходимо для 64-битных массивов.

Интерфейсная процедура mexFunction имеет четыре аргумента. Два типа mwPointer (указатели на массивы типаmxArray) и два типа inte-

ger:

nrhs – число входных аргументов, с которыми будет происходить обращение к МЕХ-функции timestwo из среды Matlab;

prhs – массив указателей на входные аргументы timestwo;

nlhs – число выходныхаргументов,с которыми будетпроисходитьобращение к МЕХ-функции timestwo из среды Matlab;

plhs — массив указателей на выходные аргументы timestwo.

Обмен данными между процедурой timestwo и средой Matlab посредством интерфейсной процедуры mexFunction осуществляется только при помощи перечисленных аргументов. Для выделения данных из prhs и занесения результата в plhs предназначены специальные функции Matlab API.

Добавим следующий код после оператора Declarations.

C mexFunction arguments: mwPointer plhs(*), prhs(*) integer nlhs, nrhs

Объявление функций и локальных переменных

8

В mexFunction использованы функции MATLAB API: mxGetPr,mxCreateDoubleMatrix, mxGetM, mxGetN возвращающие значения типа integer. Функция mxCreateDoubleMatrix вы-

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

первый вещественный элемент в массиве. Функции mxGetM и mxGetN возвращают количество строк и столбцов в массиве соответственно. Переменная mxIsNumeric определяет является ли массив числовым.

C Function declarations: mwPointer mxGetPr

mwPointer mxCreateDoubleMatrix integer mxIsNumeric

mwPointer mxGetM, mxGetN

Локальные переменные x_ptr и y_ptr понадобятся для хранения указателей на массивы, содержащие входной и выходной аргументы. Переменные mrows, ncols определяют размерность массива для хранения данных. В нашем случае имеем массив 1х1.

Объявим локальные переменные для аргументов mexFunction

C Pointers to input/output mxArrays: mwPointer x_ptr, y_ptr

Объявим переменные, характеризующие используемый массив

C Array information: mwPointer mrows, ncols mwSize size

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

C Arguments for computational routine: real*8 x_input, y_output

Чтение входного массива

Указатель на массив, содержащий входные данные получаем с помощью функции mxGetPr

x_ptr = mxGetPr(prhs(1))

Для создания массива Фортрана x_input, используем функцию mxCopyPtrToReal8.

9

C Get the size of the input array. mrows = mxGetM(prhs(1))

ncols = mxGetN(prhs(1)) size = mrows*ncols

C Create Fortran array from the input argument. call mxCopyPtrToReal8(x_ptr,x_input,size)

Для создания выходной величины используем функцииmxCreateDoubleMatrix и mxGetPr

C Create matrix for the return argument.

plhs(1) = mxCreateDoubleMatrix(mrows,ncols,0) y_ptr = mxGetPr(plhs(1))

Выполнение вычислений

Передадим аргументы в подпрограммуtimestwo.

C Call the computational subroutine. call timestwo(y_output, x_input)

Копирование полученного результата в выходную величину

C Load the data into y_ptr,

C which is the output to MATLAB.

call mxCopyReal8ToPtr(y_output,y_ptr,size)

В окончательном виде файлtimestwo.f будет выглядеть так

C ==================================================

C ==================================================

#include "fintrf.h"

CGateway routine

subroutine mexFunction(nlhs, plhs, nrhs, prhs)

C mexFunction arguments: mwPointer plhs(*), prhs(*) integer nlhs, nrhs

C Function declarations: mwPointer mxGetPr

mwPointer mxCreateDoubleMatrix

10

integer mxIsNumeric mwPointer mxGetM, mxGetN

C Pointers to input/output mxArrays: mwPointer x_ptr, y_ptr

C Array information: mwPointer mrows, ncols mwSize size

C Arguments for computational routine: real*8 x_input, y_output

x_ptr = mxGetPr(prhs(1))

C Get the size of the input array. mrows = mxGetM(prhs(1))

ncols = mxGetN(prhs(1)) size = mrows*ncols

C Create Fortran array from the input argument. call mxCopyPtrToReal8(x_ptr,x_input,size)

C Create matrix for the return argument.

plhs(1) = mxCreateDoubleMatrix(mrows,ncols,0) y_ptr = mxGetPr(plhs(1))

C Call the computational subroutine. call timestwo(y_output, x_input)

C Load the data into y_ptr,

C which is the output to MATLAB.

call mxCopyReal8ToPtr(y_output,y_ptr,size)

return end

C---------------------------------------------------

CComputational routine

subroutine timestwo(y_output, x_input) real*8 x_input, y_output

y_output = 2.0 * x_input return

11

end

Подготовка МЕХ-файла

Проведение вычислений

=========================================================

Еще один пример

Вычисление определённого интеграла по квадратурной формуле Гаycса

Для вычисления воспользуемся подпрограммой на Фортране QS12R. Подпрограмма входит в Библиотеку численного анализа НИВЦ МГУ. Ссылка на текст, приведенный ниже http://num-anal.srcc.msu.ru/lib_na/cat/q/qs12r.htm

Подпрограмма: QS12R

Назначение

Вычисление определенного интеграла по квадратурной формуле Гаусса.

Математическое описание

QS12R вычисляет интеграл с заданной абсолютной погрешностью на отрезке интегрирования по формуле

где yi и Ci соответственно узлы и веса Гаусса для отрезка [-1, 1].

Л.Г. Васильева. Набор стандартных программ численного интегрирования с фиксированным распределением узлов. Сб. «Численный анализ на ФОРТРАНе», вып. 8, Изд-во МГУ, 1974.

Использование

SUBROUTINE QS12R (RINT, A, В, F, E, N, IERR)

Параметры

12

EXTERNAL F

INTEGER N

REAL RINT, A, B, F, E

A = 1.

B = EXP(3.)

E = 0.0000001

N = 2

CALL QS12R {RINT, A, B, F, E, N, IERR)

Результаты:

RINT = 1.9999999999 N = 32

Для дальнейшей работы нам надо разобраться с таким интересным средством работы в Matlab как function_handle.

Массивы указателей на функции (function_handle arrays)

Одним из типов данных, заслуживающим внимания, является класс function_handle. С терминомhandle (произносится – хэндл) многие программисты сталкивались при работе с файлами при создании подпрограмм обработки событий в различных визуальных средах. Термин handle пока не получил единообразной трактовки в отечественной литературе. Под этим термином скрывается системный способ идентификации некоторых программных компонентов. Так, например, при работе с файлами вместо их идентификаторов (имен дисковых файлов или их полных спецификаций) в операторах обмена указывается тем или иным способом созданный handle. Им может быть конкретное число, которое операционная система присваивает открываемым файлам (int fl=open...), или указатель на блок управле-

ния файлом (см. в C++ описание типа FILE *fl;). Иногда этот термин именуют дескриптором (т. е. описателем). Но, видимо, наиболее точно термину handle врусскомязыке соответствуюттермины указательили ссылка.

Указатель на функцию (function_handle) обеспечивает возмож-

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

»f_h=@funl;

»f_h=str2func('funl')

Врезультате и той, и другой операции создается указатель f_h, который «смотрит» на функцию с именем funl. Первый вариант очень напоминает присвоение адреса переменной указателю соответствующего типа в C++ (int x; int *px=&x;). Указатель f_h может быть передан в каче-

стве аргумента другой функции:

» y = integral(a,b,f_h);

14

Для того чтобы функция integral могла воспользоваться значением функции funl (x), она использует системную процедуру feval:

»z = feval(f_h, x);

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

Дополнительную информацию об использовании этого типа данных

можно найти в справке Matlab, сделав поиск по ключевому слову function_handle.

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

Теперь, чтобы получить пятую степень числа, например, 2, надо запи-

сать

Этот же результат можно получить, используя функцию Matlab feval, что мы и будем делать в дальнейшем

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

соответствует типу Real*8. Так как в используемой нами процедуре на Фортране используется тип Real, то и в Matlab надо использовать соответствующий тип данных – single.

Вызовем Matlab Editor и создадим файл Test.m

15

Напишем теперь функцию-обертку, куда будем передавать объявленные аргументы. Перейдем в Фортран и создадим файл QS12R_wrapper.F.

#include "fintrf.h"

! Gateway routine

subroutine mexFunction(nlhs, plhs, nrhs, prhs)

!Declarations implicit none

!mexFunction arguments:

mwPointer plhs(*), prhs(*) ! массивы произвольного размера integer nlhs, nrhs

!Function declarations: mwPointer mxGetPr mwPointer mxGetData

mwPointer mxCreateNumericArray integer mxIsClass

integer mxClassIDFromClassName integer mexCallMATLAB

integer mexPrintf

!Variables

real RESULT, A, B, err ! MATLAB использует тип double, но мы используем single

integer N, IERR, status, classid

mwPointer input(2), output(1) ! массивы указателей

16

character*120 line integer k

! Program !-----------------------------------------------------------------------

! Check for proper number of arguments. if(nrhs .ne. 5) then

endif

!Validate inputs

write (line, *) 'A number of nodes = ', N ! итоговое количество узлов для заданной точности

k = mexPrintf (line//achar(13))

plhs(1) = mxCreateNumericArray (2, (/1,1/), classid, 0) ! в plhs(1) нет объ-

екта, надо его создать

call mxCopyReal4ToPtr (RESULT, mxGetPr(plhs(1)),1) !запишемрезуль-

тат на выход, не забываем о соответствии типов

! не забываем за собой убирать call mxDestroyArray (input (2))

ourfunc = out return

END FUNCTION

END SUBROUTINE

Теперь в рабочем каталоге Matlab должны находится следующие файлы

18