DOUG 0.2

SpMtx_distribution.F90

Go to the documentation of this file.
00001 module SpMtx_distribution_mod
00002   use SpMtx_class
00003   use Mesh_class
00004   use globals
00005   use SpMtx_util
00006   use SpMtx_arrangement
00007 
00008   implicit none
00009 
00010 #include<doug_config.h>
00011 
00012 ! "on-the-fly" real/complex picking
00013 #ifdef D_COMPLEX
00014 #define float complex
00015 #else
00016 #define float real
00017 #endif
00018 
00019 contains
00020 
00022   function SpMtx_exchange(A,sends,M,lgi,lgj) result(R)
00023     type(SpMtx), intent(in) :: A !< matrix, which elements are to be sent
00024     type(indlist), intent(in) :: sends(:) !< indices of matrix elements to send to each neighbour
00025     type(Mesh), intent(in) :: M !< mesh that contains neighbour info
00026     integer, intent(in) :: lgi(:) !< local to global for A row indices
00027     integer, intent(in) :: lgj(:) !< local to global for A column indices
00028     type(SpMtx) :: R
00029 
00030     type(indlist), allocatable :: recvs(:)
00031     
00032     integer, allocatable :: outreqs(:)
00033     type Buffer
00034        character, pointer :: data(:)       
00035     end type Buffer
00036     type(Buffer), allocatable :: outbuffers(:)
00037     integer :: bufsize, bufpos
00038     character, allocatable :: inbuffer(:)
00039     integer, allocatable :: outstatuses(:,:)
00040 
00041     integer :: ninds, nnz
00042     integer :: i, ierr, ngh, status(MPI_STATUS_SIZE)
00043 
00044     allocate(recvs(size(sends)))
00045 
00046     ! gather values
00047     ! gather number of matrix elements each node has
00048     do i=1,M%nnghbrs
00049       ngh = M%nghbrs(i)
00050       call MPI_Send(sends(i)%ninds, 1, MPI_INTEGER, ngh, &
00051            TAG_CREATE_PROLONG, MPI_COMM_WORLD, ierr)
00052     end do
00053     do i=1,M%nnghbrs
00054       ngh = M%nghbrs(i)
00055       call MPI_Recv(recvs(i)%ninds, 1, MPI_INTEGER, ngh, &
00056            TAG_CREATE_PROLONG, MPI_COMM_WORLD, status, ierr)
00057     end do
00058 
00059     ! gather matrix elements
00060     ! non-blockingly send
00061     allocate(outbuffers(M%nnghbrs))
00062     allocate(outreqs(M%nnghbrs))
00063     do i=1,M%nnghbrs
00064       ! prepare and fill buffer
00065       bufsize = calcBufferSize(sends(i)%ninds)
00066       allocate(outbuffers(i)%data(bufsize))
00067       bufpos = 0
00068       call MPI_Pack(lgi(A%indi(sends(i)%inds)), sends(i)%ninds, MPI_INTEGER,&
00069            outbuffers(i)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
00070       if (ierr/=0) call DOUG_abort("MPI Pack of matrix elements failed")
00071       call MPI_Pack(lgj(A%indj(sends(i)%inds)), sends(i)%ninds, MPI_INTEGER,&
00072            outbuffers(i)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
00073       if (ierr/=0) call DOUG_abort("MPI Pack of matrix elements failed")
00074       call MPI_Pack(A%val(sends(i)%inds), sends(i)%ninds, MPI_fkind,&
00075            outbuffers(i)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
00076       if (ierr/=0) call DOUG_abort("MPI Pack of matrix elements failed")
00077 
00078       ! start sending values
00079       call MPI_ISend(outbuffers(i)%data, bufpos, MPI_CHARACTER, M%nghbrs(i),&
00080            TAG_CREATE_PROLONG, MPI_COMM_WORLD, outreqs(i), ierr)
00081     end do
00082 
00083     ! receive values
00084     allocate(inbuffer(calcBufferSize(maxval(recvs%ninds))))
00085     nnz = sum(recvs%ninds)
00086     R = SpMtx_newInit(nnz)
00087     nnz = 0
00088     do i=1,M%nnghbrs
00089       ! prepare and fill buffer
00090       bufsize = calcBufferSize(recvs(i)%ninds)
00091       call MPI_Recv(inbuffer, bufsize, MPI_CHARACTER, M%nghbrs(i),&
00092            TAG_CREATE_PROLONG, MPI_COMM_WORLD, status, ierr)
00093 
00094       ! do not even try to unpack if empty
00095       ninds = recvs(i)%ninds
00096       if (ninds==0) cycle
00097 
00098       bufpos = 0
00099       call MPI_Unpack(inbuffer, bufsize, bufpos, R%indi(nnz+1), ninds,&
00100            MPI_INTEGER, MPI_COMM_WORLD, ierr)
00101       if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
00102       call MPI_Unpack(inbuffer, bufsize, bufpos, R%indj(nnz+1), ninds,&
00103            MPI_INTEGER, MPI_COMM_WORLD, ierr)
00104       if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
00105       call MPI_Unpack(inbuffer, bufsize, bufpos, R%val(nnz+1), ninds,&
00106            MPI_fkind, MPI_COMM_WORLD, ierr)
00107       if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
00108       nnz = nnz+ninds
00109     end do
00110 
00111     ! wait for all data to be sent
00112     allocate(outstatuses(MPI_STATUS_SIZE, M%nnghbrs))
00113     call MPI_Waitall(M%nnghbrs, outreqs, outstatuses, ierr)
00114 
00115     ! deallocate
00116     do i=1,M%nnghbrs
00117       deallocate(outbuffers(i)%data)
00118     end do
00119     deallocate(outbuffers)
00120     deallocate(outreqs)
00121     deallocate(inbuffer)
00122     deallocate(recvs)
00123 
00124   contains
00126     function calcBufferSize(ninds) result(bufferSize)
00127       integer, intent(in) :: ninds
00128       integer :: bufferSize
00129       integer :: size
00130       call MPI_Pack_size(ninds, MPI_INTEGER, MPI_COMM_WORLD, size, ierr)
00131       bufferSize = 2*size
00132       call MPI_Pack_size(ninds, MPI_fkind, MPI_COMM_WORLD, size, ierr)
00133       bufferSize = bufferSize+size
00134     end function calcBufferSize
00135 
00136   end function SpMtx_exchange
00137 
00139   subroutine SpMtx_localize(A,A_ghost,b,M)
00140     type(SpMtx),intent(inout)           :: A,A_ghost
00141     float(kind=rk),dimension(:),pointer :: b
00142     type(Mesh)                          :: M
00143 
00144     integer,dimension(:),pointer       :: clrorder,clrstarts
00145     integer, dimension(:), allocatable :: ccount !count colors
00146     float(kind=rk),dimension(:),pointer :: b_tmp
00147     integer :: ol,i,k,n
00148 
00149     ol=max(sctls%overlap,sctls%smoothers)    
00150 
00151     !========= count color elements ============
00152     if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
00153       n=A%nrows
00154     elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
00155       n=A%ncols
00156     else
00157       call DOUG_abort('[SpMtx_DistributeAssembled] : matrix not arranged')
00158     endif
00159 
00160     allocate(ccount(numprocs))
00161     ccount=0
00162     do i=1,n
00163       ccount(M%eptnmap(i))=ccount(M%eptnmap(i))+1
00164     enddo
00165     allocate(clrstarts(numprocs+1))
00166     clrstarts(1)=1
00167     do i=1,numprocs
00168       clrstarts(i+1)=clrstarts(i)+ccount(i)
00169     end do
00170     allocate(clrorder(n))
00171     ccount(1:numprocs)=clrstarts(1:numprocs)
00172     do i=1,n
00173       clrorder(ccount(M%eptnmap(i)))=i
00174       ccount(M%eptnmap(i))=ccount(M%eptnmap(i))+1
00175     enddo
00176     if (sctls%verbose>3.and.A%nrows<200) then 
00177       do i=1,numprocs                                                     !
00178         write(stream,*)'partition ',i,' is in:', &                        !
00179           clrorder(clrstarts(i):clrstarts(i+1)-1)                     !
00180       enddo                                                               !
00181     endif
00182     deallocate(ccount)
00183 
00184     !-------------------------------------------------------------------+
00185     if (sctls%verbose>3.and.A%nrows<200) then 
00186       write(stream,*)'A after arrange:'
00187       call SpMtx_printRaw(A)
00188     endif
00189     call SpMtx_build_ghost(myrank+1,ol,&
00190                              A,A_ghost,M,clrorder,clrstarts) 
00191     if (sctls%verbose>3.and.A%nrows<300) then 
00192       write(stream,*)'A interf(1,1):'
00193       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,1),endnz=A%mtx_bbe(1,1))
00194       write(stream,*)'A interf(1,2):'
00195       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,2),endnz=A%mtx_bbe(1,2))
00196       write(stream,*)'A interf(2,1):'
00197       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,1),endnz=A%mtx_bbe(2,1))
00198       write(stream,*)'A inner:'
00199       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,2),endnz=A%mtx_bbe(2,2))
00200       if (ol>0) then
00201         write(stream,*)'A ghost:'
00202         call SpMtx_printRaw(A_ghost)
00203       endif
00204       if (A%nnz>A%mtx_bbe(2,2)) then
00205         write(stream,*)'A additional in case of ol==0:'
00206         call SpMtx_printRaw(A=A,startnz=A%mtx_bbe(2,2)+1,endnz=A%ol0nnz)
00207       endif
00208     endif
00209     ! print neighbours
00210     if (sctls%verbose>=1) then 
00211        write(stream,"(A,I0)") "N neighbours: ", M%nnghbrs
00212        do i=1,M%nnghbrs
00213           !write(stream,"(A,I0,A,I0,A,I0)") "neighbour ", i, ": ", M%nghbrs(i), ", overlap: ", M%ol_solve(i)%ninds
00214           write(stream,"(I0,A,I0,A)",advance="no") M%nghbrs(i)+1, " ", M%ol_inner(i)%ninds+M%ol_outer(i)%ninds, " "
00215        end do
00216        write(stream,*) ""
00217     end if
00218 
00219     ! Localise A:
00220     if (ol<=0) then
00221       M%ninonol=M%ntobsent
00222       M%indepoutol=M%ninner
00223     endif
00224     call SpMtx_Build_lggl(A,A_ghost,M)
00225     if (sctls%verbose>3) then 
00226       write(stream,*)'tobsent:',M%lg_fmap(1:M%ntobsent)
00227       write(stream,*)'...nintol:',M%lg_fmap(M%ntobsent+1:M%ninonol)
00228       write(stream,*)'...nninner:',M%lg_fmap(M%ninonol+1:M%ninner)
00229       write(stream,*)'...indepoutol:',M%lg_fmap(M%ninner+1:M%indepoutol)
00230       write(stream,*)'...ghost-freds:',M%lg_fmap(M%indepoutol+1:M%nlf)
00231     endif
00232     ! Rebuild RHS vector to correspond to local freedoms
00233     allocate(b_tmp(M%nlf))
00234     do i=1,M%nlf
00235       b_tmp(i)=b(M%lg_fmap(i))
00236     end do
00237     deallocate(b)
00238     b=>b_tmp
00239     ! Localise matrices and communication arrays
00240     do k=1,M%nnghbrs
00241       M%ax_recvidx(k)%inds=M%gl_fmap(M%ax_recvidx(k)%inds)
00242       M%ax_sendidx(k)%inds=M%gl_fmap(M%ax_sendidx(k)%inds)
00243     enddo
00244     do k=1,M%nnghbrs
00245       M%ol_inner(k)%inds=M%gl_fmap(M%ol_inner(k)%inds)
00246       M%ol_outer(k)%inds=M%gl_fmap(M%ol_outer(k)%inds)
00247       M%ol_solve(k)%inds=M%gl_fmap(M%ol_solve(k)%inds)
00248     enddo
00249     do i=1,A%ol0nnz
00250       A%indi(i)=M%gl_fmap(A%indi(i))
00251       A%indj(i)=M%gl_fmap(A%indj(i))
00252     enddo
00253     A%nrows=max(0, maxval(A%indi(1:A%nnz)))
00254     A%ncols=max(0, maxval(A%indj))
00255     A%arrange_type=D_SpMTX_ARRNG_NO
00256     if(associated(A%m_bound)) deallocate(A%m_bound) ! without this A_tmp got wrong size of M_bound in pcg()
00257     if(associated(A%strong)) deallocate(A%strong)
00258     if (ol>0) then
00259       do i=1,A_ghost%nnz
00260         A_ghost%indi(i)=M%gl_fmap(A_ghost%indi(i))
00261         A_ghost%indj(i)=M%gl_fmap(A_ghost%indj(i))
00262       enddo
00263       A_ghost%nrows=max(0, maxval(A_ghost%indi))
00264       A_ghost%ncols=max(0, maxval(A_ghost%indj))
00265       call SpMtx_arrange(A_ghost,D_SpMtx_ARRNG_ROWS,sort=.true.)
00266     endif
00267     if (sctls%verbose>3.and.A%nrows<200) then 
00268       write(stream,*)'Localised A interf(1,1):'
00269       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,1),endnz=A%mtx_bbe(1,1))
00270       write(stream,*)'Localised A interf(1,2):'
00271       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,2),endnz=A%mtx_bbe(1,2))
00272       write(stream,*)'Localised A interf(2,1):'
00273       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,1),endnz=A%mtx_bbe(2,1))
00274       write(stream,*)'Localised A inner:'
00275       call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,2),endnz=A%mtx_bbe(2,2))
00276       if (ol>0) then
00277         write(stream,*)'Localised A ghost:'
00278         call SpMtx_printRaw(A_ghost)
00279       endif
00280       if (A%nnz>A%mtx_bbe(2,2)) then
00281         write(stream,*)'localised A additional in case of ol==0:'
00282         call SpMtx_printRaw(A=A,startnz=A%mtx_bbe(2,2)+1,endnz=A%ol0nnz)
00283       endif
00284       write(stream,*)'gl_fmap:',M%gl_fmap
00285       write(stream,*)'gl_fmap(lg_fmap):',M%gl_fmap(M%lg_fmap)
00286       write(stream,*)'lg_fmap:',M%lg_fmap
00287       !call MPI_BARRIER(MPI_COMM_WORLD,ierr)
00288       !call DOUG_abort('testing nodal graph partitioning',0)
00289     endif
00290   end subroutine SpMtx_localize
00291 
00292 end module SpMtx_distribution_mod