1

Итак - x64-версия библиотеки OpenBLAS.dll прекрасно подключилась в x64 c#-проект (Подключение фортрановской LAPACK x64 в x64 проект на C#). Как же теперь с её помощью решать стационарные (и даже нестационарные, и даже неодномерные) уравнения теплопроводности?

1 ответ 1

2

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

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

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.