DOUG 0.2

Distribution.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 
00023 module Distribution_mod
00024   use Distribution_base_mod
00025   use globals
00026   use DOUG_utils
00027   use SpMtx_operation
00028 
00029   use Distribution_elem_mod
00030   use Distribution_assm_mod
00031   use Distribution_struct_mod
00032   
00033   implicit none
00034 
00035 #include<doug_config.h>
00036 
00037 #ifdef D_COMPLEX
00038 #define float complex
00039 #else
00040 #define float real
00041 #endif
00042 
00043   private
00044   public :: Distribution_NewInit, Distribution_pmvm, Distribution_addoverlap
00045 
00046 contains
00047   !----------------------------------------------------------------
00050   function Distribution_NewInit(input_type, nparts, part_opts) result(D)
00051     implicit none
00052 
00053     integer,        intent(in)     :: input_type !< Input Type
00054     type(Distribution) :: D
00055     ! Partitioning
00056     integer, intent(in) :: nparts !< number of parts to partition a mesh
00057     integer, dimension(6), intent(in) :: part_opts !< partition options (see METIS manual)
00058 
00059     D = Distribution_New()
00060 
00061     select case (input_type)
00062     case (DISTRIBUTION_TYPE_ELEMENTAL)
00063        ! ELEMENTAL
00064        call parallelAssembleFromElemInput(D%mesh,D%A,D%rhs,nparts,part_opts,D%A_ghost)
00065     case (DISTRIBUTION_TYPE_ASSEMBLED)
00066        ! ASSEMBLED
00067        call parallelDistributeAssembledInput(D%mesh,D%A,D%rhs,D%A_ghost)
00068     case (DISTRIBUTION_TYPE_STRUCTURED)
00069        ! GENERATED
00070       if (sctls%grid_size<1) sctls%grid_size=100 
00071       D = Distribution_struct_NewInit(sctls%grid_size,max(sctls%overlap,sctls%smoothers))
00072     case default
00073        call DOUG_abort('[DOUG main] : Unrecognised input type.', -1)
00074     end select
00075   end function Distribution_NewInit
00076 
00078   subroutine Distribution_pmvm(D,y,x)
00079     type(Distribution),intent(inout) :: D !< Mesh and system matrix
00080     float(kind=rk),dimension(:),pointer          :: x ! Vector
00081     float(kind=rk),dimension(:),pointer          :: y ! Result
00082     integer :: i, n, p
00083 
00084     if (numprocs==1) then
00085       call SpMtx_Ax(y,D%A,x,dozero=.true.)
00086       return
00087     endif
00088     ! Initialise auxiliary data structures
00089     ! to assist with pmvm
00090     if (.not.D%cache%D_PMVM_AUXARRS_INITED) &
00091          call pmvmCommStructs_init(D%A,D%mesh,D%cache)
00092 
00093     if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED.or.&
00094          sctls%input_type==DISTRIBUTION_TYPE_STRUCTURED) then
00095       if (D%A%mtx_bbe(2,2)<D%A%nnz) then
00096         call SpMtx_pmvm_assembled_ol0(y,D%A,x,D%mesh,D%cache)
00097       else
00098         call SpMtx_pmvm_assembled(y,D%A,x,D%mesh,D%cache)
00099       endif
00100     else
00101       call SpMtx_pmvm_elemental(y,D%A,x,D%mesh,D%cache)
00102     endif
00103   end subroutine Distribution_pmvm
00104 
00105   subroutine Distribution_addoverlap(D,x)
00106     type(Distribution),intent(in) :: D
00107     float(kind=rk),dimension(:),intent(in out)   :: x ! Vector
00108 
00109     if (numprocs==1) then
00110       return
00111     endif
00112     if (sctls%input_type==DISTRIBUTION_TYPE_ELEMENTAL) then
00113       call Distribution_elem_addoverlap(D,x)
00114     else
00115       call Distribution_assm_addoverlap(D,x)
00116     endif
00117   end subroutine Distribution_addoverlap
00118 
00119 
00120   !-------------------------------------------------------
00121   ! Allocate and initialise data structures used to assist
00122   ! in parallel sparse matrix-vector multiplication during
00123   ! communications with neighbours
00124   !-------------------------------------------------------
00125   subroutine pmvmCommStructs_init(A, M, C)
00126     implicit none
00127 
00128     type(SpMtx), intent(in) :: A ! System matrix
00129     type(Mesh),  intent(in) :: M ! Mesh
00130     type(OperationCache),  intent(inout) :: C
00131 
00132     integer, dimension(:,:), pointer :: booked
00133     integer,   dimension(:), pointer :: counters
00134     integer :: p,j,h,lf,gf,ge,ptn,indx,n,f,mx
00135 
00136     if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED.or.&
00137          sctls%input_type==DISTRIBUTION_TYPE_STRUCTURED) then
00138         allocate(C%inbufs(M%nnghbrs), C%outbufs(M%nnghbrs))
00139         do p = 1,M%nnghbrs
00140           mx=max(M%ax_recvidx(p)%ninds,&
00141                  M%ol_solve(p)%ninds)
00142           allocate(C%inbufs(p)%arr(mx))
00143           mx=max(M%ax_sendidx(p)%ninds,&
00144                  M%ol_solve(p)%ninds)
00145           allocate(C%outbufs(p)%arr(mx))
00146         enddo
00147         call Vect_buildDotMask(M)
00148         C%D_PMVM_AUXARRS_INITED = .true.
00149       return
00150     endif
00151     if (numprocs==1) then
00152       C%D_PMVM_AUXARRS_INITED = .true.
00153       return
00154     endif
00155     ! <<<
00156     ! Fill in indexes of freedoms to be
00157     ! exchanged between processors
00158     allocate(C%fexchindx(maxval(M%nfreesend_map),M%nnghbrs))
00159     C%fexchindx = 0
00160 
00161     ! Map from processes's ids to indexes in 'C%fexchindx[:,C%pid2indx(:)]'
00162     allocate(C%pid2indx(M%nparts))
00163     C%pid2indx = 0
00164     do p = 1,M%nparts
00165        do j = 1,M%nnghbrs
00166           if (M%nghbrs(j) == p-1) then
00167              C%pid2indx(p) = j
00168              exit
00169           end if
00170        end do
00171     end do
00172 
00173     allocate(booked(M%nlf,M%nnghbrs))
00174     booked = 0
00175     allocate(counters(M%nnghbrs))
00176     counters = 0
00177     do lf = 1,M%nlf
00178        ! interface freedom
00179        if (M%inner_interf_fmask(lf) == D_FREEDOM_INTERF) then
00180           gf = M%lg_fmap(lf)
00181           h = M%hashlook(int(gf/M%hscale)+1)
00182           do while (M%hash(h,1) > 0)
00183              if (M%hash(h,1) == gf) then
00184                 ge = M%hash(h,2)
00185                 ptn = M%eptnmap(ge)
00186                 indx = C%pid2indx(ptn)
00187                 if (indx /= 0) then
00188                    if (booked(lf,indx) /= 1) then ! book this freedom
00189                       booked(lf,indx) = 1
00190                       counters(indx) = counters(indx) + 1
00191                       C%fexchindx(counters(indx),indx) = lf
00192                    end if
00193                 end if
00194              end if
00195              h = h + 1
00196           end do ! do while
00197        end if
00198     end do
00199 
00200     ! Substitute indexes according to freedoms apperence in SpMtx%perm_map
00201     n = sum(M%inner_interf_fmask) ! gives the number of interface freedoms
00202     do p = 1,M%nparts
00203        if (M%nfreesend_map(p) /= 0) then
00204           do indx = 1,M%nfreesend_map(p)
00205              do f = 1,n
00206                 if (A%perm_map(f) == C%fexchindx(indx,C%pid2indx(p))) then
00207                    C%fexchindx(indx,C%pid2indx(p)) = f
00208                 end if
00209              end do
00210           end do
00211        end if
00212     end do
00213     !write(stream, *) 'A%perm_map=',A%perm_map
00214 
00215     ! Bufers for incoming and outgoing messages
00216     allocate(C%inbufs(M%nnghbrs), C%outbufs(M%nnghbrs))
00217     do p = 1,M%nparts
00218        if (M%nfreesend_map(p) /= 0) then
00219           j = C%pid2indx(p)
00220           n = M%nfreesend_map(p)
00221           allocate( C%inbufs(j)%arr(n))
00222           allocate(C%outbufs(j)%arr(n))
00223        end if
00224     end do
00225 
00226     ! Auxiliary arrays has been initialised
00227     C%D_PMVM_AUXARRS_INITED = .true.
00228 
00229     deallocate(counters, booked)
00230 
00231   end subroutine pmvmCommStructs_init
00232 
00233   !------------------------------------
00234   ! Deallocate data structures used to
00235   ! assist with pmvm
00236   !------------------------------------
00237   subroutine pmvmCommStructs_destroy(C)
00238     type(OperationCache),  intent(inout) :: C
00239 
00240     integer :: i
00241 
00242     if (C%D_PMVM_AUXARRS_INITED) then
00243 
00244       if (associated(C%fexchindx)) deallocate(C%fexchindx)
00245       if (associated(C%pid2indx))  deallocate(C%pid2indx)
00246 
00247       ! Destroy incoming buffers
00248       if (associated(C%inbufs)) then
00249          do i = 1,size(C%inbufs)
00250             if (associated(C%inbufs(i)%arr)) deallocate(C%inbufs(i)%arr)
00251          end do
00252          deallocate(C%inbufs)
00253       end if
00254 
00255       ! Destroy outgoing buffers
00256       if (associated(C%outbufs)) then
00257          do i = 1,size(C%outbufs)
00258             if (associated(C%outbufs(i)%arr)) deallocate(C%outbufs(i)%arr)
00259          end do
00260          deallocate(C%outbufs)
00261       end if
00262 
00263       C%D_PMVM_AUXARRS_INITED = .false.
00264     end if
00265 
00266   end subroutine pmvmCommStructs_destroy
00267 
00268 end module Distribution_mod