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

10

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

Фон: Некоторое время назад я задал этот вопрос по теоретической физике 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
qubyte
источник

Ответы:

8

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

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

Джек Полсон
источник
1
2×24×4
1
Мне нравится этот ответ; несимметричный случай имеет достаточно ловушек, которые стоит рассмотреть, если бы у вашей задачи была формулировка, которая приводит к симметричным матрицам вместо несимметричных.
JM
@ MarkS.Everitt: Вы, кажется, почти там ... насколько велики матрицы? ~ 36 х 36 снова?
Джек Полсон
16×1636×36
2
@ MarkS.Everitt: Таким образом, ваша задача сейчас состоит в том, как построить матрицы 4x4. Это достаточно мало, чтобы асимптотический анализ не имел отношения к делу, поэтому ответ будет полностью зависеть от значений. Я не могу больше говорить, если вы не переведете свой связанный физический пост в линейную алгебру (что такое супероператор?!?).
Джек Полсон
7

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

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

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

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

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

Bj

Bj=γjI+Ej,

γjjEj

Джефф Оксберри
источник
1
O(n2)O(n3)
Без сомнения, они знают, что делают; Я не беспокоюсь о реализации LAPACK. Меня больше удивляет поведение компилятора Фортрана.
Джефф Оксберри
2
Да, компилятор может быть большей проблемой, чем хорошо написанный LAPACK. Это может сбить с толку, когда вы обнаружите, что ваша программа терпела неудачу только потому, что реализации абсолютного значения и деления, используемые компилятором, были испорчены ...
JM
-1

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

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

freejack
источник