Как при помощи x64 версии OpenBLAS/LAPACK решать стационарные уравнения теплопроводности?

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...

надеюсь... ))