|
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 00023 Module SpMtx_aggregation 00024 !!-------------------------------------------------------- 00025 !!Arrange elements in sparse matrix 00026 !!-------------------------------------------------------- 00027 00028 use RealKind 00029 use SpMtx_class 00030 use SpMtx_util 00031 use Mesh_class 00032 use globals 00033 use SpMtx_arrangement 00034 00035 Implicit None 00036 00037 #include<doug_config.h> 00038 00039 ! "on-the-fly" real/complex picking 00040 #ifdef D_COMPLEX 00041 #define float complex 00042 #else 00043 #define float real 00044 #endif 00045 00046 CONTAINS 00047 00049 subroutine SpMtx_aggregate(A,aggr,neighood,& 00050 minaggrsize,maxaggrsize,alpha,aggr_fine,M,plotting) 00051 use globals 00052 use CoarseAllgathers 00053 use Mesh_class 00054 use Vect_mod 00055 Implicit None 00056 Type(SpMtx),intent(in out) :: A ! our matrix 00057 type(AggrInfo),intent(out) :: aggr !< aggregates 00058 integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00059 ! node stat: >= 0 free and far enough (versions 1,2 only) 00060 ! ==-1 free but in neighood 00061 ! ==-2 aggregated 00062 integer,intent(in),optional :: minaggrsize,maxaggrsize 00063 float(kind=rk),intent(in) :: alpha 00064 Type(AggrInfo),intent(in),optional :: aggr_fine ! fine level aggregates 00065 !type(CoarseData),optional :: cdat !coarse data 00066 type(Mesh),optional :: M ! Mesh 00067 integer,optional :: plotting 00068 !----- 00069 integer,dimension(:),allocatable :: aggrneigs 00070 integer,dimension(:),pointer :: stat,distance 00071 ! Helper arrays for quich matrix row reference: 00072 integer :: nneigs 00073 integer :: nagrs 00074 integer, dimension(:), allocatable :: nodes 00075 integer :: i,j,k,kk,kkk,node,n,nn,sn,nsa,noffdels,unaggregated,dist,nisolated,ni 00076 integer :: nisoneigs,nfullisoneigs,mnstructneigs 00077 integer :: rs,re,cn,colr,fullj 00078 integer,dimension(:),pointer :: aggrnum,fullaggrnum ! aggregate # for each node 00079 integer,dimension(:),pointer :: moviecols 00080 integer,dimension(:),allocatable :: aggrstarts,aggrnodes 00081 integer :: minasize,maxasize,ngoodstarts,firstwocol,startnode,mindistance 00082 integer :: nouters,ok,next_start_layer,maxlayer 00083 integer,dimension(:),pointer :: nextgoodstart,layerlen 00084 float(kind=rk) :: maxconnweightsum,fullmaxconnweightsum 00085 logical :: reduced 00086 float(kind=rk) :: beta=1.0E-5_rk 00087 logical :: track_print=.false. 00088 !logical :: track_print=.true. 00089 logical :: aggrarefull 00090 ! for version 4: 00091 logical :: toosmall 00092 float(kind=rk),dimension(:),pointer :: connweightsums ! used in version 4 00093 integer :: ncolsaround,agrisize,maxstructlen,eater,nleft 00094 integer :: nagrs_new,full_nagrs_new,naggregatednodes,maxasizelargest 00095 integer :: ntoosmall,neaten,noccupied 00096 integer,dimension(:),pointer :: aggrsize,colsaround 00097 integer,dimension(:,:),pointer :: structnodes 00098 integer :: cnt,col,thiscol 00099 integer,dimension(:),pointer :: owner 00100 integer :: plot 00101 00102 aggr = AggrInfo_New() 00103 00104 toosmall=.false. 00105 n=A%nrows 00106 if (present(plotting)) then 00107 plot=plotting 00108 else 00109 plot=sctls%plotting 00110 endif 00111 if (plot==3) then 00112 track_print=.true. 00113 allocate(moviecols(A%nrows)) 00114 endif 00115 if (numprocs>1) then 00116 nn=M%nlf 00117 else 00118 nn=n 00119 endif 00120 allocate(aggrnum(nn)) 00121 allocate(fullaggrnum(n)) 00122 aggrnum=0 00123 sn=sqrt(1.0*n) 00124 if (alpha>=0.or.sn**2/=n) then !{ else do Cartesian aggr... 00125 if (.not.associated(A%strong_rowstart)) then 00126 call SpMtx_build_refs_symm(A,noffdels, & 00127 A%strong_rowstart,A%strong_colnrs,sortdown=.true.) 00128 endif 00129 if (plot==1.or.plot==3) then 00130 if (present(aggr_fine)) then 00131 write (stream,*) 'Coarse level aggregates:' 00132 else 00133 write (stream,*) 'Fine level aggregates:' 00134 endif 00135 endif 00136 nagrs=0 00137 allocate(stat(n)) 00138 allocate(nodes(n)) 00139 allocate(aggrneigs(n)) 00140 aggrneigs=0 00141 allocate(distance(n)) 00142 distance=0 ! 0 -- free 00143 ! >0 -- aggregated with distance DISTANCE-1 FROM SEED 00144 stat=0 ! 0 -- free 00145 !<0 -- -weight in case of finding rounders 00146 ! D_AGGREGATED -- aggregated node 00147 ! D_PENDING -- not fitting in large enough an aggregate here 00150 if (present(minaggrsize)) then 00151 minasize=minaggrsize 00152 else 00153 minasize=neighood 00154 endif 00155 if (present(maxaggrsize)) then 00156 maxasize=maxaggrsize 00157 else 00158 if (neighood==1) then 00159 maxasize=9 00160 else 00161 maxasize=(2*neighood+1)**2 00162 endif 00163 endif 00164 if (present(aggr_fine)) then 00165 !maxasizelargest=maxasize+32*(2*neighood) 00166 maxasizelargest=3*maxasize 00167 else 00168 maxasizelargest=maxasize+4*(2*neighood) 00169 !maxasizelargest=maxasize+8*(2*neighood) 00170 !maxasizelargest=maxasize+2*(2*neighood) 00171 endif 00172 !beta=alpha 00173 if (.not.present(aggr_fine)) then 00174 ! this seems to be problem-dependent: 00175 beta=alpha/4.0_rk 00176 !beta=alpha/2.0_rk 00177 !beta=alpha/8.0_rk 00178 !beta=alpha/5.0_rk 00179 !beta=alpha/3.0_rk 00180 endif 00181 ngoodstarts=0 00182 unaggregated=0 00183 allocate(nextgoodstart(n)) ! todo: could be done much less 00184 firstwocol=1 00185 colo4:do while (firstwocol<=n) 00186 ! if needed, take out holes from the nextgoodstart 00187 if (ngoodstarts>0) then 00188 j=0 00189 do i=1,ngoodstarts 00190 if (stat(nextgoodstart(i))<D_PENDING) then 00191 j=j+1 00192 if (j<i) then 00193 nextgoodstart(j)=nextgoodstart(i) 00194 endif 00195 endif 00196 enddo 00197 ngoodstarts=j 00198 endif 00199 if (ngoodstarts==0) then ! todo: try first some left-behind neighbour, 00200 do while(stat(firstwocol)>=D_PENDING) 00201 firstwocol=firstwocol+1 00202 if (firstwocol>n) exit colo4 ! => all done 00203 enddo 00204 startnode=firstwocol 00205 else 00206 !startnode=nextgoodstart(ngoodstarts) 00207 !startnode=nextgoodstart(1) 00208 ! 00209 ! let's look, if there are repeated goodstarts 00210 startnode=0 00211 ng4:do i=2,ngoodstarts 00212 do j=1,i-1 00213 if (nextgoodstart(i)==nextgoodstart(j)) then 00214 startnode=nextgoodstart(i) 00215 exit ng4 00216 endif 00217 enddo 00218 enddo ng4 00219 if (startnode==0) then 00220 startnode=nextgoodstart(1) 00221 !startnode=nextgoodstart(ngoodstarts) 00222 endif 00223 endif 00224 ok=lets_colour(startnode,neighood,minasize,maxasize,nneigs,nodes,& 00225 stat,distance,A%strong_rowstart,A%strong_colnrs) 00226 mindistance=D_MAXINT 00227 if (ok>0) then ! we can add the new aggregate { 00228 nagrs=nagrs+1 00229 nouters=0 00230 if (ok==3) then ! { then we know the outermost layer was: 2*neighood+1 00231 allocate(layerlen(2*neighood+2)) 00232 layerlen=0 00233 maxlayer=0 00234 if (.not.toosmall.and.nneigs<minasize) then 00235 toosmall=.true. 00236 endif 00237 do j=nneigs,1,-1 00238 k=stat(nodes(j)) 00239 if (k>maxlayer) maxlayer=k 00240 if (k<=neighood+1) exit 00241 layerlen(k)=layerlen(k)+1 00242 enddo 00243 next_start_layer=neighood+2 00244 k=layerlen(neighood+2) 00245 do j=neighood+3,maxlayer 00246 if (k<layerlen(j)) then 00247 k=layerlen(j) 00248 next_start_layer=j 00249 endif 00250 enddo 00251 if (track_print) then 00252 if (present(aggr_fine)) then 00253 do i=1,A%nrows 00254 if (aggrnum(i)>0) then 00255 moviecols(i)=aggrnum(i) 00256 else 00257 moviecols(i)=stat(i) 00258 endif 00259 enddo 00260 if (nagrs<=1) then 00261 call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,moviecols,overwrite=.false.) 00262 else 00263 call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,moviecols,overwrite=.true.) 00264 endif 00265 else 00266 do i=1,A%nrows 00267 if (aggrnum(i)>0) then 00268 moviecols(i)=aggrnum(i) 00269 else 00270 moviecols(i)=stat(i) 00271 endif 00272 enddo 00273 !call cursor0() 00274 if (nagrs<=1) then 00275 call color_print_aggrs(A%nrows,moviecols,overwrite=.false.) 00276 else 00277 call color_print_aggrs(A%nrows,moviecols,overwrite=.true.) 00278 endif 00279 endif 00280 !call color_print_aggrs(A%nrows,distance) 00281 endif 00282 deallocate(layerlen) 00283 elseif (ok==2) then ! }{ 00284 next_start_layer=neighood+2 00285 if (.not.toosmall.and.nneigs<minasize) then 00286 toosmall=.true. 00287 endif 00288 else ! }{ (ok==1 or 0) 00289 next_start_layer=0 00290 endif ! } 00291 do j=1,nneigs 00292 node=nodes(j) 00293 if (stat(node)<=neighood+1.and.j<=maxasize) then 00294 stat(node)=D_AGGREGATED ! mark as aggregated 00295 aggrnum(node)=nagrs 00296 elseif (next_start_layer>neighood+1) then 00297 ! find the smallest distance on the outer layer 00298 if (stat(node)==next_start_layer) then 00299 if (distance(node)<mindistance) then 00300 mindistance=distance(node) 00301 endif 00302 ! remember the outer layer: 00303 nouters=nouters+1 00304 nextgoodstart(n-nouters+1)=node ! using the tail of the arr. 00305 endif 00306 endif 00307 enddo 00308 ! now find the nodes with minimal distance on the outer layer 00309 do j=1,nouters 00310 if (distance(nextgoodstart(n-j+1))==mindistance) then 00311 ngoodstarts=ngoodstarts+1 00312 nextgoodstart(ngoodstarts)=nextgoodstart(n-j+1) 00313 endif 00314 enddo 00315 ! need to clean up: 00316 do j=1,nneigs 00317 node=nodes(j) 00318 if (stat(node)/=D_AGGREGATED) then ! could not aggregate 00319 distance(node)=0 00320 stat(node)=0 00321 endif 00322 enddo 00323 else ! ok==0 }{ 00324 print *,'something wrong!' 00325 endif !} 00326 enddo colo4 00327 deallocate(nextgoodstart) 00328 deallocate(distance) 00329 deallocate(aggrneigs) 00330 deallocate(nodes) 00331 deallocate(stat) 00332 else !}{ 00333 write(stream,*) 'Building Cartesian aggregates...' 00334 nsa=(sn+neighood-1)/neighood 00335 k=0 00336 do i=1,sn 00337 do j=1,sn 00338 k=k+1 00339 aggrnum(k)=((i+neighood-1)/neighood-1)*nsa+(j+neighood-1)/neighood 00340 enddo 00341 !write(*,'(i3\)') aggrnum(k-sn+1:k) 00342 !write(*,*)' ' 00343 enddo 00344 nagrs=maxval(aggrnum) 00345 nisolated=0 00346 endif !} 00347 if (sctls%debug==-5.and.n==65025) then ! read in the aggregate numbers: 00348 write(stream,*)'##############################################' 00349 write(stream,*)'##############################################' 00350 write(stream,*)'##############################################' 00351 write(stream,*)'Reading in aggregate numbers from aggr1.txt...' 00352 write(stream,*)'##############################################' 00353 write(stream,*)'##############################################' 00354 write(stream,*)'##############################################' 00355 open(34,file='aggrs1.txt',status='OLD',form='FORMATTED') 00356 do i=1,n 00357 read(34,FMT=*) aggrnum(i) 00358 !print *,i,aggrnum(i) 00359 enddo 00360 close(34) 00361 nagrs=maxval(aggrnum) 00362 fullaggrnum=aggrnum 00363 unaggregated=0 00364 endif 00365 if (.not.toosmall) then ! { 00366 fullaggrnum=aggrnum(1:A%nrows) 00367 full_nagrs_new=max(0, maxval(fullaggrnum)) 00368 call Form_Aggr(aggr%inner,nagrs,n,neighood,nisolated,aggrnum) 00369 elseif (toosmall) then ! }{ 00370 ! build the aggregate reference structure 00371 allocate(aggrsize(nagrs)) ! the initial aggr sizes 00372 aggrsize=0 00373 do i=1,n ! find the #nodes for each color 00374 j=aggrnum(i) 00375 aggrsize(j)=aggrsize(j)+1 00376 enddo 00377 ! We need to grow/shrink dynamically aggregates during the remake 00378 ! => needing a quick way for structure reference 00379 maxstructlen=max(maxasizelargest,maxval(aggrsize)) 00380 allocate(structnodes(maxstructlen,nagrs)) 00381 aggrsize=0 00382 do i=1,n ! fill in nodes for each color 00383 j=aggrnum(i) 00384 aggrsize(j)=aggrsize(j)+1 00385 structnodes(aggrsize(j),j)=i 00386 enddo 00387 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00388 ! the next stage, dealing with too small structures: 00389 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00390 mnstructneigs=4*maxasize 00391 call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS) 00392 allocate(colsaround(mnstructneigs)) ! lists the colors 00393 ! around too small structure 00394 allocate(connweightsums(mnstructneigs)) ! weight sums to each colour! 00395 nagrs_new=nagrs 00396 ntoosmall=0 00397 neaten=0 00398 noccupied=0 00399 do i=1,nagrs ! { 00400 if (aggrsize(i)<minasize) then !too small aggregate 00401 !print *,'too smaall aggr #:',i,aggrsize(i) 00402 ntoosmall=ntoosmall+1 00403 colsaround=0 00404 ncolsaround=0 00405 do k=1,aggrsize(i) ! { 00406 ! now look, if the node has coloured neighbours: 00407 !rs=A%strong_rowstart(k) 00408 !re=A%strong_rowstart(k+1)-1 00409 ! we now look all connections, not only strong ones...: 00410 kkk=structnodes(k,i) 00411 rs=A%M_bound(kkk) 00412 re=A%M_bound(kkk+1)-1 00413 !print *,' ',kkk,' -- rs,re:',rs,re 00414 do j=rs,re 00415 !cn=A%strong_colnrs(j) 00416 cn=A%indj(j) 00417 if (cn<=n) then 00418 colr=aggrnum(cn) 00419 if (colr>0.and.colr/=i) then ! not a node from the same structure 00420 colsearc4:do kk=1,ncolsaround ! (find where to put it) 00421 if (colr==colsaround(kk)) then ! that colour again! 00422 connweightsums(kk)=connweightsums(kk)+dabs(A%val(j)) 00423 exit colsearc4 00424 endif 00425 enddo colsearc4 00426 if (kk>ncolsaround) then ! add the colour 00427 ncolsaround=ncolsaround+1 00428 if (ncolsaround>mnstructneigs) then 00429 write(*,*)'mnstructneigs too small' 00430 stop 00431 endif 00432 connweightsums(kk)=dabs(A%val(j)) 00433 colsaround(kk)=colr 00434 endif 00435 endif 00436 endif 00437 enddo 00438 enddo ! } 00439 !print *,' ',ncolsaround,'colors around:',colsaround(1:ncolsaround),connweightsums(1:ncolsaround) 00440 maxconnweightsum=0.0_rk 00441 eater=0 00442 do j=1,ncolsaround 00443 ! first look for best coloured neighbour that could swallow it: 00444 if (maxconnweightsum<connweightsums(j)) then 00445 !if (aggrsize(colsaround(j))+aggrsize(i)<=maxasizelargest) then ! aggregate colsaround(j) 00446 if (aggrsize(colsaround(j))+aggrsize(i)<=maxasize) then ! aggregate colsaround(j) 00447 ! wants to eat it up 00448 maxconnweightsum=connweightsums(j) 00449 eater=colsaround(j) 00450 endif 00451 endif 00452 enddo 00453 if (present(aggr_fine)) then 00454 print *,'too small aggr ',i,' of size:',aggrsize(i),' maxcw_sum:',maxconnweightsum 00455 endif 00456 if (maxconnweightsum>=alpha.or.(present(aggr_fine).and.maxconnweightsum>=beta)) then ! let the eater get the nodes 00457 do j=1,aggrsize(i) 00458 ! print *, j, eater, lbound(aggrsize), ubound(aggrsize) 00459 aggrsize(eater)=aggrsize(eater)+1 00460 structnodes(aggrsize(eater),eater)=structnodes(j,i) 00461 aggrnum(structnodes(j,i))=eater 00462 !print *,i,' :::',eater,' has eaten node',structnodes(j,i) 00463 enddo 00464 aggrsize(i)=0 00465 nagrs_new=nagrs_new-1 00466 neaten=neaten+1 00467 print *,'eater of ',i,' is:',eater 00468 else ! try to give the struct away node by node to all good neighbours... 00469 if (present(aggr_fine)) then 00470 print *,' ...not eaten... ' 00471 endif 00472 nleft=aggrsize(i) 00473 reduced=.true. 00474 do while (nleft>0.and.reduced) 00475 reduced=.false. 00476 do k=aggrsize(i),1,-1 00477 kkk=structnodes(k,i) 00478 if (kkk>0) then 00479 rs=A%M_bound(kkk) 00480 re=A%M_bound(kkk+1)-1 00481 !print *,' ----- ',kkk,' -- rs,re:',rs,re 00482 looking:do j=rs,re 00483 cn=A%indj(j) 00484 if (cn<=n) then 00485 colr=aggrnum(cn) 00486 if (colr/=i) then ! not a node from the same structure 00487 !if (dabs(A%val(j))>=beta.and. & ! strongly connected 00488 !if (dabs(A%val(j))>=alpha.and. & ! strongly connected 00489 ! aggrsize(colr)<maxasizelargest) then ! and fits in 00490 ! 00491 !!if (aggrsize(colr)<maxasizelargest.and.( & 00492 !! dabs(A%val(j))>=alpha .or. & 00493 !! (present(aggr_fine).and.dabs(A%val(j))>=beta))) then 00494 ! 00495 if (aggrsize(colr)<maxasizelargest.and. & 00496 dabs(A%val(j))>=beta ) then 00497 aggrsize(colr)=aggrsize(colr)+1 00498 structnodes(aggrsize(colr),colr)=structnodes(k,i) 00499 aggrnum(structnodes(k,i))=colr 00500 structnodes(k,i)=-1 00501 nleft=nleft-1 00502 reduced=.true. 00503 !print *,i,' ########## sold ',structnodes(aggrsize(colr),colr), & 00504 ! ' to:',colr,'##########' 00505 exit looking ! this node is sold 00506 endif 00507 endif 00508 endif 00509 enddo looking 00510 endif 00511 enddo 00512 enddo ! while 00513 if (nleft>0.and.nleft<aggrsize(i)) then ! compress the list 00514 k=0 ! of left behind nodes 00515 do j=1,aggrsize(i) 00516 if (structnodes(j,i)>0) then 00517 k=k+1 00518 if (k>0.and.k<j) then 00519 structnodes(k,i)=structnodes(j,i) 00520 endif 00521 endif 00522 enddo 00523 aggrsize(i)=nleft 00524 print *,' ========== aggregate ',i,' remained as small as ',nleft,'@@@@@' 00525 elseif (nleft==0) then ! the aggregate got removed 00526 print *,' ========== aggregate ',i,' got removed node by node ============' 00527 noccupied=noccupied+1 00528 aggrsize(i)=0 00529 nagrs_new=nagrs_new-1 00530 endif 00531 endif 00532 endif 00533 enddo ! } 00534 naggregatednodes=sum(aggrsize(1:nagrs)) 00535 full_nagrs_new=nagrs_new 00536 !write(stream,*)'aaaa nagrs=',nagrs_new 00537 !write(stream,*)'aaaa aggrnum=',aggrnum 00538 ! compress the local numbers if needed 00539 cnt=maxval(aggrnum(1:n)) 00540 if (cnt>nagrs_new) then 00541 allocate(stat(cnt)) 00542 stat=0 00543 col=0 00544 do i=1,n 00545 if (aggrnum(i)>0) then 00546 if (stat(aggrnum(i))==0) then 00547 col=col+1 00548 stat(aggrnum(i))=col 00549 thiscol=col 00550 else 00551 thiscol=stat(aggrnum(i)) 00552 endif 00553 aggrnum(i)=thiscol 00554 endif 00555 enddo 00556 deallocate(stat) 00557 ! else 00558 !col=nagrs_new 00559 endif 00560 !write(stream,*)'bbbb nagrs=',col 00561 !write(stream,*)'bbbb aggrnum=',aggrnum 00562 if (n==naggregatednodes) then 00563 aggrarefull=.true. 00564 fullaggrnum=aggrnum(1:A%nrows) 00565 full_nagrs_new=maxval(fullaggrnum) 00566 else 00567 aggrarefull=.false. 00568 do i=1,n 00569 if (aggrnum(i)>0) then 00570 fullaggrnum(i)=aggrnum(i) 00571 else 00572 full_nagrs_new=full_nagrs_new+1 00573 fullaggrnum(i)=full_nagrs_new 00574 endif 00575 enddo 00576 endif 00577 call Form_Aggr(aggr=aggr%inner, & 00578 nagrs=nagrs_new, & 00579 n=n, & 00580 radius=neighood, & 00581 nisolated=n-naggregatednodes, & 00582 aggrnum=aggrnum) 00583 deallocate(connweightsums) ! weight sums to each colour! 00584 deallocate(colsaround) ! lists the colors 00585 deallocate(structnodes) 00586 deallocate(aggrsize) ! the initial aggr sizes 00587 nagrs=nagrs_new 00588 write(stream,*)'# too small aggrs:',ntoosmall,' #eaten:',neaten, & 00589 ' #occupied:',noccupied, & 00590 ' # remaining:',ntoosmall-neaten-noccupied 00591 endif !} 00592 call Form_Aggr(aggr=aggr%full, & 00593 nagrs=full_nagrs_new, & 00594 n=n, & 00595 radius=neighood, & 00596 nisolated=nisolated, & 00597 aggrnum=fullaggrnum) 00598 deallocate(fullaggrnum) 00599 deallocate(aggrnum) 00600 00601 ! if (plot==1.or.plot==3) then 00602 ! if (numprocs>1) then 00603 ! if (ismaster()) then 00604 ! allocate(aggrnum(M%ngf)) 00605 ! allocate(owner(M%ngf)) 00606 ! end if 00607 ! call Integer_Vect_Gather(aggr%inner%num,aggrnum,M,owner) 00608 ! if (ismaster()) then 00609 ! call color_print_aggrs(M%ngf,aggrnum,overwrite=.false.,owner=owner) 00610 ! deallocate(owner,aggrnum) 00611 ! endif 00612 ! else 00613 ! if (.not.present(aggr_fine)) then 00614 ! if (plot==3) then 00615 ! call color_print_aggrs(A%nrows,aggr%inner%num,overwrite=.true.) 00616 ! else 00617 ! write(stream,*)' fine aggregates:' 00618 ! call color_print_aggrs(A%nrows,aggr%inner%num) 00619 ! if (.not.aggrarefull) then 00620 ! write(stream,*)' FULL fine aggregates:' 00621 ! call color_print_aggrs(A%nrows,aggr%full%num) 00622 ! endif 00623 ! endif 00624 ! else 00625 ! if (plot==3) then 00626 ! call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%full%num,aggr%inner%num,overwrite=.true.) 00627 ! else 00628 ! write(stream,*)' coarse aggregates:' 00629 ! call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%inner%num,aggr%inner%num) 00630 ! if (.not.aggrarefull) then 00631 ! write(stream,*)' FULL coarse aggregates:' 00632 ! call color_print_aggrs(size(aggr_fine%full%nodes),aggr_fine%full%num,aggr%full%num) 00633 ! endif 00634 ! endif 00635 ! endif 00636 ! endif 00637 ! endif 00638 ! if (plot==3) then 00639 ! deallocate(moviecols) 00640 ! endif 00641 end subroutine SpMtx_aggregate 00642 00643 subroutine SpMtx_exchange_strong(A,A_ghost,M) 00644 type(SpMtx), intent(in) :: A 00645 type(SpMtx), intent(inout) :: A_ghost 00646 type(Mesh), intent(in) :: M 00647 00648 call exchange_strong() 00649 00650 contains 00651 function calcBufferSize(ninds) result(bufferSize) 00652 integer, intent(in) :: ninds 00653 integer :: bufferSize 00654 integer :: size, ierr 00655 call MPI_Pack_size(ninds, MPI_INTEGER, MPI_COMM_WORLD, size, ierr) 00656 bufferSize = 2*size 00657 end function calcBufferSize 00658 00659 function calcBufferSize2(ninds) result(bufferSize) 00660 integer, intent(in) :: ninds 00661 integer :: bufferSize 00662 integer :: size, ierr 00663 call MPI_Pack_size(ninds, MPI_LOGICAL, MPI_COMM_WORLD, size, ierr) 00664 bufferSize = size 00665 end function calcBufferSize2 00666 00667 subroutine exchange_strong() 00668 integer :: k, k2, neigh, ninds, bufsize, bufpos, ptn 00669 type(indlist), allocatable :: strong_sends(:), strong_recvs(:) 00670 logical, allocatable :: strongval_recvs(:) 00671 type Buffer 00672 character, pointer :: data(:) 00673 end type Buffer 00674 type(Buffer),allocatable :: outbuffers(:) 00675 integer, allocatable :: outreqs(:), outstatuses(:,:), indi(:), indj(:) 00676 integer :: status(MPI_STATUS_SIZE), ierr 00677 character, allocatable :: inbuffer(:) 00678 logical, allocatable :: strong(:) 00679 integer :: i,j,nnz 00680 00681 allocate(strong_sends(M%nnghbrs)) 00682 00683 ! report to the neighbours which elements we need 00684 !! how many elements 00685 strong_sends%ninds = 0 00686 do k=1,A_ghost%nnz 00687 ptn = M%eptnmap(M%lg_fmap(A_ghost%indi(k))) 00688 if (ptn/=myrank+1) then 00689 ! find neighbour number 00690 do i=1,M%nnghbrs 00691 if (M%nghbrs(i)+1==ptn) then 00692 neigh = i 00693 exit 00694 end if 00695 end do 00696 strong_sends(neigh)%ninds = strong_sends(neigh)%ninds+1 00697 end if 00698 end do 00699 !! collect elements 00700 !!! prepare indices 00701 do neigh=1,M%nnghbrs 00702 allocate(strong_sends(neigh)%inds(strong_sends(neigh)%ninds)) 00703 end do 00704 strong_sends%ninds = 0 ! reset 00705 do k=1,A_ghost%nnz 00706 ptn = M%eptnmap(M%lg_fmap(A_ghost%indi(k))) 00707 if (ptn/=myrank+1) then 00708 ! find neighbour number 00709 do i=1,M%nnghbrs 00710 if (M%nghbrs(i)+1==ptn) then 00711 neigh = i 00712 exit 00713 end if 00714 end do 00715 ninds = strong_sends(neigh)%ninds 00716 strong_sends(neigh)%ninds = ninds+1 00717 strong_sends(neigh)%inds(ninds+1) = k 00718 end if 00719 end do 00720 00721 ! gather number of matrix elements 00722 allocate(strong_recvs(M%nnghbrs)) 00723 do i=1,M%nnghbrs 00724 neigh = M%nghbrs(i) 00725 call MPI_Send(strong_sends(i)%ninds, 1, MPI_INTEGER, neigh, & 00726 TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, ierr) 00727 end do 00728 do i=1,M%nnghbrs 00729 neigh = M%nghbrs(i) 00730 call MPI_Recv(strong_recvs(i)%ninds, 1, MPI_INTEGER, neigh, & 00731 TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr) 00732 end do 00733 00734 !!! pack and send index data 00735 allocate(outbuffers(M%nnghbrs)) 00736 allocate(outreqs(M%nnghbrs)) 00737 do neigh=1,M%nnghbrs 00738 bufsize = calcBufferSize(strong_sends(neigh)%ninds) 00739 allocate(outbuffers(neigh)%data(bufsize)) 00740 bufpos = 0 00741 call MPI_Pack(M%lg_fmap(A_ghost%indi(strong_sends(neigh)%inds)), strong_sends(neigh)%ninds, MPI_INTEGER,& 00742 outbuffers(neigh)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr) 00743 call MPI_Pack(M%lg_fmap(A_ghost%indj(strong_sends(neigh)%inds)), strong_sends(neigh)%ninds, MPI_INTEGER,& 00744 outbuffers(neigh)%data, bufsize, bufpos, MPI_COMM_WORLD, ierr) 00745 call MPI_ISend(outbuffers(neigh)%data, bufsize, MPI_CHARACTER, M%nghbrs(neigh), TAG_EXCHANGE_STRONG,& 00746 MPI_COMM_WORLD, outreqs(neigh), ierr) 00747 end do 00748 00749 !!! recv index data 00750 allocate(inbuffer(calcBufferSize(maxval(strong_recvs%ninds)))) 00751 nnz = sum(strong_recvs%ninds) 00752 allocate(indi(nnz),indj(nnz)) 00753 nnz = 0 00754 do i=1,M%nnghbrs 00755 bufsize = calcBufferSize(strong_recvs(i)%ninds) 00756 call MPI_Recv(inbuffer, bufsize, MPI_CHARACTER, M%nghbrs(i),& 00757 TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr) 00758 ninds = strong_recvs(i)%ninds 00759 bufpos = 0 00760 if (ninds==0) cycle 00761 call MPI_Unpack(inbuffer, bufsize, bufpos, indi(nnz+1), ninds,& 00762 MPI_INTEGER, MPI_COMM_WORLD, ierr) 00763 if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed") 00764 call MPI_Unpack(inbuffer, bufsize, bufpos, indj(nnz+1), ninds,& 00765 MPI_INTEGER, MPI_COMM_WORLD, ierr) 00766 if (ierr/=0) call DOUG_abort("MPI UnPack of matrix elements failed") 00767 nnz = nnz+ninds 00768 end do 00769 00770 ! find strong values 00771 allocate(strong(nnz)) 00772 do k=1,nnz 00773 i = M%gl_fmap(indi(k)) 00774 j = M%gl_fmap(indj(k)) 00775 k2 = SpMtx_findElem(A, i, j) 00776 if(k2<=0) call DOUG_abort("Matrix element not found during 'strong' exchange") 00777 strong(k) = A%strong(k2) 00778 end do 00779 00780 ! wait for all data to be sent 00781 allocate(outstatuses(MPI_STATUS_SIZE, M%nnghbrs)) 00782 call MPI_Waitall(M%nnghbrs, outreqs, outstatuses, ierr) 00783 00784 ! send strong values back 00785 nnz = 0 00786 bufpos = 0 00787 do i=1,M%nnghbrs 00788 ninds = strong_recvs(i)%ninds 00789 if (ninds==0) cycle 00790 call MPI_ISend(strong(nnz+1), ninds, MPI_LOGICAL, M%nghbrs(i), & 00791 TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, outreqs(i), ierr) 00792 nnz = nnz+ninds 00793 end do 00794 00795 ! receive strong values (requested value indices are in strong_sends) 00796 allocate(A_ghost%strong(A_ghost%nnz)) 00797 nnz = sum(strong_sends%ninds) 00798 allocate(strongval_recvs(nnz)) 00799 nnz = 0 00800 do i=1,M%nnghbrs 00801 ninds = strong_sends(i)%ninds 00802 if (ninds==0) cycle 00803 call MPI_Recv(strongval_recvs(nnz+1), ninds, MPI_LOGICAL,M%nghbrs(i),& 00804 TAG_EXCHANGE_STRONG, MPI_COMM_WORLD, status, ierr) 00805 ! overwrite local strong value with remote 00806 !write (stream,*) "----", strong_sends(i)%ninds, ninds, strong_sends(i)%inds, strongval_recvs(nnz+1:nnz+ninds) 00807 do k=1,ninds 00808 k2 = strong_sends(i)%inds(k) 00809 !write(stream,*) "--k2", k2, strongval_recvs(nnz+k), size(A_ghost%strong) 00810 A_ghost%strong(k2) = strongval_recvs(nnz+k) 00811 end do 00812 00813 nnz = nnz+ninds 00814 end do 00815 00816 ! free memory 00817 do neigh=1,M%nnghbrs 00818 deallocate(outbuffers(neigh)%data) 00819 deallocate(strong_sends(neigh)%inds) 00820 end do 00821 deallocate(outbuffers) 00822 deallocate(strong_sends, strong_recvs) 00823 end subroutine exchange_strong 00824 end subroutine SpMtx_exchange_strong 00825 00826 !------------------------------------------------------ 00827 ! Finding strong connections in matrix 00828 !------------------------------------------------------ 00829 subroutine SpMtx_find_strong(A,alpha,A_ghost,symmetrise) 00830 Implicit None 00831 Type(SpMtx),intent(in out) :: A 00832 float(kind=rk), intent(in) :: alpha 00833 Type(SpMtx),intent(in out),optional :: A_ghost 00834 logical,optional :: symmetrise 00835 ! local: 00836 integer :: i,j,k,start,ending,nnz,ndiags 00837 logical :: did_scale 00838 logical :: simple=.false., symm=.false. 00839 float(kind=rk) :: maxndiag,aa 00840 did_scale=.false. 00841 if (A%scaling==D_SpMtx_SCALE_NO.or.A%scaling==D_SpMtx_SCALE_UNDEF) then 00842 call SpMtx_scale(A,A_ghost) 00843 did_scale=.true. 00844 endif 00845 if (A%mtx_bbe(2,2)>0) then 00846 nnz=A%mtx_bbe(2,2) 00847 else 00848 nnz=A%nnz 00849 endif 00850 00851 if (simple) then 00852 if (.not.associated(A%strong)) allocate(A%strong(nnz)) 00853 do i=1,nnz 00854 if (abs(A%val(i))>=alpha) then 00855 A%strong(i)=.true. 00856 else 00857 A%strong(i)=.false. 00858 endif 00859 enddo 00860 else ! not the simple case: 00861 ndiags=max(A%nrows,A%ncols) 00862 call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.false.) 00863 if (.not.associated(A%strong)) allocate(A%strong(nnz)) 00864 do i=1,A%nrows 00865 start=A%M_bound(i) 00866 ending=A%M_bound(i+1)-1 00867 maxndiag=-1.0e15 00868 do j=start,ending 00869 if (A%indj(j)/=i) then ! not on diagonal 00870 aa=abs(A%val(j)) 00871 if (maxndiag<aa) then 00872 maxndiag=aa 00873 endif 00874 endif 00875 enddo 00876 maxndiag=maxndiag*alpha 00877 do j=start,ending 00878 aa=abs(A%val(j)) 00879 if (A%indj(j)/=i) then ! not on diagonal 00880 if (aa>maxndiag) then 00881 A%strong(j)=.true. 00882 else 00883 A%strong(j)=.false. 00884 endif 00885 else 00886 if (A%diag(i)>maxndiag) then 00887 A%strong(j)=.true. 00888 else 00889 A%strong(j)=.false. 00890 endif 00891 !write(stream,*)'diag at i is:',i,A%diag(i),A%strong(j),maxndiag 00892 endif 00893 enddo 00894 enddo 00895 endif 00896 if (present(symmetrise)) then 00897 symm=symmetrise 00898 endif 00899 call SpMtx_symm_strong(A,A_ghost,symm) 00900 end subroutine SpMtx_find_strong 00901 00902 subroutine SpMtx_symm_strong(A,A_ghost,symm) 00903 Type(SpMtx),intent(in out) :: A 00904 Type(SpMtx),intent(in out),optional :: A_ghost 00905 logical,intent(in) :: symm 00906 logical :: mirror2ghost=.true. 00907 00908 integer :: i,k,nnz 00909 00910 if (A%mtx_bbe(2,2)>0) then 00911 nnz=A%mtx_bbe(2,2) 00912 else 00913 nnz=A%nnz 00914 endif 00915 00916 if (symm) then 00917 do i=1,nnz 00918 if (A%arrange_type == D_SpMtx_ARRNG_ROWS) then 00919 k=SpMtx_findElem(A,A%indi(i),A%indj(i)) 00920 else 00921 k=SpMtx_findElem(A,A%indj(i),A%indi(i)) 00922 endif 00923 if (k>0) then 00924 if (A%strong(i).and..not.A%strong(k)) then 00925 A%strong(k)=.true. 00926 write(stream,*)'symmetrising to strong:',A%indi(i),A%indj(i) 00927 elseif (A%strong(k).and..not.A%strong(i)) then 00928 A%strong(i)=.true. 00929 write(stream,*)'symmetrising to strong:',A%indi(i),A%indj(i) 00930 endif 00931 else 00932 write(stream,*) 'Warning: matrix does not have symmetric structure!' 00933 endif 00934 enddo 00935 endif 00936 if (mirror2ghost) then 00937 if (present(A_ghost).and.associated(A_ghost%indi)) then 00938 do i=1,A_ghost%nnz 00939 if(A_ghost%indi(i)/=A_ghost%indj(i)) then 00940 k=SpMtx_findElem(A,A_ghost%indi(i),A_ghost%indj(i)) 00941 if (k>0) then 00942 A_ghost%strong(i)=A%strong(k) 00943 endif 00944 endif 00945 enddo 00946 endif 00947 endif 00948 end subroutine SpMtx_symm_strong 00949 00950 !------------------------------------------------------ 00951 End Module SpMtx_aggregation 00952 !------------------------------------------------------ 00953
1.7.3-20110217