Матричная экспонента вещественной асимметричной матрицы с Fortran 95 и LAPACK

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

Фон: Некоторое время назад я задал этот вопрос по теоретической физике SE. Результат позволяет мне сформулировать основные уравнения как реальные асимметричные матрицы. В не зависящем от времени главном уравнении решается путем возведения в степень этой матрицы. В нестационарном случае это потребует интеграции. Сейчас меня интересует только независимость от времени.

После глядя на различных подпрограммах я думаю , что я должен называть ( ? Gehrd , ? Orghr , ? Hseqr ...) не ясно , если было бы проще бросить матрицу из real*8к complex*16и продолжить сложные двойные версии этих процедур, или придерживайтесь real*8и ударьте удвоение числа моих массивов и создание сложной матрицы из них позже.

Итак, какие подпрограммы я должен вызывать (и в каком порядке), и я должен использовать реальные двойные версии или сложные двойные версии? Ниже приведена попытка сделать это с настоящими двойными версиями. Я застрял в поиске собственных значений и собственных векторов L*t.

function time_indep_master(s,L,t)
  ! s is the length of a side of L, which is square.
  ! L is a real*8, asymmetric square matrix.
  ! t is a real*8 value corresponding to time.
  ! This function (will) compute expm(L*t).

  integer, intent(in)    :: s
  real*8,  intent(in)    :: L(s,s), t
  real*8                 :: tau(s-1), work(s), wr(s), wi(s), vl
  real*8, dimension(s,s) :: time_indep_master, A, H, vr
  integer                :: info, m, ifaill(2*s), ifailr(2*s)
  logical                :: sel(s)

  A = L*t
  sel = .true.

  call dgehrd(s,1,s,A,s,tau,work,s,info)
  H = A
  call dorghr(s,1,s,A,s,tau,work,s,info)
  call dhseqr('e','v',s,1,s,H,s,wr,wi,A,s,work,s,info)
  call dhsein('r','q','n',sel,H,s,wr,wi,vl,1,vr,s,2*s,m,work,ifaill,ifailr,info)

  ! Confused now...

end function

Сначала я бы очень серьезно подумал о том, является ли матрица действительно произвольной: есть ли преобразование, которое сделало бы ее эрмитовой? Гарантирует ли физика, что матрица должна быть диагонализируемой (с разумно обусловленной матрицей собственных векторов)?

Если окажется, что на самом деле не существует никакой симметрии для использования, то вам следует начать с чтения « Девятнадцати сомнительных способов вычисления экспоненциальной матрицы» , которая является стандартной ссылкой (написанной автором MATLAB и соавтором G & vL) ,

Чтобы основываться на том, что сказал Джек, стандартный подход, который, кажется, используется в программном обеспечении (например, EXPOKIT, упомянутый в вашем предыдущем вопросе), - это масштабирование и возведение в квадрат с последующим приближением Паде (методы 2 и 3) или подпространственными методами Крылова (метод 20). В частности, если вы смотрите на экспоненциальные интеграторы, вам нужно рассмотреть методы подпространства Крылова и статьи о экспоненциальных интеграторах (некоторые ссылки упоминаются вместе с методом 20 в статье Moler & van Loan).

Если вы одержимы использованием собственных векторов, рассмотрите возможность использования треугольных систем собственных векторов (метод 15); поскольку ваша матрица может быть недиагонализируемой, этот подход может быть не лучшим, но он лучше, чем попытка непосредственно вычислить собственные векторы и собственные значения (т. е. метод 14).

Приведение к форме Гессенберга не очень хорошая идея (метод 13).

Для меня неясно, будут ли вам лучше работать с реальной или сложной арифметикой, поскольку комплексная арифметика Фортрана быстрая, но может переполниться / опуститься (см. «Насколько лучше на самом деле компиляторы Фортрана?» ).

Вы можете спокойно игнорировать методы 5-7 (методы на основе решателей ODE неэффективны), методы 8-13 (дорого), метод 14 (вычисление собственных векторов больших матриц сложно без специальной структуры и склонно к численной ошибке в плохо обусловленных случаях) и метод 16 (вычисление жордановой декомпозиции матрицы численно неустойчиво). Методы 17-19 сложнее реализовать; в частности, методы 17 и 18 потребуют больше чтения. Метод 1 - это запасной вариант для масштабирования и возведения в квадрат, если приближения Паде не работают хорошо.

$B_{j}$

$\gamma_{j}$$j$$E_{j}$

У меня есть простая подпрограмма Fortran, которая вычисляет показатель произвольной матрицы. Я проверил это по команде Matlab, и это нормально. Он основан на масштабировании и возведении в квадрат. Я написал это несколько лет назад.

Я хотел бы найти другую подпрограмму, например, те, которые я загружаю с gams.nist.gov. Но пока не повезло.