teshsuite/smpi/mpich-test/topol/twod.f

   1 c**********************************************************************
   2 c   twod.f - a solution to the Poisson problem by using Jacobi
   3 c   interation on a 2-d decomposition
   4 c
   5 c   .... the rest of this is from pi3.f to show the style ...
   6 c
   7 c   Each node:
   8 c    1) receives the number of rectangles used in the approximation.
   9 c    2) calculates the areas of it's rectangles.
  10 c    3) Synchronizes for a global summation.
  11 c   Node 0 prints the result.
  12 c
  13 c  Variables:
  14 c
  15 c    pi  the calculated result
  16 c    n   number of points of integration.
  17 c    x           midpoint of each rectangle's interval
  18 c    f           function to integrate
  19 c    sum,pi      area of rectangles
  20 c    tmp         temporary scratch space for global summation
  21 c    i           do loop index
  22 c
  23 c     This code is included (without the prints) because one version of
  24 c     MPICH SEGV'ed (probably because of errors in handling send/recv of
  25 c     MPI_PROC_NULL source/destination).
  26 c
  27 c****************************************************************************
  28       program main
  29       include "mpif.h"
  30       integer maxn
  31       parameter (maxn = 128)
  32       double precision  a(maxn,maxn), b(maxn,maxn), f(maxn,maxn)
  33       integer nx, ny
  34       integer myid, numprocs, it, rc, comm2d, ierr, stride
  35       integer nbrleft, nbrright, nbrtop, nbrbottom
  36       integer sx, ex, sy, ey
  37       integer dims(2)
  38       logical periods(2)
  39       double precision diff2d, diffnorm, dwork
  40       double precision t1, t2
  41       external diff2d
  42       data periods/2*.false./
  43
  44       call MPI_INIT( ierr )
  45       call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
  46       call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
  47 c      print *, "Process ", myid, " of ", numprocs, " is alive"
  48       if (myid .eq. 0) then
  49 c
  50 c         Get the size of the problem
  51 c
  52 c          print *, 'Enter nx'
  53 c          read *, nx
  54            nx = 10
  55       endif
  56 c      print *, 'About to do bcast on ', myid
  57       call MPI_BCAST(nx,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
  58       ny   = nx
  59 c
  60 c Get a new communicator for a decomposition of the domain.  Let MPI
  61 c find a "good" decomposition
  62 c
  63       dims(1) = 0
  64       dims(2) = 0
  65       call MPI_DIMS_CREATE( numprocs, 2, dims, ierr )
  66       call MPI_CART_CREATE( MPI_COMM_WORLD, 2, dims, periods, .true.,
  67      *                    comm2d, ierr )
  68 c
  69 c Get my position in this communicator
  70 c
  71       call MPI_COMM_RANK( comm2d, myid, ierr )
  72 c      print *, "Process ", myid, " of ", numprocs, " is alive"
  73 c
  74 c My neighbors are now +/- 1 with my rank.  Handle the case of the
  75 c boundaries by using MPI_PROCNULL.
  76       call fnd2dnbrs( comm2d, nbrleft, nbrright, nbrtop, nbrbottom )
  77 c      print *, "Process ", myid, ":",
  78 c     *     nbrleft, nbrright, nbrtop, nbrbottom
  79 c
  80 c Compute the decomposition
  81 c
  82       call fnd2ddecomp( comm2d, nx, sx, ex, sy, ey )
  83 c      print *, "Process ", myid, ":", sx, ex, sy, ey
  84 c
  85 c Create a new, "strided" datatype for the exchange in the "non-contiguous"
  86 c direction
  87 c
  88       call mpi_Type_vector( ey-sy+1, 1, ex-sx+3,
  89      $                      MPI_DOUBLE_PRECISION, stride, ierr )
  90       call mpi_Type_commit( stride, ierr )
  91 c
  92 c
  93 c Initialize the right-hand-side (f) and the initial solution guess (a)
  94 c
  95       call twodinit( a, b, f, nx, sx, ex, sy, ey )
  96 c
  97 c Actually do the computation.  Note the use of a collective operation to
  98 c check for convergence, and a do-loop to bound the number of iterations.
  99 c
 100       call MPI_BARRIER( MPI_COMM_WORLD, ierr )
 101       t1 = MPI_WTIME()
 102       do 10 it=1, 100
 103         call exchng2( a, b, sx, ex, sy, ey, comm2d, stride,
 104      $                nbrleft, nbrright, nbrtop, nbrbottom )
 105         call sweep2d( b, f, nx, sx, ex, sy, ey, a )
 106         call exchng2( b, a, sx, ex, sy, ey, comm2d, stride,
 107      $                nbrleft, nbrright, nbrtop, nbrbottom )
 108         call sweep2d( a, f, nx, sx, ex, sy, ey, b )
 109         dwork = diff2d( a, b, nx, sx, ex, sy, ey )
 110         call MPI_Allreduce( dwork, diffnorm, 1, MPI_DOUBLE_PRECISION,
 111      $                     MPI_SUM, comm2d, ierr )
 112         if (diffnorm .lt. 1.0e-5) goto 20
 113 c        if (myid .eq. 0) print *, 2*it, ' Difference is ', diffnorm
 114 10     continue
 115       if (myid .eq. 0) print *, 'Failed to converge'
 116 20    continue
 117       t2 = MPI_WTIME()
 118 c      if (myid .eq. 0) then
 119 c         print *, 'Converged after ', 2*it, ' Iterations in ', t2 - t1,
 120 c     $        ' secs '
 121 c      endif
 122 c
 123 c
 124       call MPI_Type_free( stride, ierr )
 125       call MPI_Comm_free( comm2d, ierr )
 126       if (myid .eq. 0) then
 127          print *, ' No Errors'
 128       endif
 129       call MPI_FINALIZE(rc)
 130       end
 131 c
 132 c Perform a Jacobi sweep for a 2-d decomposition
 133 c
 134       subroutine sweep2d( a, f, n, sx, ex, sy, ey, b )
 135       integer n, sx, ex, sy, ey
 136       double precision a(sx-1:ex+1, sy-1:ey+1), f(sx-1:ex+1, sy-1:ey+1),
 137      +                 b(sx-1:ex+1, sy-1:ey+1)
 138 c
 139       integer i, j
 140       double precision h
 141 c
 142       h = 1.0d0 / dble(n+1)
 143       do 10 i=sx, ex
 144          do 10 j=sy, ey
 145             b(i,j) = 0.25 * (a(i-1,j)+a(i,j+1)+a(i,j-1)+a(i+1,j)) -
 146      +               h * h * f(i,j)
 147  10   continue
 148       return
 149       end
 150
 151        subroutine exchng2( a, b, sx, ex, sy, ey,
 152      $                     comm2d, stride,
 153      $                     nbrleft, nbrright, nbrtop, nbrbottom  )
 154        include "mpif.h"
 155        integer sx, ex, sy, ey, stride
 156        double precision a(sx-1:ex+1, sy-1:ey+1),
 157      $                  b(sx-1:ex+1, sy-1:ey+1)
 158        integer nbrleft, nbrright, nbrtop, nbrbottom, comm2d
 159        integer status(MPI_STATUS_SIZE), ierr, nx
 160 c
 161        nx = ex - sx + 1
 162 c  These are just like the 1-d versions, except for less data
 163         call MPI_SENDRECV( b(ex,sy),  nx, MPI_DOUBLE_PRECISION,
 164      $                    nbrtop, 0,
 165      $                    a(sx-1,sy), nx, MPI_DOUBLE_PRECISION,
 166      $                    nbrbottom, 0, comm2d, status, ierr )
 167         call MPI_SENDRECV( b(sx,sy),  nx, MPI_DOUBLE_PRECISION,
 168      $                    nbrbottom, 1,
 169      $                    a(ex+1,sy), nx, MPI_DOUBLE_PRECISION,
 170      $                    nbrtop, 1, comm2d, status, ierr )
 171 c
 172 c This uses the "strided" datatype
 173         call MPI_SENDRECV( b(sx,ey),  1, stride, nbrright,  0,
 174      $                     a(sx,sy-1), 1, stride, nbrleft,  0,
 175      $                     comm2d, status, ierr )
 176         call MPI_SENDRECV( b(sx,sy),  1, stride, nbrleft,   1,
 177      $                     a(sx,ey+1), 1, stride, nbrright, 1,
 178      $                     comm2d, status, ierr )
 179         return
 180         end
 181
 182 c
 183 c  The rest of the 2-d program
 184 c
 185       double precision function diff2d( a, b, nx, sx, ex, sy, ey )
 186       integer nx, sx, ex, sy, ey
 187       double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1)
 188 c
 189       double precision sum
 190       integer i, j
 191 c
 192       sum = 0.0d0
 193       do 10 j=sy,ey
 194          do 10 i=sx,ex
 195             sum = sum + (a(i,j) - b(i,j)) ** 2
 196  10      continue
 197 c
 198       diff2d = sum
 199       return
 200       end
 201       subroutine twodinit( a, b, f, nx, sx, ex, sy, ey )
 202       integer nx, sx, ex, sy, ey
 203       double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1),
 204      &                 f(sx-1:ex+1, sy-1:ey+1)
 205 c
 206       integer i, j
 207 c
 208       do 10 j=sy-1,ey+1
 209          do 10 i=sx-1,ex+1
 210             a(i,j) = 0.0d0
 211             b(i,j) = 0.0d0
 212             f(i,j) = 0.0d0
 213  10      continue
 214 c
 215 c    Handle boundary conditions
 216 c
 217       if (sx .eq. 1) then
 218          do 20 j=sy,ey
 219             a(0,j) = 1.0d0
 220             b(0,j) = 1.0d0
 221  20      continue
 222       endif
 223       if (ex .eq. nx) then
 224          do 21 j=sy,ey
 225             a(nx+1,j) = 0.0d0
 226             b(nx+1,j) = 0.0d0
 227  21      continue
 228       endif
 229       if (sy .eq. 1) then
 230          do 30 i=sx,ex
 231             a(i,0) = 1.0d0
 232             b(i,0) = 1.0d0
 233  30      continue
 234       endif
 235 c
 236       return
 237       end
 238
 239 c
 240 c  This file contains a routine for producing a decomposition of a 1-d array
 241 c  when given a number of processors.  It may be used in "direct" product
 242 c  decomposition.  The values returned assume a "global" domain in [1:n]
 243 c
 244       subroutine MPE_DECOMP1D( n, numprocs, myid, s, e )
 245       integer n, numprocs, myid, s, e
 246       integer nlocal
 247       integer deficit
 248 c
 249       nlocal  = n / numprocs
 250       s       = myid * nlocal + 1
 251       deficit = mod(n,numprocs)
 252       s       = s + min(myid,deficit)
 253       if (myid .lt. deficit) then
 254           nlocal = nlocal + 1
 255       endif
 256       e = s + nlocal - 1
 257       if (e .gt. n .or. myid .eq. numprocs-1) e = n
 258       return
 259       end
 260 c
 261 c This routine show how to determine the neighbors in a 2-d decomposition of
 262 c the domain. This assumes that MPI_Cart_create has already been called
 263 c
 264       subroutine fnd2dnbrs( comm2d,
 265      $                      nbrleft, nbrright, nbrtop, nbrbottom )
 266       integer comm2d, nbrleft, nbrright, nbrtop, nbrbottom
 267 c
 268       integer ierr
 269 c
 270       call MPI_Cart_shift( comm2d, 0,  1, nbrleft,   nbrright, ierr )
 271       call MPI_Cart_shift( comm2d, 1,  1, nbrbottom, nbrtop,   ierr )
 272 c
 273       return
 274       end
 275 c
 276 c Note: THIS IS A TEST PROGRAM.  THE ACTUAL VALUES MOVED ARE NOT
 277 c CORRECT FOR A POISSON SOLVER.
 278 c
 279       subroutine fnd2ddecomp( comm2d, n, sx, ex, sy, ey )
 280       integer comm2d
 281       integer n, sx, ex, sy, ey
 282       integer dims(2), coords(2), ierr
 283       logical periods(2)
 284 c
 285       call MPI_Cart_get( comm2d, 2, dims, periods, coords, ierr )
 286
 287       call MPE_DECOMP1D( n, dims(1), coords(1), sx, ex )
 288       call MPE_DECOMP1D( n, dims(2), coords(2), sy, ey )
 289 c
 290       return
 291       end