高性能并行计算.doc

cannon.f* subroutine cannon( a, lda, b, ldb, c, ldc, m, n, k, rowcomm, & colcomm, w, iw ) implicit none include mpif.h integer lda, ldb, ldc, m, n, k, rowcomm, colcomm, iw(*) real a(lda, *), b(ldb, *), c(ldc, *), w(*)* integer lma, lka, lkb, lnb, lmc, lnc, ldw, ldw1* integer nr, nc, rid, cid, ierr, res, arect, brect, nrb integer root, north, south, sta(mpi_status_size), i* call mpi_comm_size( colcomm, nr, ierr ) call mpi_comm_rank( colcomm, rid, ierr ) call mpi_comm_size( rowcomm, nc, ierr ) call mpi_comm_rank( rowcomm, cid, ierr )* lma = m/nr res = mod( m, nr ) if ( rid .lt. res ) lma = lma + 1* lka = k/nc res = mod( k, nc ) if ( cid .lt. res ) lka = lka + 1* lkb = k/nr res = mod( k, nr ) if ( rid .lt. res ) lkb = lkb + 1* lnc = n/nc res = mod( n, nc ) if ( cid .lt. res ) lnc = lnc + 1 lmc = lma lnb = lnc ldw = lma + 1 call mpi_allgather( lkb, 1, mpi_integer, iw, 1, mpi_integer, & colcomm, ierr ) nrb = iw(1) ldw1 = ldb* if ( nr .ne. nc ) return* call mpirect( lda, lma, nrb, arect ) call mpi_type_commit( arect, ierr ) call mpirect( ldb, nrb, lnb, brect ) call mpi_type_commit( brect, ierr ) call wrapinita(a, lda, lma, lka, rid, cid, nr, nc) call wrapinitb(b, ldb, lkb, lnb, rid, cid, nr, nc) call zeroc( c, ldc, lmc, lnc )* north = mod( nr+rid-1, nr ) south = mod( rid+1, nr ) root = 0* do 100 i=0, nr-1 root = mod(rid + i, nr) call mcopy( a, lda, w, ldw, lma, lka ) call mpi_bcast( w, 1, arect, root, rowcomm, ierr ) k = root + 1 call sgemm(w, ldw, b, ldb, c, ldc, lma, iw(k), lnc)* c = c+a*b if ( i .lt. nr-1 ) then call mpi_sendrecv( b, 1, brect, north, 1, w, 1, brect, & south, 1, colcomm, sta, ierr ) k = mod( root + 1, nr ) +1* call mcopy( w, ldw1, b, ldb, iw(k), lnb ) endif 100 continue* call mpi_type_free( arect, ierr ) call mpi_type_free( brect, ierr )* return end*datamove.f program main implicit none include mpif.h* integer comm, np, iam integer ierr integer m, n, sta(mpi_status_size), front, next* call mpibegin( comm, np, iam ) print *, Hello world! on Proc. , iam front = mod( iam -1 + np, np ) next = mod( iam + 1, np )* case1 m = iam goto 20 if ( mod( iam, 2 ) .eq. 0 .and. iam .ne. np-1 ) then call mpi_recv( n, 1, mpi_integer, iam+1, 1, comm, sta, ierr ) else if ( mod( iam, 2) .ne. 0 .and. iam .ne. 0 ) then call mpi_send( m, 1, mpi_integer, iam-1, 1, comm, ierr ) endif if ( mod( iam, 2 ) .eq. 0 .and. iam .ne. 0 ) then call mpi_send( m, 1, mpi_integer, iam-1, 1, comm, ierr ) else if ( mod( iam, 2) .ne. 0 .and. iam .ne. np-1 ) then call mpi_recv( n, 1, mpi_integer, iam+1, 1, comm, sta, ierr ) endif 20 continue if ( iam .eq. 0 ) then front = mpi_proc_null else if ( iam .eq. np-1 ) then next = mpi_proc_null endif call mpi_sendrecv( m, 1, mpi_integer, front, 1, n, 1, & mpi_integer, next, 1, comm, sta, ierr ) if ( iam .ne. np-1 ) m = n print *, value m = , m, on proc. , iam call mpiend( ) endg2dmesh.f* Generate a 2-d mesh mpi environment* subroutine g2dmesh( comm, np, iam, p, q, rowcomm, colcomm, & rowid, colid ) include mpif.h integer comm, np, iam, p, q, rowcomm, colcomm, rowid, colid* row major manner to make the mapping from 1-d to 2-d* integer color, key, ierr key = iam if( p*q .gt. np .or. iam .ge. p*q ) then color = mpi_undefined call mpi_comm_split( comm, color, key, rowcomm, ierr ) call mpi_comm_split( comm, color, key, colcomm, ierr ) return endif* generate row communicator color = iam / q call mpi_comm_split( comm, color, key, rowcomm, ierr ) call mpi_comm_rank( rowcomm, colid, ierr )* color = mod( iam, q ) call mpi_comm_split( comm, color, key, colcomm, ierr ) call mpi_comm_rank( colcomm, rowid, ierr ) return endgroup.f* Group function testing* program grptst implicit none include mpif.h* integer comm, iam, np, ierr, grp, grp1, grp2 integer ranks(10), newcomm, m, root, newcom2* call mpibegin( comm, np, iam ) if ( np .lt. 10 ) go to 99 call mpi_comm_group( comm, grp, ierr ) ranks(1) = 1 ranks(2) = 4 ranks(3) = 7 m = iam root = 0 call mpi_group_incl( grp, 3, ranks, grp1, ierr ) call mpi_comm_create( comm, grp1, newcomm, ierr ) if ( newcomm .ne. mpi_comm_null ) then call mpi_bcast( m, 1, mpi_integer, root, newcomm, ierr ) endif print *, newcomm, in proc , iam, m= , m call mpi_group_free( grp1, ierr ) if ( newcomm .ne. mpi_comm_null ) then call mpi_comm_free( newcomm, ierr ) endif ranks(1) = 3 ranks(2) = 8 ranks(3) = 2 call mpi_group_range_incl( grp, 1, ranks, grp2, ierr ) call mpi_comm_create( comm, grp2, newcom2, ierr ) if ( newcom2 .ne. mpi_comm_null ) then call mpi_bcast( m, 1, mpi_integer, root, newcom2, ierr ) endif print *, newcomm=, newcom2, in proc , iam, m= , m call mpi_group_free( grp2, ierr ) if ( newcom2 .ne. mpi_comm_null ) then call mpi_comm_free( newcom2, ierr ) endif* call mpi_group_free( grp, ierr ) call mpi_comm_free( comm, ierr )* 99 call mpiend( ) endlower.f* define a lower triangle matrix* subroutine mpilower( lda, m, lower, work ) include mpif.h integer lda, m, lower, work(*)* integer ct, ierr, disps, blks, i* ct = m blks = 1 disps = m+1 do 20 i=0, m-1 work( blks+i ) = m - i work(disps + i) = i*lda + i 20 continue call mpi_type_indexed( ct, work(blks), work(disps), mpi_real, & lower, ierr )* return endlowerm.f* define a lower triangle matrix for a special purpose* subroutine mpilowerm( lda, m, lowerm, locub, work ) include mpif.h integer lda, m, lowerm, locub, work(*)* integer ct, ierr, lower, disps, blks, i* ct = m blks = 1 disps = m+1 do 20 i=0, m-1 work( blks+i ) = m - i work(disps + i) = i*lda + i 20 continue call mpi_type_indexed( ct, work(blks), work(disps), mpi_real, & lower, ierr )* work(1) = 1 work(2) = 1* work(3) = 0 work(4) = locub* work(5)= lower work(6) = mpi_ub* call mpi_type_struct( 2, work, work(3), work(5), lowerm, ierr )* return endmain.f program main implicit none include mpif.h* integer comm, np, iam integer ierr, lda parameter( lda = 50 ) integer m, n, k, sta(mpi_status_size), front, next real a(lda, lda), b(lda,lda), c(lda, lda), w(lda*lda) integer lower, i, j, iw(lda), loc integer rowcomm, colcomm, rowid, colid, p, q* call mpibegin( comm, np, iam ) print *, Hello world! on Proc. , iam front = mod( iam -1 + np, np ) next = mod( iam + 1, np )* p = 2 q = 2 m = 70 k = 70 n = 65 loc = iam call g2dmesh(comm, np, iam, p, q, rowcomm, colcomm, rowid, colid) if(rowcomm.ne.mpi_comm_null.and.colcomm.ne.mpi_comm_null) then write(*, *) (, rowid, colid, ), in proc , iam call cannon(a, lda, b, lda, c, lda, m, n, k, rowcomm, colcomm, & w, iw) write(*, *) c(1,1), c(2,1), in proc , iam, rowid, colid else write( *, * ) rowcomm= , rowcomm, in proc , iam endif call mpiend( ) endmcopy.f subroutine mcopy(a, lda, t, ldt, ma, ka ) integer lda, ma, ka, ldt real a(lda, *), t(ldt, *) integer i, j do 10 j=1, ka do 10 i=1, ma t(i, j) = a(i, j)10 continue return endmpibegin.f* This file is created on March 29, 2010* For entering the mpi environment* subroutine mpibegin( comm, np, iam ) include mpif.h* integer comm, np, iam integer ierr* call mpi_init( ierr ) call mpi_comm_dup( mpi_comm_world, comm, ierr ) call mpi_comm_size( comm, np, ierr ) call mpi_comm_rank( comm, iam, ierr )* return endmpiend.f* This file is created on March 29, 2010* For exiting the mpi environment* subroutine mpiend( ) include mpif.h* integer ierr* call mpi_finalize( ierr )* return endmpipi.f program computing_pi*The header file for using MPI parallel environment,* which must be included for all mpi programs. include mpif.h*Variables declaration integer iam, np, comm, ierr integer n, i, num, is, ie real*8 pi, h, eps, xi, s*Enroll in MPI environment and get the MPI parameters call mpi_init(ierr) call mpi_comm_dup(mpi_comm_world, comm, ierr) call mpi_comm_rank(comm, iam, ierr) call mpi_comm_size(comm, np, ierr)*Read the number of digits you want for value of Pi. if(iam .eq. 0) then write(*, *) Number of digits(1-16)= read(*, *) num endif call mpi_bcast(num,1,mpi_integer,0,comm,ierr) eps = 1 do 10 i=1, num eps = eps * 0.110 continue n = sqrt(4.0/(3.0*eps) h = 1.0/n num = n/np if(iam .eq. 0) then s = 3.0 xi = 0 is = 0 ie = num elseif(iam .eq. np-1) then s = 0.0 is = iam*num ie = n - 1 xi = is * h else s = 0.0 is = iam*num ie = is + num xi = is * h endif if(np .eq. 1) ie = ie - 1 do 20 i=is+1, ie xi = xi + h s = s + 4.0/(1.0+xi*xi)20 continue call mpi_reduce(s, pi, 1, mpi_double_precision, & mpi_sum, 0, comm, ierr) if(iam .eq. 0) then pi = h*pi write(*, 99) pi endif call mpi_finalize(ierr)99 format(The pi= , f16.13) endmpisum.f* This program is to deal with the summation of a series of numbers* in each processes. a = s + s +* subroutine mpisum( comm, np, iam, s, a ) implicit none include mpif.h integer comm, np, iam real s, a integer miam, lens, i, n, ierr, sta(mpi_status_size) a = s n = alog( np-0.1 ) / alog( 2.0 ) + 1* print *, n, in mpisum lens = 1 do 100 i = 1, n miam = mod( iam, 2*lens ) if ( miam .eq. 0 .and. iam+lens .lt. np ) then call mpi_recv( s, 1, mpi_real, iam+lens, i, comm, sta, ierr ) a = a + s else if ( miam .eq. lens .and. iam .ne. 0 ) then call mpi_send( a, 1, mpi_real, iam-lens, i, comm, ierr ) endif lens = 2*lens 100 continue* return endrect.f* define a rectangle matrix* subroutine mpirect( lda, m, n, rect ) include mpif.h integer lda, m, n, rect* integer ct, ierr, stride, blk* ct = n blk = m stride = lda call mpi_type_vector( ct, blk, stride, mpi_real, rect, ierr )* return endsgemm.f subroutine sgemm( t, ldt, b, ldb, c, ldc, ma, & kt, nb) integer ldt, ldb, ldc, ma, kt, nb