LXF138:GNU Scientific Library
|
|
|
- GNU Scientific Library
Научная библиотека GNU
- Василий Олоничев напоминает, что существительное «компьютер» происходит от глагола «вычислять», и предлагает подходящую библиотеку для этих самых вычислений.
Будучи аспирантом, я бережно брал с полки и нес к дисководу размером с хорошую тумбочку пакет дисков емкостью 7 МБ и весом килограммов эдак 5. На нем находилась «Научная библиотека Фортрана» – как мне тогда казалось, самое главное, ради чего и сделаны компьютеры. С тех пор многое изменилось: кое-кто уже и не подозревает, что слова «компьютер» и «калькулятор» – синонимы. Они думают, что калькулятор – чтобы считать, а компьютер – чтобы бродить по Интернету, играть в игрушки и иногда набирать что-нибудь в текстовом процессоре. Хотя, как ни странно, и сегодня встречаются чудаки, пытающиеся использовать компьютер по его первоначальному назначению. Именно для них и предназначается библиотека GSL – GNU Scientific Library (http://www.gnu.org/software/gsl), правда, написана она уже на C.
GSL создана, поддерживается и развивается в основном учеными-физиками из Национальной лаборатории в Лос-Аламосе (США) и в первую очередь для своих собственных нужд. Но распространяется она по лицензии GPL и поэтому доступна всем. Обновления происходят приблизительно раз в полгода. Документация к библиотеке просто прекрасная, но только на английском языке. Приведем из нее пару абзацев, дающих общее представление о GSL:
- «Библиотека спроектирована так, что различные алгоритмы могут быть подключены или изменены во время исполнения программы, без ее перекомпиляции.
- Интерфейс достаточно простой, и его можно использовать в языках более высокого уровня, таких как GNU Guile или Python.
- Библиотека может быть использована в многопоточных приложениях без всяких ограничений.
- Там, где это возможно, функции основаны на надежных общедоступных пакетах Фортрана, таких как FFTPACK и QUADPACK.
- Библиотека не имеет зависимостей от других пакетов и поэтому без проблем компилируется и устанавливается».
- В GSL реализованы практически все алгоритмы вычислительной математики: работа с полиномами, линейная алгебра, численное интегрирование и дифференцирование, поиск корней, решение систем дифференциальных уравнений, интерполяция, аппроксимация, статистика, оптимизация и так далее.
Сплайны
Зачем, спрашивается, в наши дни самому писать расчетные программы на C, если есть Scilab, Octave, R и прочие? Во-первых, физики не шутят, когда принимаются за вычисления, и поэтому всегда есть потребность в эффективных программах. Во-вторых, есть задачи, где, помимо интенсивных расчетов, требуется доступ к системным ресурсам – в частности, средствам межпроцессного взаимодействия. Такими задачами могут быть приложения реального времени, сделанные в стиле «Unix way» и выполняющиеся на промышленных компьютерах. И, наконец, есть превеликое множество толковых инженеров и экономистов, предпочитающих все делать своими руками, и для них GSL подходит идеально. Да и задачи встречаются самые разнообразные: для одних хорошо подходит Scilab, а для других лучше написать программу.
Рассмотрим пример. Пусть управляющая система снимает показания с датчика, подключенного через преобразователь–усилитель–АЦП. Статическая характеристика датчика нелинейна. В этом случае удобно использовать интерполяцию сплайнами. GSL содержит большой набор средств для этих целей, буквально на все случаи жизни. В простейшем случае, как наш, можно ограничиться использованием следующих функций:
- подготовительные
gsl_interp_accel * gsl_interp_accel_alloc (void) gsl_spline * gsl_spline_alloc (const gsl interp type * T, size t size) int gsl_spline_init (gsl spline *spline, const double xa[], const double ya[], size t size)
- основная рабочая
double gsl_spline_eval (const gsl spline * spline, double x, gsl interp accel * acc)
- завершающие
void gsl_interp_accel_free (gsl interp accel* acc) void gsl_spline_free (gsl spline * spline)
Первая из них создает в памяти объект-акселератор. Это своего рода итератор, сохраняющий предыдущие обращения и позволяющий ускорить последующие вычисления. Следующая функция создает в памяти сам сплайн. Сплайны могут быть самых разнообразных типов. Кстати, этот тот самый случай, когда алгоритм может быть изменен во время работы программы.
Данные функции сами просятся, чтобы их инкапсулировали в класс, конструктор которого получает в качестве параметра спецификацию текстового файла, содержащего градуировочную кривую в виде «вход АЦП – значение физической величины, измеряемой датчиком». Он, а также простенькая программа, его использующая, имеются на LXFDVD. А вот фрагмент кода этого класса:
class Sensor{ protected: int N; double *x; double *y; gsl_interp_accel *acc; gsl_spline *spline; public: Sensor(char* f_calibr){ // Чтение данных из текстового файла с гра дуировочной кривой ... x=new double[N*sizeof(double)]; y=new double[N*sizeof(double)]; ... acc = gsl_interp_accel_alloc(); spline = gsl_spline_alloc(gsl_interp_cspline, N); gsl_spline_init(spline, x, y, N); } ~Sensor(){ gsl_spline_free(spline); gsl_interp_accel_free(acc); delete[] x; delete[] y; } double inline getVal(double xi){ return gsl_spline_eval(spline, xi, acc); } };
Здесь используется сплайн типа gsl_interp_cspline, соответствующий «классическому» кубическому сплайну, у которого вторая производная в первой и последней точке интервала равна 0.
В проекте необходимо указать на библиотеки GSL, которые в моей системе размещаются в каталоге /usr/local/lib. Если вы работаете в Code::Blocks, то на вкладке Linker settings следует прямо указать спецификацию библиотеки, например /usr/local/lib/lbgsl.so. Если вы используете Qmake, то в файл проекта следует добавить следующую строку: LIBS += -L/usr/local/lib -lgsl. Кроме этой основной библиотеки, если вы прямо или опосредованно используете матрицы, потребуется добавить библиотеку libgslcblas.
Метод наименьших квадратов
Приведем еще один пример из области промавтоматики, в котором покажем применение GSL для идентификации объекта управления методом наименьших квадратов. Полный текст программы также имеется на диске, прилагаемом к журналу, а здесь представлен ее фрагмент.
int M; // порядок объекта int N; // кол-во строк в используемых векторах и матрице // вспомогательные массивы для размещения // считанных из файла экспериментальных данных std::vector<double> yV,uV,tV; gsl_matrix *F; gsl_vector *T,*R,*S, *Y; // считывание данных эксперимента из текстового файла ... N=yV.size()M; // выде ление памя ти под век торы и мат рицы F=gsl_matrix_calloc(N,2*M); T = gsl_vector_calloc(2*M); R = gsl_vector_calloc(2*M); S = gsl_vector_calloc(N); Y = gsl_vector_calloc(N); // заполнение вектора Y и матрицы F данными эксперимента for(int i=0; i<N; i++) gsl_vector_set(Y,i,yV[i+M]); for(int i=0; i<N; i++){ for(int j=0;j<M;j++){ gsl_matrix_set(F,i,j,(1.0)*yV[i+Mj1]); } } for(int i=0; i<N;i++){ for(int j=0;j<M;j++){ gsl_matrix_set(F,i,j+M,uV[i+Mj1]); } } // определение параметров объекта gsl_linalg_QR_decomp(F,T); gsl_linalg_QR_lssolve(F,T,Y,R,S); // вывод результатов for(int j=0; j<M+M; j++) printf(“%f\n”,gsl_vector_get(R,j)); // освобождение динамической памяти gsl_matrix_free(F); gsl_vector_free(T); gsl_vector_free(R); gsl_vector_free(S); gsl_vector_free(Y);
Здесь функции gsl_matrix_calloc() и gsl_vector_calloc() выделяют память под матрицы и векторы заданных размеров, а gsl_matrix_free() и gsl_vector_free(), наоборот, освобождают память. Функции gsl_matrix_set() и gsl_vector_set() позволяют задавать значения элементов матриц и векторов. Нетрудно дога-даться, что существуют симметричные функции gsl_matrix_get() и gsl_vector_get(). Определение параметров модели методом наименьших квадратов осуществляется при помощи функций gsl_linalg_QR_decomp() и gsl_linalg_QR_lssolve(). Первая выполняет предварительную декомпозицию матрицы F, а вторая – решает задачу. Векторы T и S – вспомогательные. Вектор Y содержит значения на выходе объекта, а матрица F – значения на входе объекта и значения на выходе объекта в предыдущие моменты времени. Решение, содержащее коэффициенты числителя и знаменателя разностной передаточной функции, помещаются в вектор R.
Любой специалист без труда подберет десяток-другой задач из своей предметной области, требующих вычислений, для которых программу просто так с ходу не напишешь. Да даже если и напишешь, то она, скорее всего, будет неэффективна и непредсказуема во нештатных ситуациях. Здесь, как и во многих других случаях, свободный код помогает решить все проблемы.
