1

Fortran で SCALAPACK と MPI を使用して、エルミート行列の固有値と固有ベクトルを見つけようとしています。バグをつぶすために、このプログラムをできるだけ単純にしましたが、それでもセグメンテーション違反が発生します。同様の質問を持つ人々への回答によると、すべての整数を整数 * 8 に変更し、すべての実数を実数 * 8 または実数 * 16 に変更しようとしましたが、それでもこの問題が発生します。最も興味深いのは、セグメンテーション違反のバックトレースすら得られないことです。プログラムはバックトレースを提供しようとするとハングアップし、手動で中止する必要があります。

また、私の知識不足をお許しください。私はほとんどのプログラム関連のことには慣れていませんが、最善を尽くしました。これが私のコードです：

    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

問題は、2 番目の PZHEEVX 関数にあります。このコードは、正常に動作する別のより複雑なコードの単純なバージョンであるため、正しく使用していると確信しています。この目的のために、私は 1 つのプロセッサのみを使用しています。

ヘルプ！

4

2 に答える 2

2

このページの設定によると、 LWORK = -1 は PZHEEVX ルーチンにすべての作業配列の必要なサイズを返すように要求しているようです。たとえば、

LWORK = -1 の場合、LWORK はグローバル入力であり、ワークスペース クエリが想定されます。このルーチンは、すべての作業配列の最適なサイズのみを計算します。これらの各値は、対応する作業配列の最初のエントリに返され、PXERBLA によってエラー メッセージは発行されません。

LRWORK = -1 についても同様の説明が見られます。アイワークに関しては、

IWORK (ローカル ワークスペース) INTEGER 配列 返されると、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) )

このページの例も、この方法で作業配列を割り当てているようです。NPROW=INT(SQRT(REAL(NPROCS)))もう 1 つの懸念事項は、INT() がいくつかの場所で使用されていることです(たとえば、

(*) より正確には、割り当てられた配列のサイズが 0 になるため、-1 を使用した配列の割り当ては無効です (@francescalus さんに感謝)。これは、size(a) または a(:) を出力することで確認できます。この種のエラーを防ぐには、-fcheck=all (gfortran の場合) や -check (ifort の場合) などのコンパイラ オプションを付けると非常に便利です。

于 2015-08-30T06:53:24.657 に答える
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)

これはすでに不整合であり、そのサブルーチンでの計算を台無しにする可能性があります (セグメンテーション違反がなくても)。さらに、サブルーチンと scalapack は、メイン プログラムで割り当てたAサイズ(8,8)ではなくsize であると見なし(4,4)、サブルーチンが使用可能なメモリをオーバーランできるようにします。

于 2015-08-30T20:47:13.963 に答える