Итак - x64-версия библиотеки OpenBLAS.dll прекрасно подключилась в x64 c#-проект (Подключение фортрановской LAPACK x64 в x64 проект на C#). Как же теперь с её помощью решать стационарные (и даже нестационарные, и даже неодномерные) уравнения теплопроводности?
1 ответ
OpenBLAS/LAPACK x64 успешно подключен и проверен (Подключение фортрановской LAPACK x64 в x64 проект на C#). Приступим-с...
Пример, приведенный в собственном ответе на вопрос о подключении x64 LAPACK был таков - стационарное уравнение теплопроводности d2u/dx2 = 0 с граничными условиями u(0) = 0 = u0, u(1) = 1 = u1. Задача решалась на трех точках и вся система уравнений выглядела так:
v[0] = 0.0 //левое граничное условие
v[1] - 2*v[1] + v[2] = 0.0 //конечно-разностный аналог второй производной d2u/dx2 с точностью до множителя 1/dx^2
v[2] = 1.0 //правое граничное условие
Давайте же решим задачу для произвольного N (M = N - 1). После дискретизации дифференциального уравнения в частных производных мы получаем систему линейных уравнений (вот зачем нам и нужен Linear Algebra PACKage = LAPACK) с матрицей коэффициентов и вектором правых частей:
|1 0 0 0 0 0 0 | |u[0] | |u0|
|1 -2 1 0 0 0 0 | |u[1] | |0 |
|0 1 -2 1 0 0 0 | |u[i-1]| |0 |
|0 0 1 -2 1 0 0 |*|u[i] | = |0 |
|0 0 0 1 -2 1 0 | |u[i+1]| |0 |
|0 0 0 0 1 -2 1 | |u[M-1]| |0 |
|0 0 0 0 0 0 1 | |u[M] | |u1|
Матрица, как видим - трёхдиагональная. Для работы с такими матрицами в фортрановском LAPACK'е существует специальная функция DGTSV с большим набором аргументов, каждый из которых передается только по ссылке, т.е. в Фортране вообще нельзя CALL MYSUB (0) - только CALL MYSUB (x). Для совместимости фортрановские библиотеки пишутся на диалекте, который был ещё при динозаврах, а исходный код без компиляции исполняется на абаках и счётах... А тому человеку, который двинул в массы идею, что Фортран - лучший язык для математических расчётов, хочется сказать много "доброго" и "светлого"... Просто посмотрите исходники на NetLib - https://www.netlib.org/lapack/, например DGTSV (https://www.netlib.org/lapack/explore-html/d4/d62/group__double_g_tsolve_ga2bf93f2ddefa5e671866eb2191dc19d4.html).
Тем не менее, эти библиотеки работают. А написаны они вручную в бинарном коде или на Брейнфаке - значения не имеет. Главное - чтобы работали быстро и без проблем подключались.
Разберём немного подробнее алгоритм решения. Вторая производная в точке i аппроксимируется как
d2v_dx2[i] = (dv_dx[i + 1/2] - dv_dx[i - 1/2])/dx;
d2v_dx2[i] = ((v[i + 1] - v[i])/dx - (v[i] - v[i - 1])/dx)/dx;
d2v_dx2[i] = (v[i + 1] - 2*v[i] + v[i - 1])/(dx*dx);
Для стационарной задачи d2u/dx2 = 0, т.е. dx можно и сократить, но желательно оставить - он пригодится в дальнейшем для решения задач нестационарных.
LAPACK.DGTSV требует передачи в качестве аргументов трех диагоналей матрицы коэффициентов: нижнюю (поддиагональ) - DL длиной N-1, главную - D длиной N, верхнюю (наддиагональ) - DU длиной N-1. Выделим диагонали DL, D, DU в явном виде. А затем дополним их внешними элементами и сформируем трёхстрочную матрицу c[3, N]. Почему двумерный массив c[3, N], а не три одномерных c0[N], c1[N], c2[N]? А вот заодно и проверим совместимость строк двумерного массива с одномерным.
0 |1 0 0 0 0 0 0 |
|1 -2 1 0 0 0 0 |
|0 1 -2 1 0 0 0 | DU = |0 1 1 1 1 1| 0 c|2| = |0 1 1 1 1 1 0|
|0 0 1 -2 1 0 0 | -> D = |1 -2 -2 -2 -2 -2 1| -> c|1| = |1 -2 -2 -2 -2 -2 1|
|0 0 0 1 -2 1 0 | DL = 0 |1 1 1 1 1 0| c|0| = |0 1 1 1 1 1 0|
|0 0 0 0 1 -2 1 |
|0 0 0 0 0 0 1 | 0
Зачем включать дополнительные элементы? По очень простой причине - DGTSV вызывается в единственном месте, а коэффициенты могут модифицироваться в разных местах исходного кода. Поэтому вместо смещённых индексов
c[0, i - 1] = 1;
c[1, i] = -2;
c[2, i] = 1;
лучше иметь несмещённые
c[0, i] = 1;
c[1, i] = -2;
c[2, i] = 1;
Это и эстетичней, и потенциальные проблемы из-за забывчивости устраняет. А еще будет полезно для случая граничных условий четвёртого рода - циклических.
Для расширенных диагоналей c[3, N] описание и вызов функции DGTSV модифицируется. Вместо
static extern void lapack_dgtsv(ref int n, ref int nrhs, double[] dl, double[] d, double[] du, double[] v, ref int ldb, ref int info);
lapack_dgtsv(ref size, ref nrhs, c[0], c[1], c[2], v, ref size, ref info);
имеем
static extern void lapack_dgtsv(ref int n, ref int nrhs, ref double dl, ref double d, ref double du, ref double v, ref int ldb, ref int info);
lapack_dgtsv(ref size, ref nrhs, ref c[0, 1], ref c[1, 0], ref c[2, 0], ref v[0], ref size, ref info);
Еще раз. Было:
static extern void lapack_dgtsv(... double[] dl, double[] d, double[] du, double[] v ...);
lapack_dgtsv(... c[0], c[1], c[2], v ...);
Стало:
static extern void lapack_dgtsv(... ref double dl, ref double d, ref double du, ref double v ...);
lapack_dgtsv(... ref c[0, 1], ref c[1, 0], ref c[2, 0], ref v[0] ...);
Т.е. мы передаем ссылки, которые и указатели, не на массивы, а на элементы массивов. И поскольку для DGTSV первый элемент строки c[0], т.е. c[0, 0] - не нужен, мы передаем ему указатель на второй элемент - c[0, 1] (sic!!!).
Теперь можно перейти и к исходному коду:
using System;
using System.Runtime.InteropServices;
namespace Parabolic
{
class Program
{
[DllImport("openblas.dll", EntryPoint = "dgtsv_")]
static extern void lapack_dgtsv(ref int n, ref int nrhs, ref double dl, ref double d, ref double du, ref double v, ref int ldb, ref int info);
static void Main(string[] args)
{
// уравнение имеет вид d2u/dx2 = 0, u(0) = u0, u(1) = u1
int size = 10 + 1; //количество точек
int maxInd = size - 1; //максимальный индекс
double dx = 1.0 / maxInd; //шаг сетки
double dxSqrInv = 1.0 / (dx*dx); //обратный квадрат шага
int nrhs = 1; //количество правых частей - нужно для lapack
int info = 0; //код успешности выполнения - нужно для lapack
double[,] c = new double[3, size]; //матрица коэффициентов
double[] v = new double[size]; //искомая величина
double v0 = 0.0; //левое граничное условие
double v1 = 1.0; //правое граничное условие
int i = 0;
c[0, i] = 0.0;
c[1, i] = 1.0;
c[2, i] = 0.0;
v[i] = v0;
for (i = 1; i < size - 1; i++)
{
double x = i / maxInd;
c[0, i] = dxSqrInv;
c[1, i] = -2.0 * dxSqrInv;
c[2, i] = dxSqrInv;
v[i] = 0.0;
}
i = maxInd;
c[0, i] = 0.0;
c[1, i] = 1.0;
c[2, i] = 0.0;
v[i] = v1;
//важно: указан c[0, 1], а не c[0, 0]
lapack_dgtsv(ref size, ref nrhs, ref c[0, 1], ref c[1, 0], ref c[2, 0], ref v[0], ref size, ref info);
for (i = 0; i < size; i++)
{
Console.Write("{0:F} ", v[i]); //ожидаем увидеть 0.00, 0.10, 0.20, ... 0.90, 1.00
}
Console.WriteLine();
Console.WriteLine ("v[0] = {0} ", v[0]);
}
}
}
Требует пояснения ожидаемое v[0] = 0.00 и имеемое v[0] = -0.00. Я вывел в конце неокруглённое значение и оно равно -2.8421709430404013E-16. Разные версии LAPACK могут выдать разные результаты, но вы таки поверьте, что -2.84/10000000000000000 таки практически равно 0.00 и пусть таки минус вас не беспокоит... )))
to be continued...
надеюсь... ))