DOUG 0.2

ElemMtxs_base.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 !-------------------------------------------------------
00025 module ElemMtxs_base
00026 
00027   use DOUG_utils
00028   use globals
00029   use parameters
00030   use Mesh_class
00031 
00032 #include<doug_config.h>
00033 
00034 #ifdef D_COMPLEX
00035 #define float complex
00036 #else
00037 #define float real
00038 #endif
00039 
00040   implicit none
00041 
00042   !--------------------------------------------------------------------
00047   integer(kind=1), parameter :: D_ELEM_INNER  = 0
00048   integer(kind=1), parameter :: D_ELEM_INTERF = 1
00049 
00050 
00051   !--------------------------------------------------------------------
00054   integer, parameter :: D_ELEMMTXS_IN_PACKET = 1024
00055 
00056 
00057   !--------------------------------------------------------------------
00061   type ElemMtxsChunk
00062      integer :: nell !< number of elements currently used in elem_idx and elem arrays
00063 
00065      float(kind=rk), dimension(:,:,:), pointer :: elem
00066      float(kind=rk), dimension(:,:),   pointer :: elemrhs
00067 
00069      integer, dimension(:), pointer :: lg_emap
00070   end type ElemMtxsChunk
00071 
00072 
00073   !--------------------------------------------------------------------
00077   type ElemMtxsPacket
00079      type(ElemMtxsChunk)           :: chunk
00080 
00081      logical                       :: request_in_progress 
00082 
00083      integer, dimension(3)         :: request 
00084 
00085      type(ElemMtxsPacket), pointer :: next
00086   end type ElemMtxsPacket
00087 
00088 contains
00089 
00090   !=============================================
00091   ! Chunk API
00092   !=============================================
00093 
00094   !---------------------------------------------
00097   subroutine ElemMtxsChunk_Init(chunk, nell, Msh)
00098     implicit none
00099 
00100     type(ElemMtxsChunk), intent(in out) :: chunk !< chunk to initialize
00101     integer,             intent(in)     :: nell !< number of elements to allocate in chunk
00102     type(Mesh),          intent(in)     :: Msh !< mesh corresponding to elements
00103 
00104     chunk%nell = 0
00105     allocate(chunk%elem(Msh%mfrelt, Msh%mfrelt, nell))
00106     allocate(chunk%elemrhs(Msh%mfrelt, nell))
00107     allocate(chunk%lg_emap(nell))
00108   end subroutine ElemMtxsChunk_Init
00109 
00110 
00111   !---------------------------------------------
00114   subroutine ElemMtxsChunk_Destroy(chunk)
00115     implicit none
00116 
00117     type(ElemMtxsChunk), intent(in out) :: chunk
00118 
00119     if (associated(chunk%lg_emap)) deallocate(chunk%lg_emap)
00120     if (associated(chunk%elemrhs)) deallocate(chunk%elemrhs)
00121     if (associated(chunk%elem)) deallocate(chunk%elem)
00122   end subroutine ElemMtxsChunk_Destroy
00123 
00124 
00125   !---------------------------------------------
00129   subroutine ElemMtxsChunk_recv(chunk, Msh, p)
00130     implicit none
00131 
00132     type(ElemMtxsChunk), intent(in out) :: chunk
00133     type(Mesh),          intent(in)     :: Msh !< mesh associated with element matrices
00134     integer,             intent(in), optional :: p
00135 
00136     integer :: ierr
00137     integer :: bufsize, source
00138     integer :: status(MPI_STATUS_SIZE)
00139 
00140     ! First get number of elements in chunk, then receive lg_emap, rhs, matrix
00141     if (present(p)) then
00142        source = p
00143     else
00144        source = MPI_ANY_SOURCE
00145     end if
00146     call MPI_PROBE(source, D_TAG_ELEMMTXS_ELEMIDXS, MPI_COMM_WORLD, status, ierr)
00147     call MPI_GET_COUNT(status, MPI_INTEGER, chunk%nell, ierr)
00148     call MPI_RECV(chunk%lg_emap, chunk%nell, MPI_INTEGER, &
00149          source, D_TAG_ELEMMTXS_ELEMIDXS, MPI_COMM_WORLD, status, ierr)
00150     bufsize = Msh%mfrelt*chunk%nell
00151     call MPI_RECV(chunk%elemrhs, bufsize, MPI_fkind, &
00152          source, D_TAG_ELEMMTXS_ELEMRHS, MPI_COMM_WORLD, status, ierr)
00153     bufsize = Msh%mfrelt*Msh%mfrelt*chunk%nell
00154     call MPI_RECV(chunk%elem, bufsize, MPI_fkind, &
00155          source, D_TAG_ELEMMTXS_ELEMS, MPI_COMM_WORLD, status, ierr)
00156   end subroutine ElemMtxsChunk_recv
00157 
00158 
00159   !---------------------------------------------
00162   subroutine ElemMtxsChunk_fillInDense(D, chunk, M)
00163     implicit none
00164 
00165     float(kind=rk), dimension(:,:), intent(in out) :: D !< dense matrix
00166     type(ElemMtxsChunk),                intent(in) :: chunk !< element matrices
00167     type(Mesh),                         intent(in) :: M !< mesh associated with element matrices
00168 
00169     integer :: k, i, j, g, il, jl
00170     
00171     do k = 1,chunk%nell ! cycle through local elements
00172        g = chunk%lg_emap(k) ! map index of element to global numbering
00173        
00174        do j = 1,M%nfrelt(g)
00175           do i = 1,M%nfrelt(g)
00176              il = M%gl_fmap(M%mhead(i,g))
00177              jl = M%gl_fmap(M%mhead(j,g))
00178              if (chunk%elem(i,j,k) /= 0.0_rk) then
00179                 D(il,jl) = D(il,jl) + chunk%elem(i,j,k) ! k <- local index
00180              end if
00181           end do
00182        end do
00183     end do
00184   end subroutine ElemMtxsChunk_fillInDense
00185 
00186 
00187   !---------------------------------------------
00190   subroutine ElemMtxsChunk_print(chunk)
00191     implicit none
00192 
00193     type(ElemMtxsChunk), intent(in) :: chunk 
00194 
00195     integer :: nell, el, i, j
00196 
00197     write(stream,*) chunk%nell ,' element matrices:'
00198     call flush(stream)
00199 
00200     nell = chunk%nell
00201     if (chunk%nell > 25) then
00202        nell = 25
00203        write(stream,'(a,i4,a)') '[',nell,' matrices will be printed]'
00204     end if
00205 
00206     do el = 1,nell
00207        write(stream,*) 'element matrix:',el
00208        do i = 1,size(chunk%elem, dim=1)
00209           do j = 1,size(chunk%elem, dim=2)
00210              write(stream, '(f7.4)', advance='no') chunk%elem(i,j,el)
00211              if (j /= size(chunk%elem, dim=2)) write(stream,'(a)', advance='no') ', '
00212           end do
00213           write(stream,*)
00214           call flush(stream)
00215        end do
00216     end do
00217   end subroutine ElemMtxsChunk_print
00218 
00219 
00220   !=============================================
00221   ! Packet API
00222   !=============================================
00223 
00224   !---------------------------------------------
00227   subroutine ElemMtxsPacket_Init(packet, Msh, nelems)
00228     implicit none
00229 
00230     type(ElemMtxsPacket), intent(out) :: packet
00231     type(Mesh),           intent(in)  :: Msh
00232     integer, optional,    intent(in)  :: nelems
00233 
00234     if (present(nelems)) then
00235        call ElemMtxsChunk_Init(packet%chunk, nelems, Msh) 
00236     else
00237        call ElemMtxsChunk_Init(packet%chunk, D_ELEMMTXS_IN_PACKET, Msh) 
00238     end if
00239     packet%request_in_progress = .false.
00240     packet%next => NULL()
00241   end subroutine ElemMtxsPacket_Init
00242 
00243   !---------------------------------------------
00246   subroutine ElemMtxsPacket_Destroy(packet)
00247     implicit none
00248 
00249     type(ElemMtxsPacket), intent(in out) :: packet
00250     
00251     call ElemMtxsChunk_Destroy(packet%chunk)
00252   end subroutine ElemMtxsPacket_Destroy
00253 
00254 
00255   !---------------------------------------------
00258   subroutine ElemMtxsPacket_send(packet, Msh, p)
00259     implicit none
00260 
00261     type(ElemMtxsPacket), intent(in out),target :: packet
00262     type(Mesh),           intent(in)     :: Msh
00263     integer,              intent(in)     :: p
00264 
00265     integer :: ierr
00266     integer :: bufsize
00267     type(ElemMtxsChunk), pointer :: chunk
00268 
00269     ! Send lg_emap, rhs and matrix
00270     packet%request_in_progress = .true.
00271     chunk => packet%chunk
00272     call MPI_ISEND(chunk%lg_emap, chunk%nell, MPI_INTEGER, &
00273          p, D_TAG_ELEMMTXS_ELEMIDXS, MPI_COMM_WORLD, packet%request(1), ierr)
00274     bufsize = Msh%mfrelt*chunk%nell
00275     call MPI_ISEND(chunk%elemrhs, bufsize, MPI_fkind, &
00276          p, D_TAG_ELEMMTXS_ELEMRHS, MPI_COMM_WORLD, packet%request(2), ierr)
00277     bufsize = Msh%mfrelt*Msh%mfrelt*chunk%nell
00278     call MPI_ISEND(chunk%elem, bufsize, MPI_fkind, &
00279          p, D_TAG_ELEMMTXS_ELEMS, MPI_COMM_WORLD, packet%request(3), ierr)
00280   end subroutine ElemMtxsPacket_send
00281 
00282 
00283   !---------------------------------------------
00286   function ElemMtxsPacket_sendInProgress(packet) result(b)
00287     implicit none
00288 
00289     type(ElemMtxsPacket), intent(in out) :: packet
00290   
00291     integer :: ierr
00292     logical :: b, flag
00293 
00294     b = .false.
00295     if (packet%request_in_progress) then
00296        call MPI_TESTALL(size(packet%request), packet%request, flag, MPI_STATUSES_IGNORE, ierr)
00297        if (.not.flag) then
00298           b = .true.
00299        else
00300           packet%request_in_progress = .false.
00301        end if
00302     end if
00303   end function ElemMtxsPacket_sendInProgress
00304 
00305 
00306   !---------------------------------------------
00309   subroutine ElemMtxsPacket_wait(packet)
00310     implicit none
00311 
00312     type(ElemMtxsPacket), intent(in out) :: packet
00313 
00314     integer :: ierr
00315 
00316     if (packet%request_in_progress) then
00317        call MPI_WAITALL(size(packet%request), packet%request, MPI_STATUSES_IGNORE, ierr)
00318        packet%request_in_progress = .false.
00319     end if
00320   end subroutine ElemMtxsPacket_wait
00321   
00322 
00323 end module ElemMtxs_base