2016-07-19 6 views
1

В следующей программе, если я изменил «L» (прочитайте нижнюю треугольную часть матрицы) на «U» (верхний треугольный вход) в zheev, Я обнаружил, что выходные собственные векторы различны. Кто-нибудь поможет?Получены различные собственные векторы для нижнего и верхнего входных сигналов в zheev в лапаке

program zheev_test 
    INTEGER   N 
    PARAMETER  (N = 4) 
    INTEGER   LDA 
    PARAMETER  (LDA = N) 
    INTEGER   LWMAX 
    PARAMETER  (LWMAX = 1000) 

    INTEGER   INFO, LWORK 

    DOUBLE PRECISION W(N), RWORK(3*N-2) 
    COMPLEX*16  A(LDA, N), WORK(LWMAX),vect(n,n) 


a(1,:)=[(9.14,0.00),(-4.37,-9.22),(-1.98,-1.72),(-8.96,-9.50)] 
a(2,:)=[(-4.37,9.22),(-3.35,0.00),(2.25,-9.51),(2.57,2.40)] 
a(3,:)=[(-1.98,1.72),(2.25,9.51),(-4.82,0.00),(-3.24,2.04)] 
a(4,:)=[(-8.96,9.50),(2.57,-2.40),(-3.24,-2.04),(8.44,0.00)] 

WRITE(*,*)'ZHEEV Example Program Results' 

LWORK = -1 
CALL ZHEEV('V', 'L', N, A, LDA, W, WORK, LWORK, RWORK,INFO) 
LWORK = MIN(LWMAX, INT(WORK(1))) 

    CALL ZHEEV('V', 'L', N, A, LDA, W, WORK, LWORK, RWORK,INFO) 

    IF(info>0) stop 'The algorithm failed to compute eigenvalues.' 

    write(*,*) w 
    write(*,*) a(:,1) 

    end program 
+1

Интересное сочетание Fortran 2003 и некоторые старые 77 стиль. –

ответ

2

Я предполагаю, что это, вероятно, связано с различными фазами собственных векторов (которые определяются только до сложных мультипликативных факторов). Чтобы избежать этого, может быть полезно разделить весь собственный вектор, например, первым элементом таким образом, что

write(*,*) a(:,1)/a(1,1) 

Тогда мы получим (с ifort14)

For mode = 'L': 
eigenvalues = -16.0047467285727 -6.76497034056390 6.66571163878665 25.5140052777621  

a(:,1) = (0.344765913149177, 0.000000000000000E+000) 
      (0.441821763159905, -0.538931962007287) 
      (-0.479504340748548, -0.374404626123115) 
      (0.100522573509652, -0.123589281135427) 

a(:,1)/a(1,1) = (1.00000000000000, 0.000000000000000E+000) 
        (1.28151231403417, -1.56318226788881) 
        (-1.39081133737566, -1.08596764309792) 
        (0.291567610589614, -0.358473028863359) 

For mode = 'U': 
eigenvalues = -16.0047467285727 -6.76497034056390 6.66571163878665 25.5140052777620 

a(:,1) = (0.217545360791378,  0.267465046067210) 
      (0.696883676007484,  2.696699658290103E-003) 
      (-1.210616849344158E-002, -0.608240641147479) 
      (0.159308186219588,  0.000000000000000E+000) 

a(:,1)/a(1,1) = (1.00000000000000, 0.000000000000000E+000) 
        (1.28151231403417, -1.56318226788881) 
        (-1.39081133737566, -1.08596764309792) 
        (0.291567610589614, -0.358473028863359) 

, что согласуется с результатом, полученным из другое программное обеспечение (но только в пределах одинарной точности, потому что ваша a матрица установлена ​​с одинарной точностью литералов!)

eigenvalues = -16.004746472094745 -6.764970154793344 6.665711453507093 25.514005173380916 

V[1,1]/V[1,1] = 1.0     - 0.0im 
V[2,1]/V[1,1] = 1.281512342601922 - 1.5631822174403054im 
V[3,1]/V[1,1] = -1.3908112850832561 - 1.0859676556672477im 
V[4,1]/V[1,1] = 0.29156759974571633 - 0.35847302874950293im 
+0

Спасибо. Да, эти векторы на самом деле одинаковы. В лапке я считаю, что собственные векторы должны быть нормировочными, а именно: nq = sum (conjg (a (:, 1)) * a (:, 1)) = (1.0,0.0). Таким образом, фазы собственных векторов могут быть фиксированными. Но я обнаружил, что в «Нижнем» случае и «Верхнем» случае подчиняются этому правителю. Итак, почему все еще проблемы фазы? – Orders

+0

Если u (:) - нормированный собственный вектор (A * u = lambda * u), то v (:) = exp (i * phi) * u (:) с произвольным вещественным phi также является нормированным собственным вектором A. Итак, в зависимости от способа вычисления собственного вектора (здесь, выбор режима «L» или «U»), я думаю, что результирующий собственный вектор отличается этим коэффициентом exp (i * phi). Этот произвольный фазовый множитель не может быть устранен нормализацией каждого собственного вектора (насколько я понимаю ...) – roygvib

+0

[Чтобы сделать фазу более ясной (и сохранить нормализацию), лучше разделить на (1,1)/abs (a (1,1)). Это делает первый элемент вектора вещественным. Аналогично, мы можем разделить на a (k, 1)/abs (a (k, 1)) для любого k (до тех пор, пока a (k, 1) не равно 0), что соответствует другому выбору «phi »(фазовый коэффициент).] – roygvib