Неверный ответ с использованием sgeev по сравнению с matlab

Мне было интересно, есть ли у кого-нибудь опыт использования sgeev для вычисления e-vale / e-vecs в фортране. В настоящее время у меня проблема с диагонализацией матрицы, и я не уверен, почему. матрица

   1    2   4   4   22  -3  22
   3    3   8   -3  -22 -2  14
   8  -2.3  16  2.5 22  1   7
  -6   17   22  -9  22  17  -6
   7   1    22  2.5 16  -2.3 8
  14  -2   -22  -3  8   3   3
  22    -3  22  4   4   2   1 

Это определенно диагональю, так как он отлично работает в Matlab, но я не могу заставить его работать в Fortran и понятия не имею, почему,

мой призыв к использованию sgeev является правильным, так как он был протестирован с другими матрицами и работает нормально, возвращая масштабированные результаты и т. д.

Я знаю, что матрица обладает тем свойством, что первый столбец является обратным последнему столбцу и т. Д., Но я подумал, что использование общей формы матрицы в фортране было бы хорошо. Если бы кто-то мог пролить свет на это, это было бы очень ценно

    'program trial
    implicit none
    integer, parameter :: M=7
    integer, parameter :: N=6
    real :: qqq(7,7), ttt(7,7)
    character*1 :: jobvl, jobvr
    real :: wr(M)
    real :: diag(M,M)
    real :: wi(M)
    real :: vl(M,M)
    integer :: LDVL=M
    integer :: IHI, ILO
    real :: vr(M,M)
    integer :: LDVR=M
    real :: work(4*M)
    integer :: lwork=4*M
    integer :: info, infonow, check
    character (len=40) :: print_file
    integer :: filenumber=1
    integer :: r, rr, rrr

    qqq(1,1)=1
    qqq(1,2)=2
    qqq(1,3)=4
    qqq(1,4)=4
    qqq(1,5)=22
    qqq(1,6)=-3
    qqq(1,7)=22

    qqq(2,1)=3
    qqq(2,2)=3
    qqq(2,3)=8
    qqq(2,4)=-3
    qqq(2,5)=-22
    qqq(2,6)=-2
    qqq(2,7)=14

    qqq(3,1)=8
    qqq(3,2)=-2.3
    qqq(3,3)=16
    qqq(3,4)=2.5
    qqq(3,5)=22
    qqq(3,6)=1
    qqq(3,7)=7

    qqq(4,1)=-6
    qqq(4,2)=17
    qqq(4,3)=22
    qqq(4,4)=-9
    qqq(4,5)=22
    qqq(4,6)=17
    qqq(4,7)=-6

    qqq(5,1)=7
    qqq(5,2)=1
    qqq(5,3)=22
    qqq(5,4)=2.5
    qqq(5,5)=16
    qqq(5,6)=-2.3
    qqq(5,7)=8

    qqq(6,1)=14
    qqq(6,2)=-2
    qqq(6,3)=-22
    qqq(6,4)=-3
    qqq(6,5)=8
    qqq(6,6)=3
    qqq(6,7)=3

    qqq(7,1)=22
    qqq(7,2)=-3
    qqq(7,3)=22
    qqq(7,4)=4
    qqq(7,5)=4
    qqq(7,6)=2
    qqq(7,7)=1


    do rr=1,7
    do r=1,7
    ttt(r,rr)=qqq(r,rr)
    end do
    end do

    jobvl='V'
    jobvr='V'

    call SGEEV(jobvl,jobvr,M,qqq,M,wr,wi,vl,LDVL,vr&
    ,LDVR,work,lwork,info)

    do rr=1,N+1
    do r=1,N+1
    if (r==rr) then
    diag(r,rr)=wr(r)
    else
    diag(r,rr)=0
    end if
    end do
    end do


    write(*,*) wr




    vl=transpose(vl)`

1 ответ

Решение

Есть несколько проблем, которые необходимо решить:

1) Последние четыре собственных значения вашей матрицы сложны. Поскольку вы игнорируете их, вы не можете получить правильный результат.

2) Собственные векторы для сопряжения собственных значений также являются сложными и представляют собой комбинацию сообщенных собственных векторов собственных значений.

Обе проблемы могут быть решены путем перехода на сложный вариант CGEEV,

Последнее, что диагональная форма qqq = vr' * diag(wr) * vr только в общем случае, если ваша матрица qqq эрмитова или реальная симметричная.

В вашем случае вы должны рассчитать обратную матрицу vr, Другая возможность - это генерация ортонормированной системы из ваших собственных векторов, которая matlab может сделать для вас по умолчанию.

Другие вопросы по тегам