К вопросу о числодробильнях

gsl, lapack, fftw3

Функции лапака можно дергать из С (С++), главное правильно описать прототипы в .h файле. Остальные легко дергаются из С, и, если очень нужно, из фортрана (но тут уже обертки надо написать, но тоже тривиальные)

Я использую Си версии библиотек (blas, lapack). Если Си версия отсутствует (типичные примеры - arpack, spherepack), то создаю её с помощью f2c. Нужно чтобы собирался весь проект со всеми либами одной командой на всех платформах в том числе и под msvs.

Раньше использовал фортрановские либы и трахался в венде с mingw и линковкой этого дела под студией.

Фортрановские либы для си использовал? Обертку сам писал? Написал двойственно - запутал. (:

А как с кросскомпильность?

На каком языке работали с числодробильней?

С чем вы сталкивались/работали?

Какие впечатления и истории успехов?

Какие связки лучше всего работают?

• Только на С.
• fftw3, GSL, кое-что брал из «Numerical recipes», для двумерных вейвлетов пришлось «выпиливать» нужный фукнционал из чужого проекта
• Фурье быстрее всех делает fftw3, у GSL хорошо получаются одномерные вейвлеты (но медленновато), быстрая сортировка для медианной фильтрации - только в «numerical recipes»...

imsl + lapack + fortran.

Впечатления двоякие: описание показалось корявым, сами работают без нареканий. Заявленный функционал обеспечивают.

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

Вот это очень нужно.

Я использую Си версии библиотек (blas, lapack).

Как правильно они называются, не найду в репах.

imsl + lapack + fortran.

imsl с фортраном еще в винде юзал. Все хорошо, кроме фортрановской убогости. Причем imsl для c++ менее богат. Да и еще платный.

Обертка пишется тривиально. Функцию в h'нике обзываешь с подчеркиванием, а все данные в функцию передаешь по указателю. Например, если в фортране было void f (integer), то в Си пишешь void f_(int *)

clapack, cblas. я вообще их с собой таскаю в сорцах, чтобы геморроя со сборкой меньше было

>На каком языке работали с числодробильней?

На Common Lisp. На чистом Си ещё.

С чем вы сталкивались/работали?

Везде - мои кавайные велосипедики. Чтобы лучше знать матчасть (с)

Какие впечатления и истории успехов?

Common Lisp вполне себе подходит для числодробилок. Си - тоже, если у вас много времени. :)

Common Lisp вполне себе подходит для числодробилок.

Угу, конечно. Покажи как на нем будет выглядеть какой-нибудь GMRES ?

Как композиция s-выражений, очевидно ж.

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

Могу обучение нейронной сети показать. Мой собственный былдокод, писал ещё нубышом. Всего-то меньше 200 строчек.

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

На Common Lisp тоже можно писать один в один как в pdf-как. Вот реализация алгоритма из книжки Кнута, точь в точь как дядя завещал:

http://gigamonkeys.com/book/the-special-operators.html

Algorithm S (Selection sampling technique). To select n records at random from a set of N, where 0 < n <= N.

S1. [Initialize.] Set t <-- 0, m <-- 0. (During this algorithm, m represents the number of records selected so far, and t is the total number of input records that we have dealt with.)

S2. [Generate U.] Generate a random number U, uniformly distributed between zero and one.

S3. [Test.] If (N - t)U >= n - m, go to step S5.

S4. [Select.] Select the next record for the sample, and increase m and t by 1. If m < n, go to step S2; otherwise the sample is complete and the algorithm terminates.

S5. [Skip.] Skip the next record (do not include it in the sample), increase t by 1, and go back to step S2.

This description can be easily translated into a Common Lisp function, after renaming a few variables, as follows:

(defun algorithm-s (n max) ; max is N in Knuth's algorithm
  (let (seen               ; t in Knuth's algorithm
        selected           ; m in Knuth's algorithm
        u                  ; U in Knuth's algorithm
        (records ()))      ; the list where we save the records selected
    (tagbody
     s1
       (setf seen 0)
       (setf selected 0)
     s2
       (setf u (random 1.0))
     s3
       (when (>= (* (- max seen) u) (- n selected)) (go s5))
     s4
       (push seen records)
       (incf selected)
       (incf seen)
       (if (< selected n)
           (go s2)
           (return-from algorithm-s (nreverse records)))
     s5
       (incf seen)
       (go s2))))

Правда это не совсем числодробилка, но ведь алгоритм :)

Сокращается до

(defun algorithm-s (n max)
  (loop for seen from 0
     when (< (* (- max seen) (random 1.0)) n)
     collect seen and do (decf n)
     until (zerop n)))

Это плохой пример. Ты линейную алгебру покажи на лиспе. Вот, например, GMRES на Си

http://pastebin.com/5P6EjRfq

Интересно увидеть это на лиспе. Можешь показать что попроще, например умножение матриц на MPI :)

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

(defmacro matrix-do ((row-index rows col-index cols) &body body)
  `(dotimes (,row-index ,rows)
     (dotimes (,col-index ,cols)
       ,@body)))

(defmacro matrix-tabulate ((matrix row-arg col-arg) &body body)
  `(progn
     (matrix-do (,row-arg (matrix-rows ,matrix) ,col-arg (matrix-cols ,matrix))
       (setf (matrix-ref ,matrix ,row-arg ,col-arg) ,@body))
     ,matrix))

(defmacro matrix-create-tabulated ((row-arg rows-form col-arg cols-form) &body body)
  (with-gensyms (result)
    `(let ((,result (make-instance 'matrix :rows ,rows-form :cols ,cols-form)))
       (matrix-tabulate (,result ,row-arg ,col-arg) ,@body))))

(defmethod matrix* ((left matrix) (right matrix))
  (if (equal (matrix-cols left) (matrix-rows right))
      (matrix-create-tabulated (result-i (matrix-rows left) result-j (matrix-cols right))
        (do ((i 0 (1+ i))
             (s 0 (+ s (* (matrix-ref left result-i i)
                          (matrix-ref right i result-j))))) 
            ((= i (matrix-cols left)) s)))
      (error "Operation can not be performed")))

Образец gmres-а на CL я сразу не нашёл, а вкуривать в пдфку после суток аптайма как то не по силам сейчас.

Я же говорил, что вынос мозга будет. Умножение матриц математически описывается как nxn суммирований, что на Си записывается один в один тремя циклами. А у тебя получилась какая-то %уета в которой никто не разберется.

Раньше писал числодробильню на с. Последние пару лет - на фортране. Несравненно удобнее. Никакой особой убогости в стандартах, начиная с fortran90, не замечал. Необходимое для работы подмножество языка учится за три вечера, код получается весьма читабельный (если не юродствовать конечно), по производительности, как минимум, не медленнее Си в большинстве случаев.

Что скажут «питонцы»?

> На каком языке работали с числодробильней?

C (раньше - и на фортране). Часто использую PETSc, обеспечивающую высокоуровневый доступ к lapack и многим другим библиотеками.

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

>быстрая сортировка для медианной фильтрации

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

По теме: исключительно Ъ-Фортран :), у которого нет конкурентов для быстрых ортогональных преобразований (Фурье, Хаар, Адамар, Уолш, ...)

Никакой особой убогости в стандартах, начиная с fortran90, не замечал.

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

2. Перегрузки функций как делать?

Я же говорил, что вынос мозга будет. Умножение матриц математически описывается как nxn суммирований, что на Си записывается один в один тремя циклами. А у тебя получилась какая-то %уета в которой никто не разберется.

Да у него умножение в трёх строчках и есть. Ах, пардон, четыре: одну строчку эстетики ради на две разбили.

Да у него умножение в трёх строчках и есть. Ах, пардон, четыре: одну строчку эстетики ради на две разбили.

одну строчку эстетики ради на две разбили.

OMFG. И это эстетично О_о

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

OMFG. И это эстетично О_о

Нормально. Я бы по-другому написал, но этот код вполне читается.

Как подключать clapack, cblas - где почитать?

Всмысле как? Ключ -l у компилятора.

на мой взгляд очень простой и понятный код... :)

и эти люди потом ругают синтаксис шаблонов Си++ ... :)

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

А я и использовал алгоритм quicksort - он не выполняет полную сортировку, а лишь ищет медиану.

и эти люди потом ругают синтаксис шаблонов Си++ ... :)

А что страшным показалось, если не секрет? :)

Например для cblas:

Подключаем f2c.h, blaswrap.h И cобираем с -l... ?

слишком много слов за которыми ускользает суть

f2c.h нужен только для библиотек, которые конвертировали автоматически с помощью f2c да и то не всегда. -l нужно для линковки, подключать ничего не надо

Не собирается, можешь пример привести, ошибку не найду.

Так что пишет?

А откуда у тебя f2c_dswap в main ? Покажи код.

Intel MKL и ACML идут с готовыми обертками

питонцы могут сказать только одно: numpy

Intel MKL и ACML идут с готовыми обертками

MKL ещё и оптимизирована очень прилично, что большой плюс.

blaswrap убери нахер и сам сделай определение dswap_. можно прям скопипастить из dswap.c

да, f2c таки нужен

#include "f2c.h"

int dswap_(integer *n, doublereal *dx, integer *incx,
        doublereal *dy, integer *incy)


int main()
{
// твой код
}

// gcc test.c -lcblas -lf2c

numpy - обертка вокруг lapack. Причем там не все. Еще есть scipy. Как оно не знаю.

>MKL ещё и оптимизирована очень прилично, что большой плюс.

ACML тоже, только для AMD. Еще есть Goto BLAS, но лапака в нем нет. все остлальное фуфло

>кое-что брал из «Numerical recipes»

ай-яй-яй, нехорошо лицензию нарушать