source: src/datatypes/Distribution.F90 @ e9a521c

Revision e9a521c, 11.8 KB checked in by Oleg Batrashev <ogbash@…>, 2 years ago (diff)

write out neighbour info
  • Property mode set to 100644
Line 
1module Distribution_mod
2  use SpMtx_aggregation
3
4  Implicit none
5
6#include<doug_config.h>
7
8! "on-the-fly" real/complex picking
9#ifdef D_COMPLEX
10#define float complex
11#else
12#define float real
13#endif
14
15contains
16  subroutine SpMtx_DistributeAssembled(A,b,A_ghost,M,aggr)
17    use Graph_class
18    use Mesh_class
19    implicit none
20
21    type(SpMtx),intent(inout)           :: A,A_ghost
22    float(kind=rk),dimension(:),pointer :: b
23    type(Mesh)                          :: M
24    type(AggrInfo),intent(out)          :: aggr
25    !-----------------------
26    integer :: i,j,k,ierr,n,ol
27    integer, dimension(:), pointer :: xadj
28    integer, dimension(:), pointer :: adjncy
29    integer                        :: nedges
30    type(Graph) :: G
31    integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/)
32    integer,dimension(:),pointer       :: clrorder,clrstarts
33    integer, dimension(:), allocatable :: ccount !count colors
34    integer,dimension(4)           :: buf
35    float(kind=rk),dimension(:),pointer :: b_tmp
36    integer, parameter :: ONLY_METIS = 1 ! graph partitioner
37    integer, parameter :: AGGRND_METIS = 2 ! rough aggregaion based on n random seeds
38    integer,dimension(:),pointer       :: tmpnum,owner
39                                    !   going back for
40    integer :: partitioning
41    float(kind=rk) :: strong_conn1
42    integer :: aggr_radius1,min_asize1,max_asize1
43
44    if (numprocs>1) then
45      if (sctls%radius1>=0) then
46        partitioning=AGGRND_METIS
47      else
48        sctls%radius1=-sctls%radius1
49        partitioning=ONLY_METIS
50      endif
51    else
52      partitioning=ONLY_METIS
53    endif
54    ol=max(sctls%overlap,sctls%smoothers)
55    if (ismaster()) then ! Here master simply splits the matrix into pieces
56                         !   using METIS
57      if (partitioning==AGGRND_METIS) then
58        if (sctls%strong1/=0.0_rk) then
59          strong_conn1=sctls%strong1
60        else
61          strong_conn1=0.67_rk
62        endif
63        call SpMtx_find_strong(A,strong_conn1)
64        if (sctls%radius1>0) then
65          aggr_radius1=sctls%radius1+1
66        else
67          aggr_radius1=3
68        endif
69        if (sctls%minasize1>0) then
70          min_asize1=sctls%minasize1
71        else
72           min_asize1=0.5_rk*(2*aggr_radius1+1)**2
73        endif
74        if (sctls%maxasize1>0) then
75          max_asize1=sctls%maxasize1
76        else
77          max_asize1=(2*aggr_radius1+1)**2
78        endif
79        call SpMtx_roughly_aggregate(A=A,      &
80                         aggr=aggr, &
81                         neighood=aggr_radius1,&
82                      maxaggrsize=max_asize1,  &
83                            alpha=strong_conn1)
84        call SpMtx_buildAggrAdjncy(A,aggr,max_asize1,nedges,xadj,adjncy)
85        call SpMtx_unscale(A) !todo -- check
86        G=Graph_newInit(aggr%inner%nagr,nedges,xadj,adjncy,D_GRAPH_NODAL)
87if (sctls%plotting==1.or.sctls%plotting==3) then
88 allocate(owner(A%nrows))
89 do i=1,A%nnz
90  if (A%indi(i)==A%indj(i)) then
91    if (A%val(i)<1.0d0) then
92      owner(A%indi(i))=1
93    else
94      owner(A%indi(i))=2
95    endif
96  endif
97 enddo
98 allocate(tmpnum(A%nrows))
99 tmpnum=(/(i,i=1,A%nrows)/)
100 write(stream,*)'Initial Matrix structure -- 1: a(i,i)<1; 2: a(i,i)>=1 '
101 call color_print_aggrs(n=A%nrows,aggrnum=owner)
102 write(stream,*)'...with numbering as follows: '
103 call color_print_aggrs(n=A%nrows,aggrnum=tmpnum,owner=owner)
104 deallocate(tmpnum,owner)
105endif
106      elseif (partitioning==ONLY_METIS) then
107        call SpMtx_buildAdjncy(A,nedges,xadj,adjncy)
108        G=Graph_newInit(A%nrows,nedges,xadj,adjncy,D_GRAPH_NODAL)
109      endif
110      ! Deallocate temporary arrays
111      if (associated(xadj))   deallocate(xadj)
112      if (associated(adjncy)) deallocate(adjncy)
113      call Graph_Partition(G,numprocs,D_PART_VKMETIS,part_opts)
114      if (partitioning==AGGRND_METIS) then
115        do i=1,A%nrows
116          aggr%inner%num(i)=G%part(aggr%inner%num(i))
117        enddo
118        if (sctls%debug==1234) then
119          open(77,FILE='domnums.txt',FORM='FORMATTED',STATUS='new')
120          do i=1,A%nrows
121            write(77,*) aggr%inner%num(i)
122          enddo
123          close(77)
124        endif
125        if (sctls%debug==4321) then
126          open(77,FILE='domnums.txt',FORM='FORMATTED',STATUS='old')
127          do i=1,A%nrows
128            read(77,FMT=*) aggr%inner%num(i)
129          enddo
130          close(77)
131        endif
132        if (sctls%plotting==1.or.sctls%plotting==3) then
133          allocate(tmpnum(A%nrows))
134          tmpnum=(/(i,i=1,A%nrows)/)
135          call color_print_aggrs(n=A%nrows,aggrnum=tmpnum,owner=aggr%inner%num)
136          deallocate(tmpnum)
137          call flush(stream)
138        endif
139      elseif (partitioning==ONLY_METIS) then
140        if (sctls%plotting==1.or.sctls%plotting==3) then
141          write(stream,*)'Rough aggregates:'
142          call color_print_aggrs(A%nrows,G%part,overwrite=.false.)
143        endif
144      endif
145      call flush(stream)
146      buf(1)=A%nrows
147      buf(2)=A%ncols
148      buf(3)=A%nnz
149      buf(4)=numprocs
150      ! Save result in Mesh object
151      M=Mesh_newInit(nell=A%nrows,ngf=A%nrows,nsd=-2,mfrelt=-1,nnode=A%nrows)
152      M%parted  = G%parted
153      M%nparts  = G%nparts
154      !call Mesh_allocate(M,eptnmap=.true.)
155      allocate(M%eptnmap(A%nrows))
156      if (partitioning==AGGRND_METIS) then
157        M%eptnmap(1:A%nrows) = aggr%inner%num(1:A%nrows)
158      elseif (partitioning==ONLY_METIS) then
159        M%eptnmap(1:A%nrows) = G%part(1:A%nrows)
160      endif
161      call Graph_Destroy(G)
162    endif
163    ! Distribute assembled matrix; first distribute essential matrix parameters and then distribute contents
164    call MPI_BCAST(buf,4,MPI_INTEGER,D_MASTER,MPI_COMM_WORLD,ierr)
165    if (.not.ismaster()) then
166      A = SpMtx_newInit(nnz=buf(3),nblocks=sctls%number_of_blocks, &
167                        nrows=buf(1),                              &
168                        ncols=buf(2),                              &
169                        symmstruct=sctls%symmstruct,               &
170                        symmnumeric=sctls%symmnumeric              &
171                       )
172      allocate(b(A%nrows))
173      M=Mesh_newInit(nell=A%nrows,ngf=A%nrows,nsd=-2,mfrelt=-1,    &
174                     nnode=A%nrows)
175      allocate(M%eptnmap(A%nrows))
176      M%parted  = .true.
177      M%nparts  = buf(4)
178    endif
179    call MPI_BCAST(M%eptnmap,A%nrows,MPI_INTEGER,D_MASTER,&
180                   MPI_COMM_WORLD,ierr)
181    if (sctls%verbose>3.and.A%nrows<200) then
182      write(stream,*)'A orig:'
183      call SpMtx_printRaw(A)
184    endif
185    call SpMtx_distributeWithOverlap(A, b, M, ol)
186
187    !========= count color elements ============
188    if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
189      n=A%nrows
190    elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
191      n=A%ncols
192    else
193      call DOUG_abort('[SpMtx_DistributeAssembled] : matrix not arranged')
194    endif
195
196    allocate(ccount(numprocs))
197    ccount=0
198    do i=1,n
199      ccount(M%eptnmap(i))=ccount(M%eptnmap(i))+1
200    enddo
201    allocate(clrstarts(numprocs+1))
202    clrstarts(1)=1
203    do i=1,numprocs
204      clrstarts(i+1)=clrstarts(i)+ccount(i)
205    end do
206    allocate(clrorder(n))
207    ccount(1:numprocs)=clrstarts(1:numprocs)
208    do i=1,n
209      clrorder(ccount(M%eptnmap(i)))=i
210      ccount(M%eptnmap(i))=ccount(M%eptnmap(i))+1
211    enddo
212    if (sctls%verbose>3.and.A%nrows<200) then
213      do i=1,numprocs                                                     !
214        write(stream,*)'partition ',i,' is in:', &                        !
215          clrorder(clrstarts(i):clrstarts(i+1)-1)                     !
216      enddo                                                               !
217    endif
218    deallocate(ccount)
219
220    !-------------------------------------------------------------------+
221    if (sctls%verbose>3.and.A%nrows<200) then
222      write(stream,*)'A after arrange:'
223      call SpMtx_printRaw(A)
224    endif
225    call SpMtx_build_ghost(myrank+1,ol,&
226                             A,A_ghost,M,clrorder,clrstarts)
227    if (sctls%verbose>3.and.A%nrows<300) then
228      write(stream,*)'A interf(1,1):'
229      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,1),endnz=A%mtx_bbe(1,1))
230      write(stream,*)'A interf(1,2):'
231      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,2),endnz=A%mtx_bbe(1,2))
232      write(stream,*)'A interf(2,1):'
233      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,1),endnz=A%mtx_bbe(2,1))
234      write(stream,*)'A inner:'
235      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,2),endnz=A%mtx_bbe(2,2))
236      if (ol>0) then
237        write(stream,*)'A ghost:'
238        call SpMtx_printRaw(A_ghost)
239      endif
240      if (A%nnz>A%mtx_bbe(2,2)) then
241        write(stream,*)'A additional in case of ol==0:'
242        call SpMtx_printRaw(A=A,startnz=A%mtx_bbe(2,2)+1,endnz=A%ol0nnz)
243      endif
244    endif
245    ! print neighbours
246    if (sctls%verbose>=1) then
247       write(stream,"(A,I0)") "N neighbours: ", M%nnghbrs
248       do i=1,M%nnghbrs
249          !write(stream,"(A,I0,A,I0,A,I0)") "neighbour ", i, ": ", M%nghbrs(i), ", overlap: ", M%ol_solve(i)%ninds
250          write(stream,"(I0,A,I0,A)",advance="no") M%nghbrs(i)+1, " ", M%ol_inner(i)%ninds+M%ol_outer(i)%ninds, " "
251       end do
252       write(stream,*) ""
253    end if
254
255    ! Localise A:
256    if (ol<=0) then
257      M%ninonol=M%ntobsent
258      M%indepoutol=M%ninner
259    endif
260    call SpMtx_Build_lggl(A,A_ghost,M)
261    if (sctls%verbose>3) then
262      write(stream,*)'tobsent:',M%lg_fmap(1:M%ntobsent)
263      write(stream,*)'...nintol:',M%lg_fmap(M%ntobsent+1:M%ninonol)
264      write(stream,*)'...nninner:',M%lg_fmap(M%ninonol+1:M%ninner)
265      write(stream,*)'...indepoutol:',M%lg_fmap(M%ninner+1:M%indepoutol)
266      write(stream,*)'...ghost-freds:',M%lg_fmap(M%indepoutol+1:M%nlf)
267    endif
268    ! Rebuild RHS vector to correspond to local freedoms
269    allocate(b_tmp(M%nlf))
270    do i=1,M%nlf
271      b_tmp(i)=b(M%lg_fmap(i))
272    end do
273    deallocate(b)
274    b=>b_tmp
275    ! Localise matrices and communication arrays
276    do k=1,M%nnghbrs
277      M%ax_recvidx(k)%inds=M%gl_fmap(M%ax_recvidx(k)%inds)
278      M%ax_sendidx(k)%inds=M%gl_fmap(M%ax_sendidx(k)%inds)
279    enddo
280    do k=1,M%nnghbrs
281      M%ol_inner(k)%inds=M%gl_fmap(M%ol_inner(k)%inds)
282      M%ol_outer(k)%inds=M%gl_fmap(M%ol_outer(k)%inds)
283      M%ol_solve(k)%inds=M%gl_fmap(M%ol_solve(k)%inds)
284    enddo
285    do i=1,A%ol0nnz
286      A%indi(i)=M%gl_fmap(A%indi(i))
287      A%indj(i)=M%gl_fmap(A%indj(i))
288    enddo
289    A%nrows=max(0, maxval(A%indi(1:A%nnz)))
290    A%ncols=max(0, maxval(A%indj))
291    A%arrange_type=D_SpMTX_ARRNG_NO
292    if(associated(A%m_bound)) deallocate(A%m_bound) ! without this A_tmp got wrong size of M_bound in pcg()
293   
294    if (ol>0) then
295      do i=1,A_ghost%nnz
296        A_ghost%indi(i)=M%gl_fmap(A_ghost%indi(i))
297        A_ghost%indj(i)=M%gl_fmap(A_ghost%indj(i))
298      enddo
299      A_ghost%nrows=max(0, maxval(A_ghost%indi))
300      A_ghost%ncols=max(0, maxval(A_ghost%indj))
301      call SpMtx_arrange(A_ghost,D_SpMtx_ARRNG_ROWS,sort=.true.)
302    endif
303    if (sctls%verbose>3.and.A%nrows<200) then
304      write(stream,*)'Localised A interf(1,1):'
305      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,1),endnz=A%mtx_bbe(1,1))
306      write(stream,*)'Localised A interf(1,2):'
307      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(1,2),endnz=A%mtx_bbe(1,2))
308      write(stream,*)'Localised A interf(2,1):'
309      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,1),endnz=A%mtx_bbe(2,1))
310      write(stream,*)'Localised A inner:'
311      call SpMtx_printRaw(A=A,startnz=A%mtx_bbs(2,2),endnz=A%mtx_bbe(2,2))
312      if (ol>0) then
313        write(stream,*)'Localised A ghost:'
314        call SpMtx_printRaw(A_ghost)
315      endif
316      if (A%nnz>A%mtx_bbe(2,2)) then
317        write(stream,*)'localised A additional in case of ol==0:'
318        call SpMtx_printRaw(A=A,startnz=A%mtx_bbe(2,2)+1,endnz=A%ol0nnz)
319      endif
320      write(stream,*)'gl_fmap:',M%gl_fmap
321      write(stream,*)'gl_fmap(lg_fmap):',M%gl_fmap(M%lg_fmap)
322      write(stream,*)'lg_fmap:',M%lg_fmap
323      !call MPI_BARRIER(MPI_COMM_WORLD,ierr)
324      !call DOUG_abort('testing nodal graph partitioning',0)
325    endif
326  end subroutine SpMtx_DistributeAssembled
327
328end module Distribution_mod
Note: See TracBrowser for help on using the repository browser.