线性系统的标量随处理器数变化的结果

Question

我使用过Jonathan的code并试图解决A*X=B问题，但是我的结果根据处理器的数量而变化。在这方面，有谁能帮我一下吗？

! Use MPI-IO to read a diagonal matrix distributed block-cycliCALLy,
! use Scalapack to solve ax=b.
!
Program array
  USE mpi
  IMPLICIT NONE

  INTEGER, PARAMETER :: mk=SELECTED_REAL_KIND(15,307)
  INTEGER :: N,Nb    
  ! problem size and block size
  INTEGER :: myArows, myAcols, myBcols   
  ! size of local subset of global array
  REAL(mk), DIMENSION(:), ALLOCATABLE :: myA,myB
  INTEGER,  DIMENSION(:), ALLOCATABLE :: ipiv

  INTEGER, EXTERNAL :: numroc   
  ! blacs routine
  INTEGER :: me,nprocs,ictxt,prow,pcol,myrow,mycol  
  ! blacs data
  INTEGER :: info    
  ! scalapack return value
  INTEGER, DIMENSION(9) :: desca,descb 
  ! scalapack array desc
  INTEGER :: clock
  REAL(mk) :: calctime,iotime

  CHARACTER(LEN=128) :: matsize,fnamea,fnameb
  INTEGER :: myrank
  INTEGER :: ierr
  INTEGER, DIMENSION(2) :: pdims,dims,distribs,dargs
  INTEGER :: infile
  INTEGER, DIMENSION(MPI_STATUS_SIZE) :: mpistatus
  INTEGER :: darraya,darrayb
  !New data type (handle)
  INTEGER :: locsizea, nelementsa
  INTEGER :: locsizeb, nelementsb,nelementsw
  INTEGER(KIND=MPI_ADDRESS_KIND) :: lba,locextenta
  INTEGER(KIND=MPI_ADDRESS_KIND) :: lbb,locextentb
  !INTEGER(KIND=MPI_OFFSET_KIND) :: disp
  INTEGER :: nargs
  INTEGER :: m,IPACPY,ipa
  INTEGER :: mpik
  !MPI real kind precision

  CALL MPI_Init(ierr)
  CALL MPI_Comm_size(MPI_COMM_WORLD,nprocs,ierr)
  CALL MPI_Comm_rank(MPI_COMM_WORLD,myrank,ierr)
  ! Initialize MPI (for MPI-IO)

  pdims = 0
  CALL MPI_Dims_create(nprocs, 2, pdims, ierr)
  prow = pdims(1)
  pcol = pdims(2)
  !Get the process grid from MPI_Dims_create

  nargs = command_argument_count()
  IF  (nargs /= 3) THEN
      PRINT *,'Usage: mpirun -np np ./a.out  N filenamea filenameb'
      PRINT *,'       Where np is the number of processors'
      PRINT *,'       Where N is the size of Matrix A'
      PRINT *,'       Where filenamea = name of A matrix file data.'
      PRINT *,'       Where filenameb = name of B matrix file data.'
      CALL MPI_Abort(MPI_COMM_WORLD,1,ierr)
  ENDIF

  CALL get_command_argument(1, matsize)
  CALL get_command_argument(2, fnamea)
  CALL get_command_argument(3, fnameb)
  ! get filename

  READ(matsize(1:LEN_TRIM(matsize)),*) N
  !The size of the array we'll be using

  SELECT CASE (mk)
  CASE (SELECTED_REAL_KIND(6,37))
     mpik=MPI_REAL
  CASE (SELECTED_REAL_KIND(15,307))
     mpik=MPI_DOUBLE_PRECISION
  CASE DEFAULT
     PRINT*,"The input data type is not defined!"
     CALL MPI_Abort(MPI_COMM_WORLD,1,ierr)
  END SELECT

  Nb = 64 
  !Block size (for big matrix 32 or 64 are the suitable choice
  IF (Nb > (N/prow)) Nb = N/prow
  IF (Nb > (N/pcol)) Nb = N/pcol

  dims = [N,N]
  distribs = [MPI_DISTRIBUTE_CYCLIC, MPI_DISTRIBUTE_CYCLIC]
  dargs = [Nb,Nb]
  CALL MPI_Type_create_darray(nprocs,myrank,2,dims,distribs,dargs, &
  &    pdims,MPI_ORDER_FORTRAN,mpik,darraya,ierr)
  !Creates a distributed array datatype
  CALL MPI_Type_commit(darraya,ierr)
  CALL MPI_Type_size(darraya,locsizea,ierr)
  nelementsa = locsizea/mk
  CALL MPI_Type_get_extent(darraya,lba,locextenta,ierr)

  dims = [N,1]
  dargs = [Nb,1]
  CALL MPI_Type_create_darray(nprocs,myrank,2,dims,distribs,dargs, &
  &    pdims,MPI_ORDER_FORTRAN,mpik,darrayb,ierr)
  !Creates a distributed array datatype
  CALL MPI_Type_commit(darrayb,ierr)
  CALL MPI_Type_size(darrayb,locsizeb,ierr)
  nelementsb = locsizeb/mk
  CALL MPI_Type_get_extent(darrayb,lbb,locextentb,ierr)


  ALLOCATE(myA(nelementsa))
  ALLOCATE(myB(nelementsb))
  ! Initialize local arrays    

  CALL MPI_File_open(MPI_COMM_WORLD,TRIM(fnamea),MPI_MODE_RDONLY, &
  &    MPI_INFO_NULL,infile,ierr)
  CALL MPI_File_read_all(infile,myA,nelementsa,mpik,mpistatus,ierr)
  CALL MPI_File_close(infile,ierr)
  ! read in the data


  CALL MPI_File_open(MPI_COMM_WORLD,TRIM(fnameb),MPI_MODE_RDONLY, &
  &    MPI_INFO_NULL,infile,ierr)
  CALL MPI_File_read_all(infile,myB,nelementsb,mpik,mpistatus,ierr)
  CALL MPI_File_close(infile,ierr)
  ! read in the data

  iotime = tock(clock)
  IF (myrank == 0) THEN
     PRINT *,'I/O time      = ', iotime
  ENDIF

  ! Initialize blacs processor grid
  CALL tick(clock)
  CALL blacs_pinfo(me,nprocs)
  CALL blacs_get     (-1,0,ictxt)
  CALL blacs_gridinit(ictxt,'R',prow,pcol)
  CALL blacs_gridinfo(ictxt,prow,pcol,myrow,mycol)

  myArows = numroc(N,Nb,myrow,0,prow)
  myAcols = numroc(N,Nb,mycol,0,pcol)
  myBcols = numroc(1,Nb,mycol,0,pcol)
  !NUMROC computes the NUMber of Rows Or Columns of a distributed
  !matrix owned by the process indicated by the 3rd input.
  !the 4th input element is the coordinate of the process that possesses the
  !first row or column of the distributed matrix

  ! Construct local arrays
  ! Global structure:  matrix A of n rows and n columns
  CALL descinit(desca,N,N,Nb,Nb,0,0,ictxt,max(1,myArows),info)
  CALL descinit(descb,N,1,Nb, 1,0,0,ictxt,max(1,myArows),info)
  ! Prepare array descriptors for ScaLAPACK 


  nelementsw=descb(9)*myBcols
  ALLOCATE(ipiv(nelementsw))

  !CALL THE SCALAPACK ROUTINE
  !Solve the linear system A * X = B
  SELECT CASE (mk)
  CASE (SELECTED_REAL_KIND(6,37))
     CALL PSGESV(N,1,myA,1,1,desca,ipiv,myB,1,1,descb,info)
  CASE (SELECTED_REAL_KIND(15,307)) 
     CALL PDGESV(N,1,myA,1,1,desca,ipiv,myB,1,1,descb,info)
  END SELECT

  IF (info /= 0) THEN
      PRINT *, myrank,'Error -- info = ', info
      CALL MPI_Abort(MPI_COMM_WORLD,1,ierr)
  ENDIF 

  WRITE(fnameb,'(A)')'sol.dat'
  CALL PDLAWRITE(trim(fnameb),N,1,myB,1,1,DESCB,0,0,ipiv)
  !CALL MPI_File_open(MPI_COMM_WORLD,TRIM(fnameb),MPI_MODE_CREATE+ &
  !&    MPI_MODE_WRONLY,MPI_INFO_NULL,infile,ierr)
  !CALL MPI_File_write_all(infile,myB,nelementsb,mpik,mpistatus,ierr)
  !CALL MPI_File_close(infile,ierr)

  ! DEALLOCATE the local arrays
  DEALLOCATE(myA,ipiv,myB)

  ! End blacs for processors that are used
  CALL blacs_gridexit(ictxt)
  calctime = tock(clock)

  ! PRINT results
  IF (me == 0) THEN
    IF (info /= 0) THEN
         PRINT *, 'Error -- info = ', info
    ENDIF
    PRINT *,'Compute time = ', calctime
  ENDIF
  CALL MPI_Finalize(ierr)
CONTAINS
  SUBROUTINE tick(t)
      INTEGER, INTENT(  OUT) :: t
      CALL system_clock(t)
   END SUBROUTINE tick
  ! returns time in seconds from now to time described by t
  FUNCTION tock(t)
    INTEGER, INTENT(IN   ) :: t
    REAL(mk)               :: tock
    INTEGER :: now, clock_rate
    CALL system_clock(now,clock_rate)
    tock = REAL(now-t,mk)/REAL(clock_rate,mk)
  END FUNCTION tock
END PROGRAM array

矩阵A和B是标量的例子。用二进制格式重写。因此，它们可以被上面的程序使用。

 19.0  3.0  1.0  12.0  1.0  16.0  1.0  3.0  11.0 
-19.0  3.0  1.0  12.0  1.0  16.0  1.0  3.0  11.0
-19.0 -3.0  1.0  12.0  1.0  16.0  1.0  3.0  11.0 
-19.0 -3.0 -1.0  12.0  1.0  16.0  1.0  3.0  11.0 
-19.0 -3.0 -1.0 -12.0  1.0  16.0  1.0  3.0  11.0 
-19.0 -3.0 -1.0 -12.0 -1.0  16.0  1.0  3.0  11.0 
-19.0 -3.0 -1.0 -12.0 -1.0 -16.0  1.0  3.0  11.0 
-19.0 -3.0 -1.0 -12.0 -1.0 -16.0 -1.0  3.0  11.0
-19.0 -3.0 -1.0 -12.0 -1.0 -16.0 -1.0 -3.0  11.0

和B

0.0
0.0
1.0
0.0
0.0
0.0
0.0
0.0
0.0

例如，对于一个处理器：

mpirun -np 1 ./array 9 A.dat B.dat

           9           1
      0.000000000000000000D+00
     -0.166666666666666657D+00
      0.500000000000000000D+00
      0.000000000000000000D+00
      0.000000000000000000D+00
      0.000000000000000000D+00
      0.000000000000000000D+00
      0.000000000000000000D+00
      0.000000000000000000D+00

而对于2个处理器

           9           1
     -0.432625129877109577D-01
     -0.000000000000000000D+00
     -0.677331518039483299D-01
     -0.250000000000000000D+00
      0.250000000000000056D+00
      0.769230769230769273D-01
      0.117743386633788305D+00
      0.122377622377622383D+00
      0.157721108027438911D+00

对于6个处理器，我会收到一个错误

mpirun -np 6 ./a.out 9 A.dat B.dat

           0 Error -- info =            3
           1 Error -- info =            3
           2 Error -- info =            3
           3 Error -- info =            3
           4 Error -- info =            3
           5 Error -- info =            3

使用scalapack示例不会发生这种情况。

Answer 1

问题解决了。

首先，在打开文件后，我修复了进程在文件中的数据视图

CALL MPI_File_set_view(infile,0_MPI_OFFSET_KIND,mpik,darraya,"native",MPI_INFO_NULL,ierr)

其次，我修正了在为ScaLAPACK准备数组描述符时所犯的错误

CALL descinit(descb,N,1,Nb, 1,0,0,ictxt,max(1,myArows),info)

至

CALL descinit(descb,N,1,Nb,Nb,0,0,ictxt,max(1,myArows),info)

这些变化是我遗漏的问题的一部分。