2015-08-30 8 views
1

Я пытаюсь найти собственные значения и собственные векторы эрмитовой матрицы, используя SCALAPACK и MPI в Fortran. Для устранения ошибок я сделал эту программу максимально простой, но я все еще получаю ошибку сегментации. По ответам, предоставленным людям с похожими вопросами, я попытался изменить все мои целые числа на integer * 8, а все мои реалы на реальные * 8 или реальные * 16, но я все еще получаю эту проблему. Самое интересное, что я даже не получаю обратную трассировку для ошибки сегментации: программа зависает, пытаясь дать мне обратную трассировку и ее нужно прервать вручную.Ошибка сегментации с использованием SCALAPACK в Фортране? Нет обратной линии?

Кроме того, пожалуйста, простите мне недостаток знаний - я не знаком с большинством программных вещей, но я сделал все возможное. Вот мой код:

PROGRAM easydiag 
    IMPLICIT NONE 
    INCLUDE 'mpif.h' 
    EXTERNAL BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT, BLACS_GRIDINFO 
    EXTERNAL BLACS_GRIDINIT, BLACS_PINFO,BLACS_SETUP, DESCINIT 
    INTEGER,EXTERNAL::NUMROC,ICEIL 
    REAL*8,EXTERNAL::PDLAMCH 

    INTEGER,PARAMETER::XNDIM=4 ! MATRIX WILL BE XNDIM BY XNDIM 
    INTEGER,PARAMETER::EXPND=XNDIM 
    INTEGER,PARAMETER::NPROCS=1 

    INTEGER COMM,MYID,ROOT,NUMPROCS,IERR,STATUS(MPI_STATUS_SIZE) 
    INTEGER NUM_DIM 
    INTEGER NPROW,NPCOL 
    INTEGER CONTEXT, MYROW, MYCOL 

    COMPLEX*16,ALLOCATABLE::HH(:,:),ZZ(:,:),MATTODIAG(:,:) 
    REAL*8:: EIG(2*XNDIM) ! EIGENVALUES 
    CALL MPI_INIT(ierr) 
    CALL MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr) 
    CALL MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr) 
    ROOT=0 

    NPROW=INT(SQRT(REAL(NPROCS))) 
    NPCOL=NPROCS/NPROW 
    NUM_DIM=2*EXPND/NPROW 

    CALL SL_init(CONTEXT,NPROW,NPCOL) 
    CALL BLACS_GRIDINFO(CONTEXT, NPROW, NPCOL, MYROW, MYCOL) 

    ALLOCATE(MATTODIAG(XNDIM,XNDIM),HH(NUM_DIM,NUM_DIM),ZZ(NUM_DIM,NUM_DIM)) 
    MATTODIAG=0.D0 

    CALL MAKEHERMMAT(XNDIM,MATTODIAG) 

    CALL MPIDIAGH(EXPND,MATTODIAG,ZZ,MYROW,MYCOL,NPROW,NPCOL,NUM_DIM,CONTEXT,EIG) 


    DEALLOCATE(MATTODIAG,HH,ZZ) 




    CALL MPI_FINALIZE(IERR) 


END 

!**************************************************** 
SUBROUTINE MAKEHERMMAT(XNDIM,MATTODIAG) 
    IMPLICIT NONE 
    INTEGER:: XNDIM, I, J, COUNTER 
    COMPLEX*16:: MATTODIAG(XNDIM,XNDIM) 
    REAL*8:: RAND 

    COUNTER = 1 
    DO J=1,XNDIM 
    DO I=J,XNDIM 
     MATTODIAG(I,J)=COUNTER 
     COUNTER=COUNTER+1 
    END DO 
    END DO 




END 
!**************************************************** 
SUBROUTINE MPIDIAGH(EXPND,A,Z,MYROW,MYCOL,NPROW,NPCOL,NUM_DIM,CONTEXT,W) 
    IMPLICIT NONE 
    EXTERNAL DESCINIT 
    REAL*8,EXTERNAL::PDLAMCH 

    INTEGER EXPND,NUM_DIM 
    INTEGER CONTEXT 
    INTEGER MYCOL,MYROW,NPROW,NPCOL 
    COMPLEX*16 A(NUM_DIM,NUM_DIM), Z(NUM_DIM,NUM_DIM) 
    REAL*8 W(2*EXPND) 

    INTEGER N 
    CHARACTER JOBZ, RANGE, UPLO 
    INTEGER IL,IU,IA,JA,IZ,JZ 
    INTEGER LIWORK,LRWORK,LWORK 
    INTEGER M, NZ, INFO 

    REAL*8 ABSTOL, ORFAC, VL, VU 

    INTEGER DESCA(50), DESCZ(50) 
    INTEGER IFAIL(2*EXPND), ICLUSTR(2*NPROW*NPCOL) 
    REAL*8 GAP(NPROW*NPCOL) 
    INTEGER,ALLOCATABLE:: IWORK(:) 
    REAL*8,ALLOCATABLE :: RWORK(:) 
    COMPLEX*16,ALLOCATABLE::WORK(:) 

    N=2*EXPND 
    JOBZ='V' 
    RANGE='I' 
    UPLO='U' ! This should be U rather than L 
    VL=0.d0 
    VU=0.d0 
    IL=1 ! EXPND/2+1 
    IU=2*EXPND ! EXPND+(EXPND/2) ! HERE IS FOR THE CUTTING OFF OF THE STATE 
    M=IU-IL+1 
    ORFAC=-1.D0 
    IA=1 
    JA=1 
    IZ=1 
    JZ=1 


    ABSTOL=PDLAMCH(CONTEXT, 'U') 
    CALL DESCINIT(DESCA, N, N, NUM_DIM, NUM_DIM, 0, 0, CONTEXT, NUM_DIM, INFO) 
    CALL DESCINIT(DESCZ, N, N, NUM_DIM, NUM_DIM, 0, 0, CONTEXT, NUM_DIM, INFO) 



    LWORK = -1 
    LRWORK = -1 
    LIWORK = -1 
    ALLOCATE(WORK(LWORK)) 
    ALLOCATE(RWORK(LRWORK)) 
    ALLOCATE(IWORK(LIWORK)) 


    CALL PZHEEVX(JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, VL, & 
       VU, IL, IU, ABSTOL, M, NZ, W, ORFAC, Z, IZ, & 
       JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, & 
       LIWORK, IFAIL, ICLUSTR, GAP, INFO) 

    LWORK = INT(ABS(WORK(1))) 
    LRWORK = INT(ABS(RWORK(1))) 
    LIWORK =INT (ABS(IWORK(1))) 

    DEALLOCATE(WORK) 
    DEALLOCATE(RWORK) 
    DEALLOCATE(IWORK) 

    ALLOCATE(WORK(LWORK)) 
    ALLOCATE(RWORK(LRWORK)) 
    ALLOCATE(IWORK(LIWORK)) 


     PRINT*, LWORK, LRWORK, LIWORK 

    CALL PZHEEVX(JOBZ, RANGE, UPLO, N, A, IA, JA, DESCA, VL, & 
       VU, IL, IU, ABSTOL, M, NZ, W, ORFAC, Z, IZ, & 
       JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, & 
       LIWORK, IFAIL, ICLUSTR, GAP, INFO) 





    RETURN 
END 

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

Помощь!

+0

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

+0

И если вам удастся исправить segfault, не будет ли 'PZHEEVX' грустно, что входная матрица является неэрмитовой? По крайней мере, я не вижу, как ваш пример 'MAKEHERMMAT' создаст (реальную) симметричную матрицу. –

+0

Любые обновления по этой проблеме? :) –

ответ

2

В соответствии с this page установки LWORK = -1, кажется, чтобы запросить процедуру PZHEEVX, чтобы вернуть необходимый размер всех рабочих массивов, например,

Если LWORK = -1, то LWORK является глобальным ввода и предполагается запрос рабочей области ; процедура только вычисляет оптимальный размер для всех рабочих массивов. Каждое из этих значений возвращается в первой записи соответствующего рабочего массива, и сообщение об ошибке не равно , выпущенному PXERBLA.

Аналогичные объяснения могут быть найдены для LRWORK = -1. Что же касается IWORK,

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

, но в вашей программе рабочие массивы распределяются

LWORK = -1 
LRWORK = -1 
LIWORK = -1 
ALLOCATE(WORK(LWORK)) 
ALLOCATE(RWORK(LRWORK)) 
ALLOCATE(IWORK(LIWORK)) 

и после первого вызова PZHEEVX, размеры рабочих массивов получаются как

LWORK = INT(ABS(WORK(1))) 
LRWORK = INT(ABS(RWORK(1))) 
LIWORK =INT (ABS(IWORK(1))) 

который выглядит непоследовательным (-1 против 1). Так будет лучше изменить распределение, как (*)

allocate(WORK(1), RWORK(1), IWORK(1)) 

Пример в this page также, кажется, выделить работу массивов таким образом. Еще одна проблема состоит в том, что INT() используется в нескольких местах (например, NPROW=INT(SQRT(REAL(NPROCS))), но я думаю, что лучше использовать NINT(), чтобы избежать эффекта ошибок округления.

(*) Подробнее точно, выделение массива с -1 недопустимо, потому что размер выделенного массива становится 0 (благодаря @francescalus). Вы можете проверить это, напечатав размер (a) или (:). Чтобы предотвратить такую ​​ошибку , очень полезно прикрепить параметры компилятора, такие как -fcheck = all (для gfortran) или -check (для ifort).

+0

Исходный код выделяет работу (:) и т. Д. С отрицательными скалярами, например allocate (work (-1)). Хотя gfortran принимает это, я не уверен, что это хорошо для других компиляторов ... – roygvib

+0

'allocate (a (-1))' имеет определенный эффект выделения массива нулевого размера. – francescalus

+0

@francescalus Поскольку я никогда не выделял массивы с -1, я тестировал это с помощью простой программы. Тогда как gfortran, так и ifort не делали «beep» при распределении с -1 и просто отключались во время выполнения (когда не была дана опция проверки границ). Хуже того, gfortran корректно печатает значение (-1) перед segfault. Я думаю, что компилятор должен остановить программу, когда массивы выделены с -1 ... (не -1: -1). – roygvib

1

В вашем коде есть своя часть, которая может легко отвечать за segfault. В главной программе вы установите

EXPND=XNDIM=4 
NUM_DIM=2*EXPND !NPROW==1 for a single-process test 
ALLOCATE(MATTODIAG(XNDIM,XNDIM)) ! MATTODIAG(4,4) 

Затем вы передаете свой MATTODIAG, эрмитову матрицу, чтобы

CALL MPIDIAGH(EXPND,MATTODIAG,ZZ,MYROW,...) 

, который в свою очередь определяется как

SUBROUTINE MPIDIAGH(EXPND,A,Z,MYROW,...) 

COMPLEX*16 A(NUM_DIM,NUM_DIM) ! A(8,8) 

Это уже несоответствие, который может испортить вычисления в этой подпрограмме (даже без segfault). Кроме того, подпрограмма наряду с scalapack считает, что A имеет размер (8,8), а не (4,4), который вы выделили в основной программе, что позволяет подпрограмме переполнять доступную память.