ФОРТРАН ЗГЕЕВ, все 0 собственных значений - PullRequest
1 голос
/ 06 марта 2012

Я пытаюсь заставить рутину ZGEEV в Лапаке работать над тестовой проблемой и у меня возникают некоторые трудности.Я только начал кодировать на Фортране неделю назад, так что я думаю, что, возможно, я упускаю что-то очень тривиальное.

Я должен диагонализировать довольно большие сложные симметричные матрицы.Для начала, используя Matlab, я создал матрицу 200 на 200, которую, как я убедился, можно диагонализировать.Когда я запускаю код, он не выдает ошибок и INFO = 0, что говорит об успехе.Тем не менее, все собственные значения (0,0), которые я знаю, неверны.

Прикреплен мой код.

PROGRAM speed_zgeev
  IMPLICIT NONE
  INTEGER(8)  :: N
  COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: MAT
  INTEGER(8) :: INFO, I, J
  COMPLEX*16, DIMENSION(:), ALLOCATABLE :: RWORK
  COMPLEX*16, DIMENSION(:), ALLOCATABLE :: D
  COMPLEX*16, DIMENSION(1,1) :: VR, VL
  INTEGER(8) :: LWORK = -1
  COMPLEX*16, DIMENSION(:), ALLOCATABLE :: WORK
  DOUBLE PRECISION :: RPART, IPART

  EXTERNAL ZGEEV
  N = 200

  ALLOCATE(D(N))
  ALLOCATE(RWORK(2*N))
  ALLOCATE(WORK(N))
  ALLOCATE(MAT(N,N))

  OPEN(UNIT = 31, FILE = "newmat.txt")
  OPEN(UNIT = 32, FILE = "newmati.txt")
  DO J = 1,N
     DO I = 1,N
    READ(31,*) RPART
    READ(32,*) IPART
    MAT(I,J) = CMPLX(RPART, IPART)
     END DO
  END DO

  CLOSE(31)
  CLOSE(32)

  CALL ZGEEV('N','N', N, MAT, N, D, VL, 1, VR, 1, WORK, LWORK, RWORK, INFO)
  INFO = WORK(1)

  DEALLOCATE(WORK)
  ALLOCATE(WORK(INFO))

  CALL ZGEEV('N','N', N, MAT, N, D, VL, 1, VR, 1, WORK, LWORK, RWORK, INFO)

  IF (INFO .EQ. 0) THEN
     PRINT*, D(1:10)
  ELSE
     PRINT*, INFO
  END IF

  DEALLOCATE(MAT)
  DEALLOCATE(D)
  DEALLOCATE(RWORK)
  DEALLOCATE(WORK)


END PROGRAM speed_zgeev

Я пробовал тот же код на меньших матрицах, размер 30к 30 и они работают нормально.Любая помощь будет оценена!Спасибо.

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

1 Ответ

3 голосов
/ 06 марта 2012

Может быть LWORK = WORK (1) вместо INFO = WORK(1)?Также измените ALLOCATE(WORK(INFO)).

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...