[cuda] Решение систем ОДУ

Про куды к Eddy_Em, в жизни на них не писал;-)

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

Я только начинаю с CUDA разбираться, переношу постепенно алгоритмы на GPU.

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

В CUDA БПФ распараллелено.

dave, советую для начала посмотреть простенькие примеры на CUDA. Вы поймете, что не так уж там все и страшно. Вместо switch у вас будет номер блока/номер потока в блоке. Синхронизация внутренняя в CUDA есть. Для обмена данных между потоками (если данных немного) можно использовать shm или (если данных много) - глобальную память.

Если бы у вас было устройство с архитектурой 2.X, было бы интереснее. А так - придется все равно довольно много работы делать вручную. И жалеть, что нет атомарных операций.

Меня больше интересует автогенерация кода, который потом можно было бы запустить с использованием CUDA. Да взять ту же систему ОДУ. Человек задает ее в обычном виде на специализированном языке (типа моего http://sourceforge.net/projects/mapsim). Такая система затем обрабатывается и на выходе имеем код, который можно запустить на GPU.

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

А можешь примеры на языке на какой-нибудь pastebin.com выложить ?

Ну задать её в моём случае скорее всего либо матрицей (коэффициенты), либо значение N функций f(x) по точкам..... я пока не знаю... эх. сегодня надо было спросить.

Сначала почитаю тех ссылок, что мне накидали... может, просветление и наступит :)

Пример 1.

starttime = 0;
stoptime = 10;
dt = 0.02;

A = integ (- ka * A, 100);
B = integ (ka * A - kb * B, 0);
C = integ (kb * B, 0);

ka = 1;
kb = 1;

Пример 2 (с массивами).

starttime = 0.0;
stoptime = 10.0;
dt = 0.02;

A = ((0.1, 0.5, 0.3, 0.1),
     (0.2, 0.2, 0.4, 0.2),
     (0.6, 0.1, 0.2, 0.1));

i = low (A,0) .. high (A,0);
j = low (A,1) .. high (A,1);

A1 = min (max (A [i,j]: i): j);    // minimax: min_j max_i {a_ij}
A2 = max (min (A [i,j]: j): i);    // maximin: max_i min_j {a_ij}

Пример 3 (с конвейером).

starttime = 0;
stoptime = 25;
dt = 1;

conveyor C {
	
	inflow of Input;
	outflow of Output;
	overflow of Over;
	init = 0;
}

Input = pulse (1, 1, 3);

outflow Output from C {
	transit_time = 10;
}

overflow Over from C {
	capacity = 3;
}

Пример 4 (из жизни).

stoptime = 13;

model Fish_Bank_Model {
    
    Annual_Profit = Profit;
    
    Area = 100;
    
    Carrying_Capacity = 1000;
    
    catch_per_Ship = table ((0, -0.048), (1.2, 10.875), (2.4, 17.194), (3.6, 20.548), (4.8, 22.086), (6, 23.344), (7.2, 23.903), (8.4, 24.462), (9.6, 24.882), (10.8, 25.301), (12, 25.86)) (Density);
    
    Death_Fraction = table ((0, 5.161), (0.1, 5.161), (0.2, 5.161), (0.3, 5.161), (0.4, 5.161), (0.5, 5.161), (0.6, 5.118), (0.7, 5.247), (0.8, 5.849), (0.9, 6.151), (1, 6.194)) (Fish/Carrying_Capacity);
    
    Density = Fish/Area;
    
    reservoir Fish {
        
        inflow of Fish_Hatch_Rate;
        outflow of Fish_Death_Rate, Total_catch_per_Year;
        init = init (1000);
    }
    
    Fish_Death_Rate = max (0, Fish*Death_Fraction);
    
    Fish_Hatch_Rate = max (0, Fish*Hatch_Fraction);
    
    Fish_Price = 20;
    
    Fraction_Invested = .2;
    
    Hatch_Fraction = 6;
    
    Operating_Cost = Ships*250;
    
    Profit = Revenue-Operating_Cost;
    
    Revenue = Total_catch_per_Year*Fish_Price;
    
    Ship_building_Rate = max (0, Profit*Fraction_Invested/Ship_Cost);
    
    Ship_Cost = 300;
    
    reservoir Ships {
        
        inflow of Ship_building_Rate;
        init = init (10);
    }
    
    Total_catch_per_Year = max (0, Ships*catch_per_Ship);
    
    reservoir Total_Profit {
        
        inflow of Annual_Profit;
        init = init (0);
    }
}

Резервуар - это другое название интеграла, специализированное для стоков.

А уравнение теплопроводности в квадрате задать можно?

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

  (define (void calc_op ((vec n) u1) ((vec n) u))
   (let* (((vec n) w [L u])
          ((vec n) const_rp (w / dt 
                           + mu * (1 - theta) * L * w - sigma * (1 - theta) * w 
                           - k1 * [J (0.5 * (u + u)) (0.5 * (w + w))] 
                           - k2 * [J (0.5 * (u + u)) (l + h)] + [f x y]))
          ((vec n) w' w)
          ((vec n) u' u)
          )
     [for ((int it 0) (it < 1000) (it ++))
       (let* ((vec jac [J (0.5 * (u + u_)) (0.5 k1 * (w + w_) + k2 * (l + h))])
              (vec rp (const_rp - jac))
              )
         (begin
           (renew w_ rp)
           (L^-1 u1 w_)
           (when ((dist un u) < 1e-10)
               (break)
             )
           (swap u1 u_)
           )
         )
       ]))

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

Кстати, скобки неправильно заканчиваешь. Не тот стиль :)

Да... для тех, кто читал эту тему... Я понял, что мне нужно.

http://www.exponenta.ru/educat/class/courses/ode/theme16/theory.asp вот подходящий метод.

А именно - надо посчитать преобразование Лапласа (вроде как сгодятся и некоторые другие интегральные преобразования), тем самым перейдя к алгебраической системе, потом решить её, ну и найти ответ, выполним обратное преобразование.

Только у меня ступор в том, как закодить преобразование лапласа... численное.

Напишу в отдельную тему.