Context Navigation

SpMtx_aggregation.F90 @ b378423

Revision b378423, 33.9 KB checked in by Oleg Batrashev <ogbash@…>, 2 years ago (diff)
Factor out Aggregates from sparse matrix structure.
Property mode set to `100644`

Line
1	! DOUG - Domain decomposition On Unstructured Grids
2	! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
3	! Department of Mathematics, University of Bath
4	!
5	! This library is free software; you can redistribute it and/or
6	! modify it under the terms of the GNU Lesser General Public
7	! License as published by the Free Software Foundation; either
8	! version 2.1 of the License, or (at your option) any later version.
9	!
10	! This library is distributed in the hope that it will be useful,
11	! but WITHOUT ANY WARRANTY; without even the implied warranty of
12	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13	! Lesser General Public License for more details.
14	!
15	! You should have received a copy of the GNU Lesser General Public
16	! License along with this library; if not, write to the Free Software
17	! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18	! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
19	! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
20	! mailto:info(at)dougdevel.org)
21
22	!> Aggregation procedures
23	Module SpMtx_aggregation
24	!!--------------------------------------------------------
25	!!Arrange elements in sparse matrix
26	!!--------------------------------------------------------
27
28	use RealKind
29	use SpMtx_class
30	use SpMtx_util
31	use Mesh_class
32	use globals
33	use SpMtx_arrangement
34
35	Implicit None
36
37	#include<doug_config.h>
38
39	! "on-the-fly" real/complex picking
40	#ifdef D_COMPLEX
41	#define float complex
42	#else
43	#define float real
44	#endif
45
46	CONTAINS
47
48	!> Finding aggregates
49	subroutine SpMtx_aggregate(A,aggr,neighood,&
50	minaggrsize,maxaggrsize,alpha,aggr_fine,M,plotting)
51	use globals
52	use CoarseAllgathers
53	use Mesh_class
54	use Vect_mod
55	Implicit None
56	Type(SpMtx),intent(in out) :: A ! our matrix
57	type(AggrInfo),intent(out) :: aggr !< aggregates
58	integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood...
59	! node stat: >= 0 free and far enough (versions 1,2 only)
60	! ==-1 free but in neighood
61	! ==-2 aggregated
62	integer,intent(in),optional :: minaggrsize,maxaggrsize
63	float(kind=rk),intent(in) :: alpha
64	Type(AggrInfo),intent(in),optional :: aggr_fine ! fine level aggregates
65	!type(CoarseData),optional :: cdat !coarse data
66	type(Mesh),optional :: M ! Mesh
67	integer,optional :: plotting
68	!-----
69	integer,dimension(:),allocatable :: aggrneigs
70	integer,dimension(:),pointer :: stat,distance
71	! Helper arrays for quich matrix row reference:
72	integer :: nneigs
73	integer :: nagrs
74	integer, dimension(:), allocatable :: nodes
75	integer :: i,j,k,kk,kkk,node,n,nn,sn,nsa,noffdels,unaggregated,dist,nisolated,ni
76	integer :: nisoneigs,nfullisoneigs,mnstructneigs
77	integer :: rs,re,cn,colr,fullj
78	integer,dimension(:),pointer :: aggrnum,fullaggrnum ! aggregate # for each node
79	integer,dimension(:),pointer :: moviecols
80	integer,dimension(:),allocatable :: aggrstarts,aggrnodes
81	integer :: minasize,maxasize,ngoodstarts,firstwocol,startnode,mindistance
82	integer :: nouters,ok,next_start_layer,maxlayer
83	integer,dimension(:),pointer :: nextgoodstart,layerlen
84	float(kind=rk) :: maxconnweightsum,fullmaxconnweightsum
85	logical :: reduced
86	float(kind=rk) :: beta=1.0E-5_rk
87	logical :: track_print=.false.
88	!logical :: track_print=.true.
89	logical :: aggrarefull
90	! for version 4:
91	logical :: toosmall
92	float(kind=rk),dimension(:),pointer :: connweightsums ! used in version 4
93	integer :: ncolsaround,agrisize,maxstructlen,eater,nleft
94	integer :: nagrs_new,full_nagrs_new,naggregatednodes,maxasizelargest
95	integer :: ntoosmall,neaten,noccupied
96	integer,dimension(:),pointer :: aggrsize,colsaround
97	integer,dimension(:,:),pointer :: structnodes
98	integer :: cnt,col,thiscol
99	integer,dimension(:),pointer :: owner
100	integer :: plot
101
102	aggr = AggrInfo_New()
103
104	toosmall=.false.
105	n=A%nrows
106	if (present(plotting)) then
107	plot=plotting
108	else
109	plot=sctls%plotting
110	endif
111	if (plot==3) then
112	track_print=.true.
113	allocate(moviecols(A%nrows))
114	endif
115	if (numprocs>1) then
116	nn=M%nlf
117	else
118	nn=n
119	endif
120	allocate(aggrnum(nn))
121	allocate(fullaggrnum(n))
122	aggrnum=0
123	sn=sqrt(1.0*n)
124	if (alpha>=0.or.sn**2/=n) then !{ else do Cartesian aggr...
125	if (.not.associated(A%strong_rowstart)) then
126	call SpMtx_build_refs_symm(A,noffdels, &
127	A%strong_rowstart,A%strong_colnrs,sortdown=.true.)
128	endif
129	if (plot==1.or.plot==3) then
130	if (present(aggr_fine)) then
131	write (stream,*) 'Coarse level aggregates:'
132	else
133	write (stream,*) 'Fine level aggregates:'
134	endif
135	endif
136	nagrs=0
137	allocate(stat(n))
138	allocate(nodes(n))
139	allocate(aggrneigs(n))
140	aggrneigs=0
141	allocate(distance(n))
142	distance=0 ! 0 -- free
143	! >0 -- aggregated with distance DISTANCE-1 FROM SEED
144	stat=0 ! 0 -- free
145	!<0 -- -weight in case of finding rounders
146	! D_AGGREGATED -- aggregated node
147	! D_PENDING -- not fitting in large enough an aggregate here
148	!>0 -- shows LAYER_NUMBER+1
149	! in general, if >0 then considered as in an aggregate
150	if (present(minaggrsize)) then
151	minasize=minaggrsize
152	else
153	minasize=neighood
154	endif
155	if (present(maxaggrsize)) then
156	maxasize=maxaggrsize
157	else
158	if (neighood==1) then
159	maxasize=9
160	else
161	maxasize=(2neighood+1)*2
162	endif
163	endif
164	if (present(aggr_fine)) then
165	!maxasizelargest=maxasize+32(2neighood)
166	maxasizelargest=3*maxasize
167	else
168	maxasizelargest=maxasize+4(2neighood)
169	!maxasizelargest=maxasize+8(2neighood)
170	!maxasizelargest=maxasize+2(2neighood)
171	endif
172	!beta=alpha
173	if (.not.present(aggr_fine)) then
174	! this seems to be problem-dependent:
175	beta=alpha/4.0_rk
176	!beta=alpha/2.0_rk
177	!beta=alpha/8.0_rk
178	!beta=alpha/5.0_rk
179	!beta=alpha/3.0_rk
180	endif
181	ngoodstarts=0
182	unaggregated=0
183	allocate(nextgoodstart(n)) ! todo: could be done much less
184	firstwocol=1
185	colo4:do while (firstwocol<=n)
186	! if needed, take out holes from the nextgoodstart
187	if (ngoodstarts>0) then
188	j=0
189	do i=1,ngoodstarts
190	if (stat(nextgoodstart(i))<D_PENDING) then
191	j=j+1
192	if (j<i) then
193	nextgoodstart(j)=nextgoodstart(i)
194	endif
195	endif
196	enddo
197	ngoodstarts=j
198	endif
199	if (ngoodstarts==0) then ! todo: try first some left-behind neighbour,
200	do while(stat(firstwocol)>=D_PENDING)
201	firstwocol=firstwocol+1
202	if (firstwocol>n) exit colo4 ! => all done
203	enddo
204	startnode=firstwocol
205	else
206	!startnode=nextgoodstart(ngoodstarts)
207	!startnode=nextgoodstart(1)
208	!
209	! let's look, if there are repeated goodstarts
210	startnode=0
211	ng4:do i=2,ngoodstarts
212	do j=1,i-1
213	if (nextgoodstart(i)==nextgoodstart(j)) then
214	startnode=nextgoodstart(i)
215	exit ng4
216	endif
217	enddo
218	enddo ng4
219	if (startnode==0) then
220	startnode=nextgoodstart(1)
221	!startnode=nextgoodstart(ngoodstarts)
222	endif
223	endif
224	ok=lets_colour(startnode,neighood,minasize,maxasize,nneigs,nodes,&
225	stat,distance,A%strong_rowstart,A%strong_colnrs)
226	mindistance=D_MAXINT
227	if (ok>0) then ! we can add the new aggregate {
228	nagrs=nagrs+1
229	nouters=0
230	if (ok==3) then ! { then we know the outermost layer was: 2*neighood+1
231	allocate(layerlen(2*neighood+2))
232	layerlen=0
233	maxlayer=0
234	if (.not.toosmall.and.nneigs<minasize) then
235	toosmall=.true.
236	endif
237	do j=nneigs,1,-1
238	k=stat(nodes(j))
239	if (k>maxlayer) maxlayer=k
240	if (k<=neighood+1) exit
241	layerlen(k)=layerlen(k)+1
242	enddo
243	next_start_layer=neighood+2
244	k=layerlen(neighood+2)
245	do j=neighood+3,maxlayer
246	if (k<layerlen(j)) then
247	k=layerlen(j)
248	next_start_layer=j
249	endif
250	enddo
251	if (track_print) then
252	if (present(aggr_fine)) then
253	do i=1,A%nrows
254	if (aggrnum(i)>0) then
255	moviecols(i)=aggrnum(i)
256	else
257	moviecols(i)=stat(i)
258	endif
259	enddo
260	if (nagrs<=1) then
261	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,moviecols,overwrite=.false.)
262	else
263	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,moviecols,overwrite=.true.)
264	endif
265	else
266	do i=1,A%nrows
267	if (aggrnum(i)>0) then
268	moviecols(i)=aggrnum(i)
269	else
270	moviecols(i)=stat(i)
271	endif
272	enddo
273	!call cursor0()
274	if (nagrs<=1) then
275	call color_print_aggrs(A%nrows,moviecols,overwrite=.false.)
276	else
277	call color_print_aggrs(A%nrows,moviecols,overwrite=.true.)
278	endif
279	endif
280	!call color_print_aggrs(A%nrows,distance)
281	endif
282	deallocate(layerlen)
283	elseif (ok==2) then ! }{
284	next_start_layer=neighood+2
285	if (.not.toosmall.and.nneigs<minasize) then
286	toosmall=.true.
287	endif
288	else ! }{ (ok==1 or 0)
289	next_start_layer=0
290	endif ! }
291	do j=1,nneigs
292	node=nodes(j)
293	if (stat(node)<=neighood+1.and.j<=maxasize) then
294	stat(node)=D_AGGREGATED ! mark as aggregated
295	aggrnum(node)=nagrs
296	elseif (next_start_layer>neighood+1) then
297	! find the smallest distance on the outer layer
298	if (stat(node)==next_start_layer) then
299	if (distance(node)<mindistance) then
300	mindistance=distance(node)
301	endif
302	! remember the outer layer:
303	nouters=nouters+1
304	nextgoodstart(n-nouters+1)=node ! using the tail of the arr.
305	endif
306	endif
307	enddo
308	! now find the nodes with minimal distance on the outer layer
309	do j=1,nouters
310	if (distance(nextgoodstart(n-j+1))==mindistance) then
311	ngoodstarts=ngoodstarts+1
312	nextgoodstart(ngoodstarts)=nextgoodstart(n-j+1)
313	endif
314	enddo
315	! need to clean up:
316	do j=1,nneigs
317	node=nodes(j)
318	if (stat(node)/=D_AGGREGATED) then ! could not aggregate
319	distance(node)=0
320	stat(node)=0
321	endif
322	enddo
323	else ! ok==0 }{
324	print *,'something wrong!'
325	endif !}
326	enddo colo4
327	deallocate(nextgoodstart)
328	deallocate(distance)
329	deallocate(aggrneigs)
330	deallocate(nodes)
331	deallocate(stat)
332	else !}{
333	write(stream,*) 'Building Cartesian aggregates...'
334	nsa=(sn+neighood-1)/neighood
335	k=0
336	do i=1,sn
337	do j=1,sn
338	k=k+1
339	aggrnum(k)=((i+neighood-1)/neighood-1)*nsa+(j+neighood-1)/neighood
340	enddo
341	!write(*,'(i3\)') aggrnum(k-sn+1:k)
342	!write(,)' '
343	enddo
344	nagrs=maxval(aggrnum)
345	nisolated=0
346	endif !}
347	if (sctls%debug==-5.and.n==65025) then ! read in the aggregate numbers:
348	write(stream,*)'##############################################'
349	write(stream,*)'##############################################'
350	write(stream,*)'##############################################'
351	write(stream,*)'Reading in aggregate numbers from aggr1.txt...'
352	write(stream,*)'##############################################'
353	write(stream,*)'##############################################'
354	write(stream,*)'##############################################'
355	open(34,file='aggrs1.txt',status='OLD',form='FORMATTED')
356	do i=1,n
357	read(34,FMT=*) aggrnum(i)
358	!print *,i,aggrnum(i)
359	enddo
360	close(34)
361	nagrs=maxval(aggrnum)
362	fullaggrnum=aggrnum
363	unaggregated=0
364	endif
365	if (.not.toosmall) then ! {
366	fullaggrnum=aggrnum(1:A%nrows)
367	full_nagrs_new=max(0, maxval(fullaggrnum))
368	call Form_Aggr(aggr%inner,nagrs,n,neighood,nisolated,aggrnum)
369	! communicate the neighbours' aggregate numbers and renumber:
370	if (numprocs>1) then
371	write(stream,*) "aggrnum", n, aggrnum
372	call setup_aggr_cdat(nagrs,n,aggrnum,M)
373	write(stream,*) "aggrnum", nn, aggrnum
374	call Form_Aggr(aggr%expanded,nagrs,nn,neighood,nisolated,aggrnum)
375	endif
376	elseif (toosmall) then ! }{
377	! build the aggregate reference structure
378	allocate(aggrsize(nagrs)) ! the initial aggr sizes
379	aggrsize=0
380	do i=1,n ! find the #nodes for each color
381	j=aggrnum(i)
382	aggrsize(j)=aggrsize(j)+1
383	enddo
384	! We need to grow/shrink dynamically aggregates during the remake
385	! => needing a quick way for structure reference
386	maxstructlen=max(maxasizelargest,maxval(aggrsize))
387	allocate(structnodes(maxstructlen,nagrs))
388	aggrsize=0
389	do i=1,n ! fill in nodes for each color
390	j=aggrnum(i)
391	aggrsize(j)=aggrsize(j)+1
392	structnodes(aggrsize(j),j)=i
393	enddo
394	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
395	! the next stage, dealing with too small structures:
396	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
397	mnstructneigs=4*maxasize
398	call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS)
399	allocate(colsaround(mnstructneigs)) ! lists the colors
400	! around too small structure
401	allocate(connweightsums(mnstructneigs)) ! weight sums to each colour!
402	nagrs_new=nagrs
403	ntoosmall=0
404	neaten=0
405	noccupied=0
406	do i=1,nagrs ! {
407	if (aggrsize(i)<minasize) then !too small aggregate
408	!print *,'too smaall aggr #:',i,aggrsize(i)
409	ntoosmall=ntoosmall+1
410	colsaround=0
411	ncolsaround=0
412	do k=1,aggrsize(i) ! {
413	! now look, if the node has coloured neighbours:
414	!rs=A%strong_rowstart(k)
415	!re=A%strong_rowstart(k+1)-1
416	! we now look all connections, not only strong ones...:
417	kkk=structnodes(k,i)
418	rs=A%M_bound(kkk)
419	re=A%M_bound(kkk+1)-1
420	!print *,' ',kkk,' -- rs,re:',rs,re
421	do j=rs,re
422	!cn=A%strong_colnrs(j)
423	cn=A%indj(j)
424	if (cn<=n) then
425	colr=aggrnum(cn)
426	if (colr>0.and.colr/=i) then ! not a node from the same structure
427	colsearc4:do kk=1,ncolsaround ! (find where to put it)
428	if (colr==colsaround(kk)) then ! that colour again!
429	connweightsums(kk)=connweightsums(kk)+dabs(A%val(j))
430	exit colsearc4
431	endif
432	enddo colsearc4
433	if (kk>ncolsaround) then ! add the colour
434	ncolsaround=ncolsaround+1
435	if (ncolsaround>mnstructneigs) then
436	write(,)'mnstructneigs too small'
437	stop
438	endif
439	connweightsums(kk)=dabs(A%val(j))
440	colsaround(kk)=colr
441	endif
442	endif
443	endif
444	enddo
445	enddo ! }
446	!print *,' ',ncolsaround,'colors around:',colsaround(1:ncolsaround),connweightsums(1:ncolsaround)
447	maxconnweightsum=0.0_rk
448	eater=0
449	do j=1,ncolsaround
450	! first look for best coloured neighbour that could swallow it:
451	if (maxconnweightsum<connweightsums(j)) then
452	!if (aggrsize(colsaround(j))+aggrsize(i)<=maxasizelargest) then ! aggregate colsaround(j)
453	if (aggrsize(colsaround(j))+aggrsize(i)<=maxasize) then ! aggregate colsaround(j)
454	! wants to eat it up
455	maxconnweightsum=connweightsums(j)
456	eater=colsaround(j)
457	endif
458	endif
459	enddo
460	if (present(aggr_fine)) then
461	print *,'too small aggr ',i,' of size:',aggrsize(i),' maxcw_sum:',maxconnweightsum
462	endif
463	if (maxconnweightsum>=alpha.or.(present(aggr_fine).and.maxconnweightsum>=beta)) then ! let the eater get the nodes
464	do j=1,aggrsize(i)
465	! print *, j, eater, lbound(aggrsize), ubound(aggrsize)
466	aggrsize(eater)=aggrsize(eater)+1
467	structnodes(aggrsize(eater),eater)=structnodes(j,i)
468	aggrnum(structnodes(j,i))=eater
469	!print *,i,' :::',eater,' has eaten node',structnodes(j,i)
470	enddo
471	aggrsize(i)=0
472	nagrs_new=nagrs_new-1
473	neaten=neaten+1
474	print *,'eater of ',i,' is:',eater
475	else ! try to give the struct away node by node to all good neighbours...
476	if (present(aggr_fine)) then
477	print *,' ...not eaten... '
478	endif
479	nleft=aggrsize(i)
480	reduced=.true.
481	do while (nleft>0.and.reduced)
482	reduced=.false.
483	do k=aggrsize(i),1,-1
484	kkk=structnodes(k,i)
485	if (kkk>0) then
486	rs=A%M_bound(kkk)
487	re=A%M_bound(kkk+1)-1
488	!print *,' ----- ',kkk,' -- rs,re:',rs,re
489	looking:do j=rs,re
490	cn=A%indj(j)
491	if (cn<=n) then
492	colr=aggrnum(cn)
493	if (colr/=i) then ! not a node from the same structure
494	!if (dabs(A%val(j))>=beta.and. & ! strongly connected
495	!if (dabs(A%val(j))>=alpha.and. & ! strongly connected
496	! aggrsize(colr)<maxasizelargest) then ! and fits in
497	!
498	!!if (aggrsize(colr)<maxasizelargest.and.( &
499	!! dabs(A%val(j))>=alpha .or. &
500	!! (present(aggr_fine).and.dabs(A%val(j))>=beta))) then
501	!
502	if (aggrsize(colr)<maxasizelargest.and. &
503	dabs(A%val(j))>=beta ) then
504	aggrsize(colr)=aggrsize(colr)+1
505	structnodes(aggrsize(colr),colr)=structnodes(k,i)
506	aggrnum(structnodes(k,i))=colr
507	structnodes(k,i)=-1
508	nleft=nleft-1
509	reduced=.true.
510	!print *,i,' ########## sold ',structnodes(aggrsize(colr),colr), &
511	! ' to:',colr,'##########'
512	exit looking ! this node is sold
513	endif
514	endif
515	endif
516	enddo looking
517	endif
518	enddo
519	enddo ! while
520	if (nleft>0.and.nleft<aggrsize(i)) then ! compress the list
521	k=0 ! of left behind nodes
522	do j=1,aggrsize(i)
523	if (structnodes(j,i)>0) then
524	k=k+1
525	if (k>0.and.k<j) then
526	structnodes(k,i)=structnodes(j,i)
527	endif
528	endif
529	enddo
530	aggrsize(i)=nleft
531	print *,' ========== aggregate ',i,' remained as small as ',nleft,'@@@@@'
532	elseif (nleft==0) then ! the aggregate got removed
533	print *,' ========== aggregate ',i,' got removed node by node ============'
534	noccupied=noccupied+1
535	aggrsize(i)=0
536	nagrs_new=nagrs_new-1
537	endif
538	endif
539	endif
540	enddo ! }
541	naggregatednodes=sum(aggrsize(1:nagrs))
542	full_nagrs_new=nagrs_new
543	!write(stream,*)'aaaa nagrs=',nagrs_new
544	!write(stream,*)'aaaa aggrnum=',aggrnum
545	! compress the local numbers if needed
546	cnt=maxval(aggrnum(1:n))
547	if (cnt>nagrs_new) then
548	allocate(stat(cnt))
549	stat=0
550	col=0
551	do i=1,n
552	if (aggrnum(i)>0) then
553	if (stat(aggrnum(i))==0) then
554	col=col+1
555	stat(aggrnum(i))=col
556	thiscol=col
557	else
558	thiscol=stat(aggrnum(i))
559	endif
560	aggrnum(i)=thiscol
561	endif
562	enddo
563	deallocate(stat)
564	! else
565	!col=nagrs_new
566	endif
567	!write(stream,*)'bbbb nagrs=',col
568	!write(stream,*)'bbbb aggrnum=',aggrnum
569	if (n==naggregatednodes) then
570	aggrarefull=.true.
571	fullaggrnum=aggrnum(1:A%nrows)
572	full_nagrs_new=maxval(fullaggrnum)
573	else
574	aggrarefull=.false.
575	do i=1,n
576	if (aggrnum(i)>0) then
577	fullaggrnum(i)=aggrnum(i)
578	else
579	full_nagrs_new=full_nagrs_new+1
580	fullaggrnum(i)=full_nagrs_new
581	endif
582	enddo
583	endif
584	call Form_Aggr(aggr=aggr%inner, &
585	nagrs=nagrs_new, &
586	n=n, &
587	radius=neighood, &
588	nisolated=n-naggregatednodes, &
589	aggrnum=aggrnum)
590	if (numprocs>1) then
591	call setup_aggr_cdat(nagrs_new,n,aggrnum,M)
592	call Form_Aggr(aggr=aggr%expanded, &
593	nagrs=nagrs_new, &
594	n=nn, &
595	radius=neighood, &
596	nisolated=n-naggregatednodes, &
597	aggrnum=aggrnum)
598	endif
599	deallocate(connweightsums) ! weight sums to each colour!
600	deallocate(colsaround) ! lists the colors
601	deallocate(structnodes)
602	deallocate(aggrsize) ! the initial aggr sizes
603	nagrs=nagrs_new
604	write(stream,*)'# too small aggrs:',ntoosmall,' #eaten:',neaten, &
605	' #occupied:',noccupied, &
606	' # remaining:',ntoosmall-neaten-noccupied
607	endif !}
608	call Form_Aggr(aggr=aggr%full, &
609	nagrs=full_nagrs_new, &
610	n=n, &
611	radius=neighood, &
612	nisolated=nisolated, &
613	aggrnum=fullaggrnum)
614	deallocate(fullaggrnum)
615	deallocate(aggrnum)
616
617	if (plot==1.or.plot==3) then
618	if (numprocs>1) then
619	if (ismaster()) then
620	allocate(aggrnum(M%ngf))
621	allocate(owner(M%ngf))
622	end if
623	call Integer_Vect_Gather(aggr%inner%num,aggrnum,M,owner)
624	if (ismaster()) then
625	call color_print_aggrs(M%ngf,aggrnum,overwrite=.false.,owner=owner)
626	deallocate(owner,aggrnum)
627	endif
628	else
629	if (.not.present(aggr_fine)) then
630	if (plot==3) then
631	call color_print_aggrs(A%nrows,aggr%inner%num,overwrite=.true.)
632	else
633	write(stream,*)' fine aggregates:'
634	call color_print_aggrs(A%nrows,aggr%inner%num)
635	if (.not.aggrarefull) then
636	write(stream,*)' FULL fine aggregates:'
637	call color_print_aggrs(A%nrows,aggr%full%num)
638	endif
639	endif
640	else
641	if (plot==3) then
642	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%full%num,aggr%inner%num,overwrite=.true.)
643	else
644	write(stream,*)' coarse aggregates:'
645	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,aggr%inner%num)
646	if (.not.aggrarefull) then
647	write(stream,*)' FULL coarse aggregates:'
648	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%full%num,aggr%full%num)
649	endif
650	endif
651	endif
652	endif
653	endif
654	if (plot==3) then
655	deallocate(moviecols)
656	endif
657	end subroutine SpMtx_aggregate
658
659	subroutine SpMtx_exchange_strong(A,A_ghost,M)
660	type(SpMtx), intent(in) :: A
661	type(SpMtx), intent(inout) :: A_ghost
662	type(Mesh), intent(in) :: M
663
664	call exchange_strong()
665
666	contains
667	function calcBufferSize(ninds) result(bufferSize)
668	integer, intent(in) :: ninds
669	integer :: bufferSize
670	integer :: size, ierr
671	call MPI_Pack_size(ninds, MPI_INTEGER, MPI_COMM_WORLD, size, ierr)
672	bufferSize = 2*size
673	end function calcBufferSize
674
675	function calcBufferSize2(ninds) result(bufferSize)
676	integer, intent(in) :: ninds
677	integer :: bufferSize
678	integer :: size, ierr
679	call MPI_Pack_size(ninds, MPI_LOGICAL, MPI_COMM_WORLD, size, ierr)
680	bufferSize = size
681	end function calcBufferSize2
682
683	subroutine exchange_strong()
684	integer :: k, k2, neigh, ninds, bufsize, bufpos, ptn
685	type(indlist), allocatable :: strong_sends(:), strong_recvs(:)
686	logical, allocatable :: strongval_recvs(:)
687	type Buffer
688	character, pointer :: data(:)
689	end type Buffer
690	type(Buffer),allocatable :: outbuffers(:)
691	integer, allocatable :: outreqs(:), outstatuses(:,:), indi(:), indj(:)
692	integer :: status(MPI_STATUS_SIZE), ierr
693	character, allocatable :: inbuffer(:)
694	logical, allocatable :: strong(:)
695	integer :: i,j,nnz
696
697	allocate(strong_sends(M%nnghbrs))
698
699	! report to the neighbours which elements we need
700	!! how many elements
701	strong_sends%ninds = 0
702	do k=1,A_ghost%nnz
703	ptn = M%eptnmap(M%lg_fmap(A_ghost%indi(k)))
704	if (ptn/=myrank+1) then
705	! find neighbour number
706	do i=1,M%nnghbrs
707	if (M%nghbrs(i)+1==ptn) then
708	neigh = i
709	exit
710	end if
711	end do
712	strong_sends(neigh)%ninds = strong_sends(neigh)%ninds+1
713	end if
714	end do
715	!! collect elements
716	!!! prepare indices
717	do neigh=1,M%nnghbrs
718	allocate(strong_sends(neigh)%inds(strong_sends(neigh)%ninds))
719	end do
720	strong_sends%ninds = 0 ! reset
721	do k=1,A_ghost%nnz
722	ptn = M%eptnmap(M%lg_fmap(A_ghost%indi(k)))
723	if (ptn/=myrank+1) then
724	! find neighbour number
725	do i=1,M%nnghbrs
726	if (M%nghbrs(i)+1==ptn) then
727	neigh = i
728	exit
729	end if
730	end do
731	ninds = strong_sends(neigh)%ninds
732	strong_sends(neigh)%ninds = ninds+1
733	strong_sends(neigh)%inds(ninds+1) = k
734	end if
735	end do
736
737	! gather number of matrix elements
738	allocate(strong_recvs(M%nnghbrs))
739	do i=1,M%nnghbrs
740	neigh = M%nghbrs(i)
741	call MPI_Send(strong_sends(i)%ninds, 1, MPI_INTEGER, neigh, &
742	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, ierr)
743	end do
744	do i=1,M%nnghbrs
745	neigh = M%nghbrs(i)
746	call MPI_Recv(strong_recvs(i)%ninds, 1, MPI_INTEGER, neigh, &
747	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr)
748	end do
749
750	!!! pack and send index data
751	allocate(outbuffers(M%nnghbrs))
752	allocate(outreqs(M%nnghbrs))
753	do neigh=1,M%nnghbrs
754	bufsize = calcBufferSize(strong_sends(neigh)%ninds)
755	allocate(outbuffers(neigh)%data(bufsize))
756	bufpos = 0
757	call MPI_Pack(M%lg_fmap(A_ghost%indi(strong_sends(neigh)%inds)), strong_sends(neigh)%ninds, MPI_INTEGER,&
758	outbuffers(neigh)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
759	call MPI_Pack(M%lg_fmap(A_ghost%indj(strong_sends(neigh)%inds)), strong_sends(neigh)%ninds, MPI_INTEGER,&
760	outbuffers(neigh)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
761	call MPI_ISend(outbuffers(neigh)%data, bufsize, MPI_CHARACTER, M%nghbrs(neigh), TAG_EXCHANGE_STRONG,&
762	MPI_COMM_WORLD, outreqs(neigh), ierr)
763	end do
764
765	!!! recv index data
766	allocate(inbuffer(calcBufferSize(maxval(strong_recvs%ninds))))
767	nnz = sum(strong_recvs%ninds)
768	allocate(indi(nnz),indj(nnz))
769	nnz = 0
770	do i=1,M%nnghbrs
771	bufsize = calcBufferSize(strong_recvs(i)%ninds)
772	call MPI_Recv(inbuffer, bufsize, MPI_CHARACTER, M%nghbrs(i),&
773	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr)
774	ninds = strong_recvs(i)%ninds
775	bufpos = 0
776	call MPI_Unpack(inbuffer, bufsize, bufpos, indi(nnz+1), ninds,&
777	MPI_INTEGER, MPI_COMM_WORLD, ierr)
778	if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
779	call MPI_Unpack(inbuffer, bufsize, bufpos, indj(nnz+1), ninds,&
780	MPI_INTEGER, MPI_COMM_WORLD, ierr)
781	if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
782	nnz = nnz+ninds
783	end do
784
785	! find strong values
786	allocate(strong(nnz))
787	do k=1,nnz
788	i = M%gl_fmap(indi(k))
789	j = M%gl_fmap(indj(k))
790	k2 = SpMtx_findElem(A, i, j)
791	if(k2<=0) call DOUG_abort("Matrix element not found during 'strong' exchange")
792	strong(k) = A%strong(k2)
793	end do
794
795	! wait for all data to be sent
796	allocate(outstatuses(MPI_STATUS_SIZE, M%nnghbrs))
797	call MPI_Waitall(M%nnghbrs, outreqs, outstatuses, ierr)
798
799	! send strong values back
800	nnz = 0
801	bufpos = 0
802	do i=1,M%nnghbrs
803	ninds = strong_recvs(i)%ninds
804	call MPI_ISend(strong(nnz+1), ninds, MPI_LOGICAL, M%nghbrs(i), &
805	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, outreqs(i), ierr)
806	nnz = nnz+ninds
807	end do
808
809	! receive strong values (requested value indices are in strong_sends)
810	allocate(A_ghost%strong(A_ghost%nnz))
811	nnz = sum(strong_sends%ninds)
812	allocate(strongval_recvs(nnz))
813	nnz = 0
814	do i=1,M%nnghbrs
815	ninds = strong_sends(i)%ninds
816	call MPI_Recv(strongval_recvs(nnz+1), ninds, MPI_LOGICAL,M%nghbrs(i),&
817	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr)
818	! overwrite local strong value with remote
819	!write (stream,*) "----", strong_sends(i)%ninds, ninds, strong_sends(i)%inds, strongval_recvs(nnz+1:nnz+ninds)
820	do k=1,ninds
821	k2 = strong_sends(i)%inds(k)
822	!write(stream,*) "--k2", k2, strongval_recvs(nnz+k), size(A_ghost%strong)
823	A_ghost%strong(k2) = strongval_recvs(nnz+k)
824	end do
825
826	nnz = nnz+ninds
827	end do
828
829	! free memory
830	do neigh=1,M%nnghbrs
831	deallocate(outbuffers(neigh)%data)
832	deallocate(strong_sends(neigh)%inds)
833	end do
834	deallocate(outbuffers)
835	deallocate(strong_sends, strong_recvs)
836	end subroutine exchange_strong
837	end subroutine SpMtx_exchange_strong
838
839	!------------------------------------------------------
840	! Finding strong connections in matrix
841	!------------------------------------------------------
842	subroutine SpMtx_find_strong(A,alpha,A_ghost,symmetrise,M)
843	Implicit None
844	Type(SpMtx),intent(in out) :: A
845	float(kind=rk), intent(in) :: alpha
846	Type(SpMtx),intent(in out),optional :: A_ghost
847	logical,intent(in),optional :: symmetrise
848	type(Mesh),intent(in),optional :: M
849	! local:
850	integer :: i,j,k,start,ending,nnz,ndiags
851	logical :: did_scale
852	logical :: simple=.false.,symm=.false.,mirror2ghost=.true.
853	float(kind=rk) :: maxndiag,aa
854	did_scale=.false.
855	if (A%scaling==D_SpMtx_SCALE_NO.or.A%scaling==D_SpMtx_SCALE_UNDEF) then
856	call SpMtx_scale(A,A_ghost)
857	did_scale=.true.
858	endif
859	if (A%mtx_bbe(2,2)>0) then
860	nnz=A%mtx_bbe(2,2)
861	else
862	nnz=A%nnz
863	endif
864	if (.not.associated(A%strong)) then
865	allocate(A%strong(nnz))
866	endif
867	if (simple) then
868	do i=1,nnz
869	if (abs(A%val(i))>=alpha) then
870	A%strong(i)=.true.
871	else
872	A%strong(i)=.false.
873	endif
874	enddo
875	else ! not the simple case:
876	ndiags=max(A%nrows,A%ncols)
877	call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.false.)
878	do i=1,A%nrows
879	start=A%M_bound(i)
880	ending=A%M_bound(i+1)-1
881	maxndiag=-1.0e15
882	do j=start,ending
883	if (A%indj(j)/=i) then ! not on diagonal
884	aa=abs(A%val(j))
885	if (maxndiag<aa) then
886	maxndiag=aa
887	endif
888	endif
889	enddo
890	maxndiag=maxndiag*alpha
891	do j=start,ending
892	aa=abs(A%val(j))
893	if (A%indj(j)/=i) then ! not on diagonal
894	if (aa>maxndiag) then
895	A%strong(j)=.true.
896	else
897	A%strong(j)=.false.
898	endif
899	else
900	if (A%diag(i)>maxndiag) then
901	A%strong(j)=.true.
902	else
903	A%strong(j)=.false.
904	endif
905	!write(stream,*)'diag at i is:',i,A%diag(i),A%strong(j),maxndiag
906	endif
907	enddo
908	enddo
909	endif
910	!if (did_scale) then
911	! call SpMtx_unscale(A)
912	!endif
913	if (present(A_ghost).and.associated(A_ghost%indi)) then
914	! this should only be called once, during local strong calculations
915	call SpMtx_exchange_strong(A,A_ghost,M)
916	end if
917	if (present(symmetrise)) then
918	symm=symmetrise
919	endif
920	if (symm) then
921	do i=1,nnz
922	if (A%arrange_type == D_SpMtx_ARRNG_ROWS) then
923	k=SpMtx_findElem(A,A%indi(i),A%indj(i))
924	else
925	k=SpMtx_findElem(A,A%indj(i),A%indi(i))
926	endif
927	if (k>0) then
928	if (A%strong(i).and..not.A%strong(k)) then
929	A%strong(k)=.true.
930	write(stream,*)'symmetrising to strong:',A%indi(i),A%indj(i)
931	elseif (A%strong(k).and..not.A%strong(i)) then
932	A%strong(i)=.true.
933	write(stream,*)'symmetrising to strong:',A%indi(i),A%indj(i)
934	endif
935	else
936	write(stream,*) 'Warning: matrix does not have symmetric structure!'
937	endif
938	enddo
939	endif
940	if (mirror2ghost) then
941	if (present(A_ghost).and.associated(A_ghost%indi)) then
942	do i=1,A_ghost%nnz
943	if(A_ghost%indi(i)/=A_ghost%indj(i)) then
944	k=SpMtx_findElem(A,A_ghost%indi(i),A_ghost%indj(i))
945	if (k>0) then
946	A_ghost%strong(i)=A%strong(k)
947	endif
948	endif
949	enddo
950	endif
951	endif
952	end subroutine SpMtx_find_strong
953
954	!------------------------------------------------------
955	End Module SpMtx_aggregation
956	!------------------------------------------------------
957

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format