DOUG 0.2

Distribution_assm.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 module Distribution_assm_mod
00023   use SpMtx_class
00024   use SpMtx_aggregation
00025   use SpMtx_distribution_mod
00026   use Graph_class
00027   use SpMtx_util
00028   use Distribution_base_mod
00029   use Vect_mod
00030 
00031   implicit none
00032 
00033 #include<doug_config.h>
00034 
00035 ! "on-the-fly" real/complex picking
00036 #ifdef D_COMPLEX
00037 #define float complex
00038 #else
00039 #define float real
00040 #endif
00041 
00042   private
00043   public :: parallelDistributeAssembledInput, Distribution_assm_addoverlap
00044 
00045 contains
00046   !----------------------------------------------------------------
00049   subroutine parallelDistributeAssembledInput(Msh, A, b, A_interf)
00050     implicit none
00051 
00052     type(Mesh),     intent(in out) :: Msh !< Mesh
00053     type(SpMtx),    intent(in out) :: A !< System matrix
00054     float(kind=rk), dimension(:), pointer :: b !< local RHS
00055     type(SpMtx),intent(in out),optional :: A_interf !< matrix at interface
00056 
00057     type(AggrInfo) :: aggr
00058     integer :: n
00059 
00060     aggr = AggrInfo_New()
00061 
00062     ! ======================
00063     ! Read matrix from file
00064     ! ======================
00065     if (ismaster()) then
00066       write(stream,'(a,a)') ' ##### Assembled input file: ##### ', &
00067             mctls%assembled_mtx_file
00068       call ReadInSparseAssembled(A,trim(mctls%assembled_mtx_file))
00069       allocate(b(A%nrows))
00070       if (len_trim(mctls%assembled_rhs_file)>0) then
00071         write(stream,'(a,a)') ' ##### Assembled RHS file: ##### ', &
00072               mctls%assembled_rhs_file
00073         call Vect_ReadFromFile(b, trim(mctls%assembled_rhs_file), mctls%assembled_rhs_format)
00074       else
00075         b=1.0_rk
00076       end if
00077     endif
00078 
00079     ! =====================
00080     ! Build mesh structure/distribute
00081     ! =====================
00082     if (numprocs==1) then
00083       n=sqrt(1.0_rk*A%nrows)
00084       if (n*n /= A%nrows) then
00085         write (stream,*) 'Not a Cartesian Mesh!!!'
00086         Msh=Mesh_New()
00087         Msh%ngf=A%nrows
00088         Msh%nlf=A%nrows
00089         Msh%ninner=Msh%ngf
00090       else
00091         write (stream,*) 'Cartesian Mesh!!!'
00092         call Mesh_BuildSquare(Msh,n)
00093         Msh%ninner=Msh%ngf
00094       endif
00095     else ! numprocs>1
00096       Msh=Mesh_New()
00097       call SpMtx_DistributeAssembled(A,b,Msh,aggr)
00098       call SpMtx_localize(A,A_interf,b,Msh)
00099     endif
00100 
00101     call AggrInfo_Destroy(aggr)
00102   end subroutine parallelDistributeAssembledInput
00103 
00105   subroutine SpMtx_DistributeAssembled(A,b,M,aggr)
00106     implicit none
00107 
00108     type(SpMtx),intent(inout)           :: A
00109     float(kind=rk),dimension(:),pointer :: b
00110     type(Mesh)                          :: M
00111     type(AggrInfo),intent(out)          :: aggr
00112     !-----------------------
00113     integer :: i,j,k,ierr,ol
00114     integer, dimension(:), pointer :: xadj
00115     integer, dimension(:), pointer :: adjncy
00116     integer                        :: nedges
00117     type(Graph) :: G
00118     integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/)
00119     integer,dimension(4)           :: buf
00120     integer, parameter :: ONLY_METIS = 1 ! graph partitioner
00121     integer, parameter :: AGGRND_METIS = 2 ! rough aggregaion based on n random seeds
00122     integer,dimension(:),pointer       :: tmpnum,owner
00123                                     !   going back for 
00124     integer :: partitioning
00125     float(kind=rk) :: strong_conn1
00126     integer :: aggr_radius1,min_asize1,max_asize1
00127 
00128     if (numprocs>1) then
00129       if (sctls%radius1>=0) then
00130         partitioning=AGGRND_METIS
00131       else
00132         sctls%radius1=-sctls%radius1
00133         partitioning=ONLY_METIS
00134       endif
00135     else
00136       partitioning=ONLY_METIS
00137     endif
00138     ol=max(sctls%overlap,sctls%smoothers)
00139     if (ismaster()) then ! Here master simply splits the matrix into pieces
00140                          !   using METIS
00141       if (partitioning==AGGRND_METIS) then
00142         if (sctls%strong1/=0.0_rk) then
00143           strong_conn1=sctls%strong1
00144         else
00145           strong_conn1=0.67_rk
00146         endif
00147         call SpMtx_find_strong(A,strong_conn1)
00148         if (sctls%radius1>0) then
00149           aggr_radius1=sctls%radius1+1
00150         else
00151           aggr_radius1=3
00152         endif
00153         if (sctls%minasize1>0) then
00154           min_asize1=sctls%minasize1
00155         else
00156            min_asize1=0.5_rk*(2*aggr_radius1+1)**2
00157         endif
00158         if (sctls%maxasize1>0) then
00159           max_asize1=sctls%maxasize1
00160         else
00161           max_asize1=(2*aggr_radius1+1)**2
00162         endif
00163         call SpMtx_roughly_aggregate(A=A,      &
00164                          aggr=aggr, &
00165                          neighood=aggr_radius1,&
00166                       maxaggrsize=max_asize1,  &
00167                             alpha=strong_conn1)
00168         call SpMtx_buildAggrAdjncy(A,aggr,max_asize1,nedges,xadj,adjncy)
00169         call SpMtx_unscale(A) !todo -- check
00170         G=Graph_newInit(aggr%inner%nagr,nedges,xadj,adjncy,D_GRAPH_NODAL)
00171 if (sctls%plotting==1.or.sctls%plotting==3) then
00172  allocate(owner(A%nrows))
00173  do i=1,A%nnz
00174   if (A%indi(i)==A%indj(i)) then
00175     if (A%val(i)<1.0d0) then
00176       owner(A%indi(i))=1
00177     else
00178       owner(A%indi(i))=2
00179     endif
00180   endif
00181  enddo
00182  allocate(tmpnum(A%nrows))
00183  tmpnum=(/(i,i=1,A%nrows)/)
00184  write(stream,*)'Initial Matrix structure -- 1: a(i,i)<1; 2: a(i,i)>=1 '
00185  call color_print_aggrs(n=A%nrows,aggrnum=owner)
00186  write(stream,*)'...with numbering as follows: '
00187  call color_print_aggrs(n=A%nrows,aggrnum=tmpnum,owner=owner)
00188  deallocate(tmpnum,owner)
00189 endif
00190       elseif (partitioning==ONLY_METIS) then
00191         call SpMtx_buildAdjncy(A,nedges,xadj,adjncy)
00192         G=Graph_newInit(A%nrows,nedges,xadj,adjncy,D_GRAPH_NODAL)
00193       endif
00194       ! Deallocate temporary arrays
00195       if (associated(xadj))   deallocate(xadj)
00196       if (associated(adjncy)) deallocate(adjncy)
00197       call Graph_Partition(G,numprocs,D_PART_VKMETIS,part_opts)
00198       if (partitioning==AGGRND_METIS) then
00199         do i=1,A%nrows
00200           aggr%inner%num(i)=G%part(aggr%inner%num(i))
00201         enddo
00202         if (sctls%debug==1234) then
00203           open(77,FILE='domnums.txt',FORM='FORMATTED',STATUS='new')
00204           do i=1,A%nrows
00205             write(77,*) aggr%inner%num(i)
00206           enddo
00207           close(77)
00208         endif
00209         if (sctls%debug==4321) then
00210           open(77,FILE='domnums.txt',FORM='FORMATTED',STATUS='old')
00211           do i=1,A%nrows
00212             read(77,FMT=*) aggr%inner%num(i)
00213           enddo
00214           close(77)
00215         endif
00216         if (sctls%plotting==1.or.sctls%plotting==3) then
00217           allocate(tmpnum(A%nrows))
00218           tmpnum=(/(i,i=1,A%nrows)/)
00219           call color_print_aggrs(n=A%nrows,aggrnum=tmpnum,owner=aggr%inner%num)
00220           deallocate(tmpnum)
00221           call flush(stream)
00222         endif
00223       elseif (partitioning==ONLY_METIS) then
00224         if (sctls%plotting==1.or.sctls%plotting==3) then
00225           write(stream,*)'Rough aggregates:'
00226           call color_print_aggrs(A%nrows,G%part,overwrite=.false.)
00227         endif
00228       endif
00229       call flush(stream)
00230       buf(1)=A%nrows
00231       buf(2)=A%ncols
00232       buf(3)=A%nnz
00233       buf(4)=numprocs
00234       ! Save result in Mesh object
00235       M=Mesh_newInit(nell=A%nrows,ngf=A%nrows,nsd=-2,mfrelt=-1,nnode=A%nrows)
00236       M%parted  = G%parted
00237       M%nparts  = G%nparts
00238       !call Mesh_allocate(M,eptnmap=.true.)
00239       allocate(M%eptnmap(A%nrows))
00240       if (partitioning==AGGRND_METIS) then
00241         M%eptnmap(1:A%nrows) = aggr%inner%num(1:A%nrows)
00242       elseif (partitioning==ONLY_METIS) then
00243         M%eptnmap(1:A%nrows) = G%part(1:A%nrows)
00244       endif
00245       call Graph_Destroy(G)
00246     endif
00247     ! Distribute assembled matrix; first distribute essential matrix parameters and then distribute contents
00248     call MPI_BCAST(buf,4,MPI_INTEGER,D_MASTER,MPI_COMM_WORLD,ierr)
00249     if (.not.ismaster()) then
00250       A = SpMtx_newInit(nnz=buf(3),nblocks=sctls%number_of_blocks, &
00251                         nrows=buf(1),                              &
00252                         ncols=buf(2),                              &
00253                         symmstruct=sctls%symmstruct,               &
00254                         symmnumeric=sctls%symmnumeric              &
00255                        )
00256       allocate(b(A%nrows))
00257       M=Mesh_newInit(nell=A%nrows,ngf=A%nrows,nsd=-2,mfrelt=-1,    &
00258                      nnode=A%nrows)
00259       allocate(M%eptnmap(A%nrows))
00260       M%parted  = .true.
00261       M%nparts  = buf(4)
00262     endif
00263     call MPI_BCAST(M%eptnmap,A%nrows,MPI_INTEGER,D_MASTER,&
00264                    MPI_COMM_WORLD,ierr)
00265     if (sctls%verbose>3.and.A%nrows<200) then 
00266       write(stream,*)'A orig:'
00267       call SpMtx_printRaw(A)
00268     endif
00269     call SpMtx_distributeWithOverlap(A, b, M, ol)
00270   end subroutine SpMtx_DistributeAssembled
00271 
00272   subroutine Distribution_assm_addoverlap(D,x)
00273     type(Distribution),intent(in) :: D
00274     float(kind=rk),dimension(:),intent(in out)   :: x ! Vector
00275     integer :: i,j,k,n,n2,p,ol,mx
00276     ! MPI
00277     integer, dimension(:), pointer :: in_reqs
00278     integer                        :: ierr, out_req, status(MPI_STATUS_SIZE)
00279     integer, parameter             :: D_TAG_FREE_OUTEROL = 778
00280     !logical :: takeaverage=.true.
00281     logical :: takeaverage=.false.
00282     float(kind=rk),dimension(:),pointer,save :: nowners
00283 
00284     if (numprocs==1) then
00285       return
00286     endif
00287     ol=max(sctls%overlap,sctls%smoothers)
00288     if (ol<1) then
00289       return
00290     endif
00291     if (takeaverage) then
00292       if (.not.associated(nowners)) then
00293         allocate(nowners(size(x)))
00294         nowners=1.0_rk
00295         do i=1,D%mesh%nnghbrs
00296           n=D%mesh%ol_solve(i)%ninds
00297           do j=1,n
00298             k=D%mesh%ol_solve(i)%inds(j)
00299             nowners(k)=nowners(k)+1.0_rk
00300           enddo
00301         enddo
00302       endif
00303     endif
00304     allocate(in_reqs(D%mesh%nnghbrs))
00305     ! initialise receives
00306     do i=1,D%mesh%nnghbrs
00307       n=D%mesh%ol_solve(i)%ninds
00308       p=D%mesh%nghbrs(i)
00309 !write(stream,*) '**** starting non-blocking recv from ',p
00310       call MPI_IRECV(D%cache%inbufs(i)%arr,n,MPI_fkind, &
00311                p,D_TAG_FREE_OUTEROL,MPI_COMM_WORLD,in_reqs(i),ierr)
00312     enddo
00313     ! non-blocking send:
00314     do i=1,D%mesh%nnghbrs
00315       n=D%mesh%ol_solve(i)%ninds
00316       p=D%mesh%nghbrs(i)
00317       D%cache%outbufs(i)%arr(1:n)=x(D%mesh%ol_solve(i)%inds)
00318       call MPI_ISEND(D%cache%outbufs(i)%arr,n,MPI_fkind, &
00319                p,D_TAG_FREE_OUTEROL,MPI_COMM_WORLD,out_req,ierr)
00320 !write(stream,*) '**** sending to ',p,D%cache%outbufs(i)%arr(1:n)
00321     enddo
00322 call MPI_Barrier(MPI_COMM_WORLD,ierr)!todo: remove
00323     do while (.true.)
00324       call MPI_WAITANY(D%mesh%nnghbrs,in_reqs,i,status,ierr)
00325       if (i/=MPI_UNDEFINED) then
00326         n=D%mesh%ol_solve(i)%ninds
00327 !write(stream,*)'**** received from ',D%mesh%nghbrs(i),D%cache%inbufs(i)%arr(1:n)
00328 !write(stream,*)i,'ol_solve%inds are(glob):',D%mesh%lg_fmap(D%mesh%ol_solve(i)%inds)
00329 !write(stream,*)'BBB before x(5):',x(11)
00330         x(D%mesh%ol_solve(i)%inds)=x(D%mesh%ol_solve(i)%inds)+D%cache%inbufs(i)%arr(1:n)
00331 !write(stream,*)'BBB after x(5):',x(11)
00332       else
00333         exit
00334       endif
00335 !write(stream,*)'=== the updated vector is:',x
00336     enddo
00337     if (takeaverage) then
00338       do i=1,size(x)
00339         if (nowners(i)>1.0_rk) then
00340           x(i)=x(i)/nowners(i)
00341         endif
00342       enddo
00343     endif
00344   end subroutine Distribution_assm_addoverlap
00345 
00346 end module Distribution_assm_mod