|
DOUG 0.2
|
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
1.7.3-20110217