Context Navigation

SpMtx_aggregation.F90 @ bac5f50

Revision bac5f50, 33.6 KB checked in by Oleg Batrashev <ogbash@…>, 2 years ago (diff)
Fixed the 'strong' handling.
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	call setup_aggr_cdat(nagrs,n,aggrnum,M)
372	endif
373	elseif (toosmall) then ! }{
374	! build the aggregate reference structure
375	allocate(aggrsize(nagrs)) ! the initial aggr sizes
376	aggrsize=0
377	do i=1,n ! find the #nodes for each color
378	j=aggrnum(i)
379	aggrsize(j)=aggrsize(j)+1
380	enddo
381	! We need to grow/shrink dynamically aggregates during the remake
382	! => needing a quick way for structure reference
383	maxstructlen=max(maxasizelargest,maxval(aggrsize))
384	allocate(structnodes(maxstructlen,nagrs))
385	aggrsize=0
386	do i=1,n ! fill in nodes for each color
387	j=aggrnum(i)
388	aggrsize(j)=aggrsize(j)+1
389	structnodes(aggrsize(j),j)=i
390	enddo
391	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
392	! the next stage, dealing with too small structures:
393	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
394	mnstructneigs=4*maxasize
395	call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS)
396	allocate(colsaround(mnstructneigs)) ! lists the colors
397	! around too small structure
398	allocate(connweightsums(mnstructneigs)) ! weight sums to each colour!
399	nagrs_new=nagrs
400	ntoosmall=0
401	neaten=0
402	noccupied=0
403	do i=1,nagrs ! {
404	if (aggrsize(i)<minasize) then !too small aggregate
405	!print *,'too smaall aggr #:',i,aggrsize(i)
406	ntoosmall=ntoosmall+1
407	colsaround=0
408	ncolsaround=0
409	do k=1,aggrsize(i) ! {
410	! now look, if the node has coloured neighbours:
411	!rs=A%strong_rowstart(k)
412	!re=A%strong_rowstart(k+1)-1
413	! we now look all connections, not only strong ones...:
414	kkk=structnodes(k,i)
415	rs=A%M_bound(kkk)
416	re=A%M_bound(kkk+1)-1
417	!print *,' ',kkk,' -- rs,re:',rs,re
418	do j=rs,re
419	!cn=A%strong_colnrs(j)
420	cn=A%indj(j)
421	if (cn<=n) then
422	colr=aggrnum(cn)
423	if (colr>0.and.colr/=i) then ! not a node from the same structure
424	colsearc4:do kk=1,ncolsaround ! (find where to put it)
425	if (colr==colsaround(kk)) then ! that colour again!
426	connweightsums(kk)=connweightsums(kk)+dabs(A%val(j))
427	exit colsearc4
428	endif
429	enddo colsearc4
430	if (kk>ncolsaround) then ! add the colour
431	ncolsaround=ncolsaround+1
432	if (ncolsaround>mnstructneigs) then
433	write(,)'mnstructneigs too small'
434	stop
435	endif
436	connweightsums(kk)=dabs(A%val(j))
437	colsaround(kk)=colr
438	endif
439	endif
440	endif
441	enddo
442	enddo ! }
443	!print *,' ',ncolsaround,'colors around:',colsaround(1:ncolsaround),connweightsums(1:ncolsaround)
444	maxconnweightsum=0.0_rk
445	eater=0
446	do j=1,ncolsaround
447	! first look for best coloured neighbour that could swallow it:
448	if (maxconnweightsum<connweightsums(j)) then
449	!if (aggrsize(colsaround(j))+aggrsize(i)<=maxasizelargest) then ! aggregate colsaround(j)
450	if (aggrsize(colsaround(j))+aggrsize(i)<=maxasize) then ! aggregate colsaround(j)
451	! wants to eat it up
452	maxconnweightsum=connweightsums(j)
453	eater=colsaround(j)
454	endif
455	endif
456	enddo
457	if (present(aggr_fine)) then
458	print *,'too small aggr ',i,' of size:',aggrsize(i),' maxcw_sum:',maxconnweightsum
459	endif
460	if (maxconnweightsum>=alpha.or.(present(aggr_fine).and.maxconnweightsum>=beta)) then ! let the eater get the nodes
461	do j=1,aggrsize(i)
462	! print *, j, eater, lbound(aggrsize), ubound(aggrsize)
463	aggrsize(eater)=aggrsize(eater)+1
464	structnodes(aggrsize(eater),eater)=structnodes(j,i)
465	aggrnum(structnodes(j,i))=eater
466	!print *,i,' :::',eater,' has eaten node',structnodes(j,i)
467	enddo
468	aggrsize(i)=0
469	nagrs_new=nagrs_new-1
470	neaten=neaten+1
471	print *,'eater of ',i,' is:',eater
472	else ! try to give the struct away node by node to all good neighbours...
473	if (present(aggr_fine)) then
474	print *,' ...not eaten... '
475	endif
476	nleft=aggrsize(i)
477	reduced=.true.
478	do while (nleft>0.and.reduced)
479	reduced=.false.
480	do k=aggrsize(i),1,-1
481	kkk=structnodes(k,i)
482	if (kkk>0) then
483	rs=A%M_bound(kkk)
484	re=A%M_bound(kkk+1)-1
485	!print *,' ----- ',kkk,' -- rs,re:',rs,re
486	looking:do j=rs,re
487	cn=A%indj(j)
488	if (cn<=n) then
489	colr=aggrnum(cn)
490	if (colr/=i) then ! not a node from the same structure
491	!if (dabs(A%val(j))>=beta.and. & ! strongly connected
492	!if (dabs(A%val(j))>=alpha.and. & ! strongly connected
493	! aggrsize(colr)<maxasizelargest) then ! and fits in
494	!
495	!!if (aggrsize(colr)<maxasizelargest.and.( &
496	!! dabs(A%val(j))>=alpha .or. &
497	!! (present(aggr_fine).and.dabs(A%val(j))>=beta))) then
498	!
499	if (aggrsize(colr)<maxasizelargest.and. &
500	dabs(A%val(j))>=beta ) then
501	aggrsize(colr)=aggrsize(colr)+1
502	structnodes(aggrsize(colr),colr)=structnodes(k,i)
503	aggrnum(structnodes(k,i))=colr
504	structnodes(k,i)=-1
505	nleft=nleft-1
506	reduced=.true.
507	!print *,i,' ########## sold ',structnodes(aggrsize(colr),colr), &
508	! ' to:',colr,'##########'
509	exit looking ! this node is sold
510	endif
511	endif
512	endif
513	enddo looking
514	endif
515	enddo
516	enddo ! while
517	if (nleft>0.and.nleft<aggrsize(i)) then ! compress the list
518	k=0 ! of left behind nodes
519	do j=1,aggrsize(i)
520	if (structnodes(j,i)>0) then
521	k=k+1
522	if (k>0.and.k<j) then
523	structnodes(k,i)=structnodes(j,i)
524	endif
525	endif
526	enddo
527	aggrsize(i)=nleft
528	print *,' ========== aggregate ',i,' remained as small as ',nleft,'@@@@@'
529	elseif (nleft==0) then ! the aggregate got removed
530	print *,' ========== aggregate ',i,' got removed node by node ============'
531	noccupied=noccupied+1
532	aggrsize(i)=0
533	nagrs_new=nagrs_new-1
534	endif
535	endif
536	endif
537	enddo ! }
538	naggregatednodes=sum(aggrsize(1:nagrs))
539	full_nagrs_new=nagrs_new
540	!write(stream,*)'aaaa nagrs=',nagrs_new
541	!write(stream,*)'aaaa aggrnum=',aggrnum
542	! compress the local numbers if needed
543	cnt=maxval(aggrnum(1:n))
544	if (cnt>nagrs_new) then
545	allocate(stat(cnt))
546	stat=0
547	col=0
548	do i=1,n
549	if (aggrnum(i)>0) then
550	if (stat(aggrnum(i))==0) then
551	col=col+1
552	stat(aggrnum(i))=col
553	thiscol=col
554	else
555	thiscol=stat(aggrnum(i))
556	endif
557	aggrnum(i)=thiscol
558	endif
559	enddo
560	deallocate(stat)
561	! else
562	!col=nagrs_new
563	endif
564	!write(stream,*)'bbbb nagrs=',col
565	!write(stream,*)'bbbb aggrnum=',aggrnum
566	if (n==naggregatednodes) then
567	aggrarefull=.true.
568	fullaggrnum=aggrnum(1:A%nrows)
569	full_nagrs_new=maxval(fullaggrnum)
570	else
571	aggrarefull=.false.
572	do i=1,n
573	if (aggrnum(i)>0) then
574	fullaggrnum(i)=aggrnum(i)
575	else
576	full_nagrs_new=full_nagrs_new+1
577	fullaggrnum(i)=full_nagrs_new
578	endif
579	enddo
580	endif
581	call Form_Aggr(aggr=aggr%inner, &
582	nagrs=nagrs_new, &
583	n=n, &
584	radius=neighood, &
585	nisolated=n-naggregatednodes, &
586	aggrnum=aggrnum)
587	if (numprocs>1) then
588	call setup_aggr_cdat(nagrs_new,n,aggrnum,M)
589	endif
590	deallocate(connweightsums) ! weight sums to each colour!
591	deallocate(colsaround) ! lists the colors
592	deallocate(structnodes)
593	deallocate(aggrsize) ! the initial aggr sizes
594	nagrs=nagrs_new
595	write(stream,*)'# too small aggrs:',ntoosmall,' #eaten:',neaten, &
596	' #occupied:',noccupied, &
597	' # remaining:',ntoosmall-neaten-noccupied
598	endif !}
599	call Form_Aggr(aggr=aggr%full, &
600	nagrs=full_nagrs_new, &
601	n=n, &
602	radius=neighood, &
603	nisolated=nisolated, &
604	aggrnum=fullaggrnum)
605	deallocate(fullaggrnum)
606	deallocate(aggrnum)
607
608	if (plot==1.or.plot==3) then
609	if (numprocs>1) then
610	if (ismaster()) then
611	allocate(aggrnum(M%ngf))
612	allocate(owner(M%ngf))
613	end if
614	call Integer_Vect_Gather(aggr%inner%num,aggrnum,M,owner)
615	if (ismaster()) then
616	call color_print_aggrs(M%ngf,aggrnum,overwrite=.false.,owner=owner)
617	deallocate(owner,aggrnum)
618	endif
619	else
620	if (.not.present(aggr_fine)) then
621	if (plot==3) then
622	call color_print_aggrs(A%nrows,aggr%inner%num,overwrite=.true.)
623	else
624	write(stream,*)' fine aggregates:'
625	call color_print_aggrs(A%nrows,aggr%inner%num)
626	if (.not.aggrarefull) then
627	write(stream,*)' FULL fine aggregates:'
628	call color_print_aggrs(A%nrows,aggr%full%num)
629	endif
630	endif
631	else
632	if (plot==3) then
633	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%full%num,aggr%inner%num,overwrite=.true.)
634	else
635	write(stream,*)' coarse aggregates:'
636	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,aggr%inner%num)
637	if (.not.aggrarefull) then
638	write(stream,*)' FULL coarse aggregates:'
639	call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%full%num,aggr%full%num)
640	endif
641	endif
642	endif
643	endif
644	endif
645	if (plot==3) then
646	deallocate(moviecols)
647	endif
648	end subroutine SpMtx_aggregate
649
650	subroutine SpMtx_exchange_strong(A,A_ghost,M)
651	type(SpMtx), intent(in) :: A
652	type(SpMtx), intent(inout) :: A_ghost
653	type(Mesh), intent(in) :: M
654
655	call exchange_strong()
656
657	contains
658	function calcBufferSize(ninds) result(bufferSize)
659	integer, intent(in) :: ninds
660	integer :: bufferSize
661	integer :: size, ierr
662	call MPI_Pack_size(ninds, MPI_INTEGER, MPI_COMM_WORLD, size, ierr)
663	bufferSize = 2*size
664	end function calcBufferSize
665
666	function calcBufferSize2(ninds) result(bufferSize)
667	integer, intent(in) :: ninds
668	integer :: bufferSize
669	integer :: size, ierr
670	call MPI_Pack_size(ninds, MPI_LOGICAL, MPI_COMM_WORLD, size, ierr)
671	bufferSize = size
672	end function calcBufferSize2
673
674	subroutine exchange_strong()
675	integer :: k, k2, neigh, ninds, bufsize, bufpos, ptn
676	type(indlist), allocatable :: strong_sends(:), strong_recvs(:)
677	logical, allocatable :: strongval_recvs(:)
678	type Buffer
679	character, pointer :: data(:)
680	end type Buffer
681	type(Buffer),allocatable :: outbuffers(:)
682	integer, allocatable :: outreqs(:), outstatuses(:,:), indi(:), indj(:)
683	integer :: status(MPI_STATUS_SIZE), ierr
684	character, allocatable :: inbuffer(:)
685	logical, allocatable :: strong(:)
686	integer :: i,j,nnz
687
688	allocate(strong_sends(M%nnghbrs))
689
690	! report to the neighbours which elements we need
691	!! how many elements
692	strong_sends%ninds = 0
693	do k=1,A_ghost%nnz
694	ptn = M%eptnmap(M%lg_fmap(A_ghost%indi(k)))
695	if (ptn/=myrank+1) then
696	! find neighbour number
697	do i=1,M%nnghbrs
698	if (M%nghbrs(i)+1==ptn) then
699	neigh = i
700	exit
701	end if
702	end do
703	strong_sends(neigh)%ninds = strong_sends(neigh)%ninds+1
704	end if
705	end do
706	!! collect elements
707	!!! prepare indices
708	do neigh=1,M%nnghbrs
709	allocate(strong_sends(neigh)%inds(strong_sends(neigh)%ninds))
710	end do
711	strong_sends%ninds = 0 ! reset
712	do k=1,A_ghost%nnz
713	ptn = M%eptnmap(M%lg_fmap(A_ghost%indi(k)))
714	if (ptn/=myrank+1) then
715	! find neighbour number
716	do i=1,M%nnghbrs
717	if (M%nghbrs(i)+1==ptn) then
718	neigh = i
719	exit
720	end if
721	end do
722	ninds = strong_sends(neigh)%ninds
723	strong_sends(neigh)%ninds = ninds+1
724	strong_sends(neigh)%inds(ninds+1) = k
725	end if
726	end do
727
728	! gather number of matrix elements
729	allocate(strong_recvs(M%nnghbrs))
730	do i=1,M%nnghbrs
731	neigh = M%nghbrs(i)
732	call MPI_Send(strong_sends(i)%ninds, 1, MPI_INTEGER, neigh, &
733	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, ierr)
734	end do
735	do i=1,M%nnghbrs
736	neigh = M%nghbrs(i)
737	call MPI_Recv(strong_recvs(i)%ninds, 1, MPI_INTEGER, neigh, &
738	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr)
739	end do
740
741	!!! pack and send index data
742	allocate(outbuffers(M%nnghbrs))
743	allocate(outreqs(M%nnghbrs))
744	do neigh=1,M%nnghbrs
745	bufsize = calcBufferSize(strong_sends(neigh)%ninds)
746	allocate(outbuffers(neigh)%data(bufsize))
747	bufpos = 0
748	call MPI_Pack(M%lg_fmap(A_ghost%indi(strong_sends(neigh)%inds)), strong_sends(neigh)%ninds, MPI_INTEGER,&
749	outbuffers(neigh)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
750	call MPI_Pack(M%lg_fmap(A_ghost%indj(strong_sends(neigh)%inds)), strong_sends(neigh)%ninds, MPI_INTEGER,&
751	outbuffers(neigh)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr)
752	call MPI_ISend(outbuffers(neigh)%data, bufsize, MPI_CHARACTER, M%nghbrs(neigh), TAG_EXCHANGE_STRONG,&
753	MPI_COMM_WORLD, outreqs(neigh), ierr)
754	end do
755
756	!!! recv index data
757	allocate(inbuffer(calcBufferSize(maxval(strong_recvs%ninds))))
758	nnz = sum(strong_recvs%ninds)
759	allocate(indi(nnz),indj(nnz))
760	nnz = 0
761	do i=1,M%nnghbrs
762	bufsize = calcBufferSize(strong_recvs(i)%ninds)
763	call MPI_Recv(inbuffer, bufsize, MPI_CHARACTER, M%nghbrs(i),&
764	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr)
765	ninds = strong_recvs(i)%ninds
766	bufpos = 0
767	call MPI_Unpack(inbuffer, bufsize, bufpos, indi(nnz+1), ninds,&
768	MPI_INTEGER, MPI_COMM_WORLD, ierr)
769	if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
770	call MPI_Unpack(inbuffer, bufsize, bufpos, indj(nnz+1), ninds,&
771	MPI_INTEGER, MPI_COMM_WORLD, ierr)
772	if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed")
773	nnz = nnz+ninds
774	end do
775
776	! find strong values
777	allocate(strong(nnz))
778	do k=1,nnz
779	i = M%gl_fmap(indi(k))
780	j = M%gl_fmap(indj(k))
781	k2 = SpMtx_findElem(A, i, j)
782	if(k2<=0) call DOUG_abort("Matrix element not found during 'strong' exchange")
783	strong(k) = A%strong(k2)
784	end do
785
786	! wait for all data to be sent
787	allocate(outstatuses(MPI_STATUS_SIZE, M%nnghbrs))
788	call MPI_Waitall(M%nnghbrs, outreqs, outstatuses, ierr)
789
790	! send strong values back
791	nnz = 0
792	bufpos = 0
793	do i=1,M%nnghbrs
794	ninds = strong_recvs(i)%ninds
795	call MPI_ISend(strong(nnz+1), ninds, MPI_LOGICAL, M%nghbrs(i), &
796	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, outreqs(i), ierr)
797	nnz = nnz+ninds
798	end do
799
800	! receive strong values (requested value indices are in strong_sends)
801	allocate(A_ghost%strong(A_ghost%nnz))
802	nnz = sum(strong_sends%ninds)
803	allocate(strongval_recvs(nnz))
804	nnz = 0
805	do i=1,M%nnghbrs
806	ninds = strong_sends(i)%ninds
807	call MPI_Recv(strongval_recvs(nnz+1), ninds, MPI_LOGICAL,M%nghbrs(i),&
808	TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr)
809	! overwrite local strong value with remote
810	!write (stream,*) "----", strong_sends(i)%ninds, ninds, strong_sends(i)%inds, strongval_recvs(nnz+1:nnz+ninds)
811	do k=1,ninds
812	k2 = strong_sends(i)%inds(k)
813	!write(stream,*) "--k2", k2, strongval_recvs(nnz+k), size(A_ghost%strong)
814	A_ghost%strong(k2) = strongval_recvs(nnz+k)
815	end do
816
817	nnz = nnz+ninds
818	end do
819
820	! free memory
821	do neigh=1,M%nnghbrs
822	deallocate(outbuffers(neigh)%data)
823	deallocate(strong_sends(neigh)%inds)
824	end do
825	deallocate(outbuffers)
826	deallocate(strong_sends, strong_recvs)
827	end subroutine exchange_strong
828	end subroutine SpMtx_exchange_strong
829
830	!------------------------------------------------------
831	! Finding strong connections in matrix
832	!------------------------------------------------------
833	subroutine SpMtx_find_strong(A,alpha,A_ghost,symmetrise)
834	Implicit None
835	Type(SpMtx),intent(in out) :: A
836	float(kind=rk), intent(in) :: alpha
837	Type(SpMtx),intent(in out),optional :: A_ghost
838	logical,optional :: symmetrise
839	! local:
840	integer :: i,j,k,start,ending,nnz,ndiags
841	logical :: did_scale
842	logical :: simple=.false., symm=.false.
843	float(kind=rk) :: maxndiag,aa
844	did_scale=.false.
845	if (A%scaling==D_SpMtx_SCALE_NO.or.A%scaling==D_SpMtx_SCALE_UNDEF) then
846	call SpMtx_scale(A,A_ghost)
847	did_scale=.true.
848	endif
849	if (A%mtx_bbe(2,2)>0) then
850	nnz=A%mtx_bbe(2,2)
851	else
852	nnz=A%nnz
853	endif
854
855	if (simple) then
856	if (.not.associated(A%strong)) allocate(A%strong(nnz))
857	do i=1,nnz
858	if (abs(A%val(i))>=alpha) then
859	A%strong(i)=.true.
860	else
861	A%strong(i)=.false.
862	endif
863	enddo
864	else ! not the simple case:
865	ndiags=max(A%nrows,A%ncols)
866	call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.false.)
867	if (.not.associated(A%strong)) allocate(A%strong(nnz))
868	do i=1,A%nrows
869	start=A%M_bound(i)
870	ending=A%M_bound(i+1)-1
871	maxndiag=-1.0e15
872	do j=start,ending
873	if (A%indj(j)/=i) then ! not on diagonal
874	aa=abs(A%val(j))
875	if (maxndiag<aa) then
876	maxndiag=aa
877	endif
878	endif
879	enddo
880	maxndiag=maxndiag*alpha
881	do j=start,ending
882	aa=abs(A%val(j))
883	if (A%indj(j)/=i) then ! not on diagonal
884	if (aa>maxndiag) then
885	A%strong(j)=.true.
886	else
887	A%strong(j)=.false.
888	endif
889	else
890	if (A%diag(i)>maxndiag) then
891	A%strong(j)=.true.
892	else
893	A%strong(j)=.false.
894	endif
895	!write(stream,*)'diag at i is:',i,A%diag(i),A%strong(j),maxndiag
896	endif
897	enddo
898	enddo
899	endif
900	if (present(symmetrise)) then
901	symm=symmetrise
902	endif
903	call SpMtx_symm_strong(A,A_ghost,symm)
904	end subroutine SpMtx_find_strong
905
906	subroutine SpMtx_symm_strong(A,A_ghost,symm)
907	Type(SpMtx),intent(in out) :: A
908	Type(SpMtx),intent(in out),optional :: A_ghost
909	logical,intent(in) :: symm
910	logical :: mirror2ghost=.true.
911
912	integer :: i,k,nnz
913
914	if (A%mtx_bbe(2,2)>0) then
915	nnz=A%mtx_bbe(2,2)
916	else
917	nnz=A%nnz
918	endif
919
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_symm_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