|
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 Aggregate_mod 00024 use RealKind 00025 use globals 00026 use Mesh_class 00027 00028 Implicit None 00029 00030 #include<doug_config.h> 00031 00032 ! "on-the-fly" real/complex picking 00033 #ifdef D_COMPLEX 00034 #define float complex 00035 #else 00036 #define float real 00037 #endif 00038 00039 Type Aggrs 00040 integer :: nagr !< #aggregates 00041 integer :: radius !< radius of aggregation 00042 integer :: nisolated !< # of isolated nodes 00043 integer, dimension(:),pointer :: num !< aggregate # for each node 00044 integer, dimension(:),pointer :: starts,nodes !< compressed storage 00045 end Type Aggrs !Aggrs 00046 00048 type AggrInfo 00049 type(Aggrs) :: inner !< aggregates (on all inner freedoms) 00050 type(Aggrs) :: full !< aggr with holes painted over 00051 end type AggrInfo 00052 00053 logical :: debu = .false. 00054 CONTAINS 00055 00057 function Aggrs_New() result (aggr) 00058 type(Aggrs) :: aggr 00059 00060 aggr%nagr=0 00061 aggr%radius=0 00062 aggr%nisolated=0 00063 nullify(aggr%num) 00064 nullify(aggr%starts) 00065 nullify(aggr%nodes) 00066 end function Aggrs_New 00067 00068 function AggrInfo_New() result (aggr) 00069 type(AggrInfo) :: aggr 00070 00071 aggr%inner = Aggrs_New() 00072 aggr%full = Aggrs_New() 00073 end function AggrInfo_New 00074 00075 subroutine AggrInfo_Destroy(aggr) 00076 type(AggrInfo) :: aggr 00077 00078 call Destruct_Aggrs(aggr%inner) 00079 call Destruct_Aggrs(aggr%full) 00080 end subroutine AggrInfo_Destroy 00081 00082 subroutine Form_Aggr(aggr,nagrs,n,radius,nisolated,aggrnum) 00083 Implicit None 00084 type(Aggrs) :: aggr 00085 integer, intent(in) :: nagrs ! number of aggregates 00086 integer, intent(in) :: n ! number of unknowns 00087 integer, intent(in) :: radius ! aggr radius (or neighood) 00088 integer, intent(in) :: nisolated ! number of isolated nodes 00089 integer, dimension(:), intent(in) :: aggrnum 00090 00091 integer,dimension(:),pointer :: aggrstarts,stat,aggrnodes 00092 integer :: i,j 00093 allocate(aggrstarts(nagrs+1)) 00094 allocate(stat(nagrs)) ! stat gets different meaning here... 00095 stat=0 00096 do i=1,n ! find the #nodes for each color 00097 j=aggrnum(i) 00098 if (j>0) then 00099 stat(j)=stat(j)+1 00100 endif 00101 enddo 00102 aggrstarts(1)=1 00103 do i=1,nagrs 00104 aggrstarts(i+1)=aggrstarts(i)+stat(i) 00105 stat(i)=aggrstarts(i) ! shows the place to fill the nodes 00106 enddo 00107 allocate(aggrnodes(aggrstarts(nagrs+1)-1)) 00108 do i=1,n ! put the node#-s in 00109 j=aggrnum(i) 00110 if (j>0) then 00111 aggrnodes(stat(j))=i 00112 stat(j)=stat(j)+1 00113 endif 00114 enddo 00115 if (sctls%verbose>=5) then 00116 do i=1,nagrs 00117 write(stream,*) & 00118 'aggregate',i,':',(aggrnodes(j),j=aggrstarts(i),aggrstarts(i+1)-1) 00119 enddo 00120 endif 00121 call Construct_Aggrs(aggr,nagrs,n,radius,nisolated,aggrnum,aggrstarts,aggrnodes) 00122 deallocate(aggrnodes) 00123 deallocate(stat) 00124 deallocate(aggrstarts) 00125 end subroutine Form_Aggr 00126 00127 subroutine Construct_Aggrs(aggr,nagr,n,radius,nisolated,num,starts,nodes) 00128 Implicit None 00129 type(Aggrs) :: aggr 00130 integer, intent(in) :: nagr ! number of aggregates 00131 integer, intent(in) :: n ! number of unknowns 00132 integer, intent(in) :: radius ! aggr radius (or neighood) 00133 integer, intent(in) :: nisolated ! number of isolated nodes 00134 integer, dimension(:), intent(in) :: num 00135 integer, dimension(:), intent(in) :: starts 00136 integer, dimension(:), intent(in) :: nodes 00137 00138 aggr%nagr=nagr 00139 aggr%radius=radius 00140 aggr%nisolated=nisolated 00141 allocate(aggr%num(n)) 00142 allocate(aggr%starts(nagr+1)) 00143 allocate(aggr%nodes(starts(nagr+1)-1)) 00144 aggr%num=num(1:n) 00145 aggr%starts=starts 00146 aggr%nodes=nodes 00147 end subroutine Construct_Aggrs 00148 00149 subroutine Destruct_Aggrs(aggr) 00150 Implicit None 00151 type(Aggrs) :: aggr 00152 00153 if (associated(aggr%num)) deallocate(aggr%num) 00154 if (associated(aggr%starts)) deallocate(aggr%starts) 00155 if (associated(aggr%nodes)) deallocate(aggr%nodes) 00156 aggr%nagr=0 00157 aggr%radius=0 00158 aggr%nisolated=0 00159 end subroutine Destruct_Aggrs 00160 00162 subroutine Get_aggregate_nodes(cAggr,cAggrs,fAggrs,maxnodes, nodes, nnodes) 00163 type(Aggrs),intent(in) :: cAggrs 00164 type(Aggrs),intent(in) :: fAggrs 00165 integer,intent(in) :: cAggr 00166 integer,intent(in) :: maxnodes 00167 integer,intent(inout) :: nodes(:) 00168 integer,intent(out) :: nnodes 00169 00170 integer,dimension(:),allocatable :: nodeLoc ! node location in nodes array 00171 integer :: ifAggr, fAggr, inode, node 00172 00173 allocate(nodeLoc(maxnodes)) 00174 nodeLoc = 0 00175 nnodes = 0 00176 00177 nodeLoc=0 00178 nnodes=0 00179 do ifAggr=cAggrs%starts(cAggr),cAggrs%starts(cAggr+1)-1 ! look through fine agrs. 00180 fAggr=cAggrs%nodes(ifAggr) ! fine aggregate numbers 00181 do inode=fAggrs%starts(fAggr),fAggrs%starts(fAggr+1)-1 ! loop over nodes in aggr. 00182 node=fAggrs%nodes(inode) ! node number 00183 if (nodeLoc(node)==0) then 00184 nnodes=nnodes+1 00185 nodeLoc(node)=nnodes 00186 nodes(nnodes)=node 00187 endif 00188 enddo 00189 enddo 00190 end subroutine Get_aggregate_nodes 00191 00192 function node_neighood_fits(innode,neighood,nneigs,nodes, & 00193 stat,rowstart,colnrs) result(ok) 00194 Implicit None 00195 logical :: ok 00196 integer, intent(in) :: innode ! the node who's neighood is being built 00197 integer, intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00198 integer, dimension(:), pointer :: stat 00199 integer, dimension(:), pointer :: rowstart,colnrs 00200 integer, intent(in out) :: nneigs 00201 integer, dimension(:) :: nodes 00202 nneigs=0 00203 stat(innode)=1 ! mark it as if on layer 1 (although it's layer is 0) 00204 if (marking_neighs_ok(0,innode,neighood,stat, & ! ok 00205 nneigs,nodes,rowstart,colnrs)) then 00206 ok=.true. 00207 else ! not ok 00208 stat(innode)=0 ! reverse to 0 00209 ok=.false. 00210 endif 00211 end function node_neighood_fits 00212 00213 recursive function marking_neighs_ok(layer,innode,neighood,stat, & 00214 nneigs,nodes,rowstart,colnrs) result(ok) 00215 00216 Implicit None 00217 logical :: ok 00218 integer,intent(in) :: layer ! the layer # of innode 00219 integer, intent(in) :: innode ! the node who's neighood is being built 00220 integer, intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00221 integer, dimension(:), pointer :: stat 00222 integer, dimension(:), pointer :: rowstart,colnrs 00223 integer, intent(in out) :: nneigs 00224 integer, dimension(:) :: nodes 00225 integer :: nneigs_in,nxtlayer ! nxtlayer -- the next layer to look at 00226 integer :: my_nneigs,radius2 00227 integer :: i,j,rs,re,cn 00228 00229 ! we assume, symmetric structure 00230 if (neighood==1) then 00231 radius2=2 00232 else 00233 !radius2=2*neighood 00234 !radius2=2*neighood-1 00235 radius2=neighood+1 00236 endif 00237 nneigs_in=nneigs 00238 nxtlayer=layer+1 00239 rs=rowstart(innode) 00240 re=rowstart(innode+1)-1 00241 ! Let's look first, if any of innode neighbours fails to fit 00242 if (nxtlayer<=neighood) then ! still r-neighood 00243 do i=rs,re ! check the neighbours: 00244 cn=colnrs(i) 00245 if (stat(cn)==-2) then ! this r-neighood will not fit in! 00246 ! reverse all stat nneigs 00247 do j=1,nneigs 00248 if (nodes(j)<0) then 00249 stat(-nodes(j))=-1 ! node in the shadow 00250 else 00251 stat(nodes(j))=0 00252 endif 00253 enddo 00254 nneigs=0 00255 ok=.false. 00256 return 00257 endif 00258 enddo 00259 endif 00260 do i=rs,re ! put the neighs into the list if needed 00261 cn=colnrs(i) 00262 if (stat(cn)==-1) then !node cn was in shade and not included yet 00263 stat(cn)=nxtlayer 00264 nneigs=nneigs+1 00265 nodes(nneigs)=-cn ! "-" marking previous shade 00266 else if (stat(cn)==0) then ! node cn to be included 00267 stat(cn)=nxtlayer 00268 nneigs=nneigs+1 00269 nodes(nneigs)=cn 00270 endif 00271 enddo 00272 my_nneigs=nneigs 00273 !if (nxtlayer<2*neighood) then 00274 if (nxtlayer<radius2) then ! if there are more layers to process 00275 do i=nneigs_in+1,my_nneigs ! recursively check the neighbours 00276 if (.not.marking_neighs_ok(nxtlayer,abs(nodes(i)),neighood,stat, & 00277 nneigs,nodes,rowstart,colnrs)) then 00278 ok=.false. 00279 return 00280 endif 00281 enddo 00282 endif 00283 ok=.true. 00284 return 00285 end function marking_neighs_ok 00286 00287 function lets_colour(innode,neighood,minasize,maxasize,nneigs,nodes, & 00288 stat,distance,rowstart,colnrs) result(ok) 00289 !ok = 0 => nothing found at all (nneigs==0) 00290 ! 1 => neighood not achieved 00291 ! 2 => neighood achieved but not the double-neighood 00292 ! 3 => double-neighood achieved 00293 Implicit None 00294 integer :: ok 00295 integer,intent(in) :: innode ! the node who's neighood is being built 00296 integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00297 integer,intent(in) :: minasize,maxasize ! aggregate limits 00298 integer,intent(inout) :: nneigs 00299 integer,dimension(:), pointer :: stat 00300 integer,dimension(:), pointer :: rowstart,colnrs,distance 00301 integer,dimension(:) :: nodes 00302 integer :: ninnodes,virtaggrsize 00303 integer,dimension(:), pointer :: inlayer,newlayer,roundernodes 00304 nneigs=0 00305 virtaggrsize=2*maxasize**2 00306 stat(innode)=1 ! mark it as if on layer 1 (although it's layer is 0) 00307 distance(innode)=0 ! seed node distance is 0 00308 allocate(inlayer(virtaggrsize)) 00309 allocate(newlayer(virtaggrsize)) 00310 allocate(roundernodes(maxasize)) 00311 nneigs=0 00312 ninnodes=1 00313 inlayer(1)=innode 00314 ok=can_add_roundlayer(0,ninnodes,inlayer,newlayer, & 00315 roundernodes,neighood,minasize,maxasize,stat,distance, & 00316 nneigs,nodes,rowstart,colnrs) 00317 deallocate(roundernodes) 00318 deallocate(newlayer) 00319 deallocate(inlayer) 00320 return 00321 end function lets_colour 00322 00323 recursive function can_add_roundlayer(layer,ninnodes, & 00324 inlayer,newlayer,roundernodes,neighood,& 00325 minasize,maxasize,stat,distance, & 00326 nneigs,nodes,rowstart,colnrs) result(ok) 00327 use globals 00328 Implicit None 00329 integer :: ok 00330 integer,intent(in) :: layer ! the layer # of innode 00331 integer,intent(inout) :: ninnodes ! number of nodes on the last layer 00332 integer,dimension(:),pointer :: inlayer ! nodes on the last outer layer 00333 integer,dimension(:),pointer :: newlayer ! storage for new layer 00334 integer,dimension(:),pointer :: roundernodes ! storage for rounding 00335 integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00336 integer,intent(in) :: minasize,maxasize ! aggregate limits 00337 integer,dimension(:),pointer :: stat 00338 integer,dimension(:),pointer :: rowstart,colnrs,distance 00339 integer,intent(inout) :: nneigs 00340 integer,dimension(:) :: nodes 00341 integer :: nxtlayer ! nxtlayer -- the next layer to look at 00342 integer :: i,j,k,d,rs,re,cn,layernode,ntoadd,nrounders,mind,maxd 00343 logical :: layerfits 00344 !logical :: checklayerfit=.false. 00345 logical :: checklayerfit=.true. 00346 !logical :: dorounding=.false. 00347 logical :: dorounding=.true. 00348 integer,save :: allmind=D_MAXINT 00349 00350 ! we assume, symmetric structure 00351 nxtlayer=layer+1 00352 layerfits=.true. 00353 ntoadd=0 00354 if (layerfits.and.layer<=2*neighood) then 00355 aaa:do k=1,ninnodes 00356 layernode=inlayer(k) 00357 rs=rowstart(layernode) 00358 re=rowstart(layernode+1)-1 00359 do i=rs,re ! put the neighs into the list if needed 00360 cn=colnrs(i) 00361 if (stat(cn)<=0) then !node not marked yet 00362 ! wheather this node has been seen from inlayer already?: 00363 if (stat(cn)==0) then ! first time 00364 !if (layer<=neighood.and.nneigs+ntoadd+1>maxasize) then 00365 ! layerfits=.false. 00366 ! exit aaa 00367 !endif 00368 ntoadd=ntoadd+1 00369 newlayer(ntoadd)=cn 00370 distance(cn)=distance(layernode)+1 00371 endif 00372 stat(cn)=stat(cn)-1 ! We use stat here for finding weights for 00373 ! rounding step... 00374 endif 00375 enddo 00376 enddo aaa 00377 endif 00378 if (layerfits.and.layer<=neighood) then !{ 00379 ! do the rounding step: 00380 nrounders=0 00381 ! find the min/max weight for the newlayer: 00382 mind=D_MAXINT 00383 maxd=0 00384 do k=1,ntoadd 00385 d=-stat(newlayer(k)) 00386 mind=min(d,mind) 00387 maxd=max(d,maxd) 00388 enddo 00389 if (mind<allmind) then 00390 allmind=mind 00391 endif 00392 if (mind<maxd.or.mind>allmind) then 00393 bbb:do k=1,ntoadd 00394 if (stat(newlayer(k))==-maxd) then ! this is a rounder node 00395 nrounders=nrounders+1 00396 if (checklayerfit.and.nneigs+ninnodes+1>maxasize) then 00397 layerfits=.false. 00398 exit bbb 00399 endif 00400 ninnodes=ninnodes+1 00401 inlayer(ninnodes)=newlayer(k) ! add this rounder node 00402 roundernodes(nrounders)=newlayer(k) 00403 newlayer(k)=-1 ! this is taken out from the new layer 00404 endif 00405 enddo bbb 00406 endif 00407 ! if rounding of the incoming layer succeeded, we can add this to the 00408 ! aggregate... 00409 ! ... done later 00410 if (.not.layerfits) then 00411 nrounders=0 00412 endif 00413 ! now still add neighbours of rounded nodes to newlayer: 00414 ! first take out holes from the newlayer 00415 !if (nrounders>0) then 00416 j=0 00417 do i=1,ntoadd 00418 if (newlayer(i)/=-1) then 00419 j=j+1 00420 if (j<i) then 00421 newlayer(j)=newlayer(i) 00422 endif 00423 endif 00424 enddo 00425 ntoadd=j 00426 !endif 00427 ! add rounder's neighbours to the newlayer: 00428 ccc:do k=1,nrounders 00429 layernode=roundernodes(k) 00430 rs=rowstart(layernode) 00431 re=rowstart(layernode+1)-1 00432 do i=rs,re ! put the neighs into the list if needed 00433 cn=colnrs(i) 00434 if (stat(cn)==0) then !node not seen yet 00435 !if (nneigs+ntoadd+1>maxasize) then 00436 ! layerfits=.false. 00437 ! exit ccc 00438 !endif 00439 ntoadd=ntoadd+1 00440 newlayer(ntoadd)=cn 00441 distance(cn)=distance(layernode)+1 00442 endif 00443 enddo 00444 enddo ccc 00445 endif !(layer<=neighood) } 00446 if (checklayerfit.and.layer<=neighood.and.nneigs+ninnodes>maxasize) then ! this could 00447 ! happen, if no rounding nodes got added... 00448 layerfits=.false. 00449 endif 00450 if (layerfits) then ! { add the layer 00451 do k=1,ninnodes 00452 stat(inlayer(k))=layer+1 00453 nneigs=nneigs+1 00454 nodes(nneigs)=inlayer(k) 00455 enddo 00456 if (layer<=2*neighood) then ! recurse to the next layer 00457 do k=1,ntoadd 00458 stat(newlayer(k))=D_PENDING-1 ! just temporary pos. marker here... 00459 enddo 00460 ninnodes=ntoadd 00461 inlayer(1:ninnodes)=newlayer(1:ntoadd) ! (this does not incl. rounders) 00462 ok=can_add_roundlayer(nxtlayer,ninnodes, & 00463 inlayer,newlayer,roundernodes,neighood, & 00464 minasize,maxasize,stat,distance, & 00465 nneigs,nodes,rowstart,colnrs) 00466 else 00467 ok=3 00468 endif 00469 return 00470 else ! }{ clean up the possible aggregate nodes and also on pending layers: 00471 do k=1,ninnodes 00472 distance(inlayer(k))=0 00473 enddo 00474 do k=1,ntoadd 00475 distance(newlayer(k))=0 00476 enddo 00477 do k=1,nneigs 00478 distance(nodes(k))=0 00479 enddo 00480 ! add the found nodes to the inlayer to be able to give them appropriate 00481 ! failing structure number outside: 00482 do k=1,ninnodes 00483 nneigs=nneigs+1 00484 nodes(nneigs)=inlayer(k) 00485 enddo 00486 if (nneigs>0) then 00487 if (layer>neighood) then 00488 ok=2 00489 else 00490 ok=1 00491 endif 00492 else 00493 ok=0 00494 endif 00495 return 00496 endif ! } 00497 end function can_add_roundlayer 00498 00499 function lets_colour3(innode,neighood,minasize,maxasize,nneigs,nodes, & 00500 stat,distance,rowstart,colnrs) result(ok) 00501 !ok = 0 => minasize not achieved 00502 ! 1 => minasize achieved but double-neighood not achieved 00503 ! 2 => minasize achieved and double-neighood achieved 00504 Implicit None 00505 integer :: ok 00506 integer,intent(in) :: innode ! the node who's neighood is being built 00507 integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00508 integer,intent(in) :: minasize,maxasize ! aggregate limits 00509 integer,intent(inout) :: nneigs 00510 integer,dimension(:), pointer :: stat 00511 integer,dimension(:), pointer :: rowstart,colnrs,distance 00512 integer,dimension(:) :: nodes 00513 integer :: ninnodes,virtaggrsize 00514 integer,dimension(:), pointer :: inlayer,newlayer,roundernodes 00515 nneigs=0 00516 virtaggrsize=2*maxasize**2 00517 stat(innode)=1 ! mark it as if on layer 1 (although it's layer is 0) 00518 distance(innode)=0 ! seed node distance is 0 00519 allocate(inlayer(virtaggrsize)) 00520 allocate(newlayer(virtaggrsize)) 00521 allocate(roundernodes(maxasize)) 00522 nneigs=0 00523 ninnodes=1 00524 inlayer(1)=innode 00525 ok=can_add_roundlayer(0,ninnodes,inlayer,newlayer, & 00526 roundernodes,neighood,minasize,maxasize,stat,distance, & 00527 nneigs,nodes,rowstart,colnrs) 00528 if (nneigs<minasize) then 00529 ok=0 00530 endif 00531 deallocate(roundernodes) 00532 deallocate(newlayer) 00533 deallocate(inlayer) 00534 return 00535 end function lets_colour3 00536 00537 recursive function can_add_roundlayer3(layer,ninnodes, & 00538 inlayer,newlayer,roundernodes,neighood,& 00539 minasize,maxasize,stat,distance, & 00540 nneigs,nodes,rowstart,colnrs) result(ok) 00541 use globals 00542 Implicit None 00543 integer :: ok 00544 integer,intent(in) :: layer ! the layer # of innode 00545 integer,intent(inout) :: ninnodes ! number of nodes on the last layer 00546 integer,dimension(:),pointer :: inlayer ! nodes on the last outer layer 00547 integer,dimension(:),pointer :: newlayer ! storage for new layer 00548 integer,dimension(:),pointer :: roundernodes ! storage for rounding 00549 integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00550 integer,intent(in) :: minasize,maxasize ! aggregate limits 00551 integer,dimension(:),pointer :: stat 00552 integer,dimension(:),pointer :: rowstart,colnrs,distance 00553 integer,intent(inout) :: nneigs 00554 integer,dimension(:) :: nodes 00555 integer :: nxtlayer ! nxtlayer -- the next layer to look at 00556 integer :: i,j,k,d,rs,re,cn,layernode,ntoadd,nrounders,mind,maxd 00557 logical :: layerfits 00558 !logical :: checklayerfit=.false. 00559 logical :: checklayerfit=.true. 00560 !logical :: dorounding=.false. 00561 logical :: dorounding=.true. 00562 integer,save :: allmind=D_MAXINT 00563 00564 ! we assume, symmetric structure 00565 nxtlayer=layer+1 00566 layerfits=.true. 00567 ntoadd=0 00568 if (layerfits.and.layer<=2*neighood) then 00569 aaa:do k=1,ninnodes 00570 layernode=inlayer(k) 00571 rs=rowstart(layernode) 00572 re=rowstart(layernode+1)-1 00573 do i=rs,re ! put the neighs into the list if needed 00574 cn=colnrs(i) 00575 if (stat(cn)<=0) then !node not marked yet 00576 ! wheather this node has been seen from inlayer already?: 00577 if (stat(cn)==0) then ! first time 00578 !if (layer<=neighood.and.nneigs+ntoadd+1>maxasize) then 00579 ! layerfits=.false. 00580 ! exit aaa 00581 !endif 00582 ntoadd=ntoadd+1 00583 newlayer(ntoadd)=cn 00584 distance(cn)=distance(layernode)+1 00585 endif 00586 stat(cn)=stat(cn)-1 ! We use stat here for finding weights for 00587 ! rounding step... 00588 endif 00589 enddo 00590 enddo aaa 00591 endif 00592 if (layerfits.and.layer<=neighood) then !{ 00593 ! do the rounding step: 00594 nrounders=0 00595 ! find the min/max weight for the newlayer: 00596 mind=D_MAXINT 00597 maxd=0 00598 do k=1,ntoadd 00599 d=-stat(newlayer(k)) 00600 mind=min(d,mind) 00601 maxd=max(d,maxd) 00602 enddo 00603 if (mind<allmind) then 00604 allmind=mind 00605 endif 00606 !if (debu) then 00607 ! print *,'mind,maxd:',mind,maxd 00608 !endif 00609 if (mind<maxd.or.mind>allmind) then 00610 bbb:do k=1,ntoadd 00611 if (stat(newlayer(k))==-maxd) then ! this is a rounder node 00612 nrounders=nrounders+1 00613 if (checklayerfit.and.nneigs+ninnodes+1>maxasize) then 00614 layerfits=.false. 00615 exit bbb 00616 endif 00617 ninnodes=ninnodes+1 00618 inlayer(ninnodes)=newlayer(k) ! add this rounder node 00619 roundernodes(nrounders)=newlayer(k) 00620 newlayer(k)=-1 ! this is taken out from the new layer 00621 endif 00622 enddo bbb 00623 endif 00624 ! if rounding of the incoming layer succeeded, we can add this to the 00625 ! aggregate... 00626 ! ... done later 00627 if (.not.layerfits) then 00628 nrounders=0 00629 endif 00630 ! now still add neighbours of rounded nodes to newlayer: 00631 ! first take out holes from the newlayer 00632 !if (nrounders>0) then 00633 j=0 00634 do i=1,ntoadd 00635 if (newlayer(i)/=-1) then 00636 j=j+1 00637 if (j<i) then 00638 newlayer(j)=newlayer(i) 00639 endif 00640 endif 00641 enddo 00642 ntoadd=j 00643 !endif 00644 ! add rounder's neighbours to the newlayer: 00645 ccc:do k=1,nrounders 00646 layernode=roundernodes(k) 00647 rs=rowstart(layernode) 00648 re=rowstart(layernode+1)-1 00649 do i=rs,re ! put the neighs into the list if needed 00650 cn=colnrs(i) 00651 if (stat(cn)==0) then !node not seen yet 00652 !if (nneigs+ntoadd+1>maxasize) then 00653 ! layerfits=.false. 00654 ! exit ccc 00655 !endif 00656 ntoadd=ntoadd+1 00657 newlayer(ntoadd)=cn 00658 distance(cn)=distance(layernode)+1 00659 endif 00660 enddo 00661 enddo ccc 00662 !elseif (ntoadd>0) then ! the layer is larger than neighood, and we have 00663 ! ! found some nodes around... 00664 ! ok=1 ! was completely wrong here... 00665 endif !(layer<=neighood) } 00666 if (checklayerfit.and.layer<=neighood.and.nneigs+ninnodes>maxasize) then ! this could 00667 ! happen, if no rounding nodes got added... 00668 layerfits=.false. 00669 endif 00670 if (layerfits) then ! add the layer 00671 do k=1,ninnodes 00672 stat(inlayer(k))=layer+1 00673 nneigs=nneigs+1 00674 nodes(nneigs)=inlayer(k) 00675 enddo 00676 if (layer<=2*neighood) then ! recurse to the next layer 00677 do k=1,ntoadd 00678 stat(newlayer(k))=D_PENDING-1 ! just temporary pos. marker here... 00679 enddo 00680 ninnodes=ntoadd 00681 inlayer(1:ninnodes)=newlayer(1:ntoadd) ! (this does not incl. rounders) 00682 ok=can_add_roundlayer(nxtlayer,ninnodes, & 00683 inlayer,newlayer,roundernodes,neighood, & 00684 minasize,maxasize,stat,distance, & 00685 nneigs,nodes,rowstart,colnrs) 00686 else 00687 ok=2 00688 endif 00689 return 00690 else ! clean up the possible aggregate nodes and also on pending layers: 00691 !print *,'layer=',layer 00692 !print *,'minasize,maxasize:',minasize,maxasize 00693 !print *,'nneigs,ninnodes,ntoadd:',nneigs,ninnodes,ntoadd 00694 !print *,'nodes:',nodes(1:nneigs) 00695 !print *,'inlayer:',inlayer(1:ninnodes) 00696 !print *,'newlayer:',newlayer(1:ntoadd) 00697 !print *,'--------------------------------------' 00698 do k=1,ninnodes 00699 distance(inlayer(k))=0 00700 enddo 00701 do k=1,ntoadd 00702 distance(newlayer(k))=0 00703 enddo 00704 do k=1,nneigs 00705 distance(nodes(k))=0 00706 enddo 00707 !if (nneigs<minasize) then 00708 ! do k=1,ninnodes 00709 ! stat(inlayer(k))=D_PENDING 00710 ! enddo 00711 ! do k=1,ntoadd 00712 ! stat(newlayer(k))=D_PENDING 00713 ! enddo 00714 ! do k=1,nneigs 00715 ! stat(nodes(k))=D_PENDING 00716 ! enddo 00717 ! ok=0 00718 !else 00719 ! do k=1,ninnodes 00720 ! stat(inlayer(k))=0 00721 ! enddo 00722 ! do k=1,ntoadd 00723 ! stat(newlayer(k))=0 00724 ! enddo 00725 ! ok=1 00726 !endif 00727 if (nneigs+ninnodes<minasize) then 00728 !!if (nneigs<minasize) then 00729 ! add the found nodes to the inlayer to be able to give them appropriate 00730 ! failing structure number outside: 00731 do k=1,ninnodes 00732 nneigs=nneigs+1 00733 nodes(nneigs)=inlayer(k) 00734 enddo 00735 !do k=1,ntoadd 00736 ! nneigs=nneigs+1 00737 ! nodes(nneigs)=newlayer(k) 00738 !enddo 00739 ok=0 00740 endif 00741 return 00742 endif 00743 end function can_add_roundlayer3 00744 00745 subroutine lets_colour2(innode,neighood,minasize,maxasize,nneigs,nodes, & 00746 stat,rowstart,colnrs,aggrnum) 00747 Implicit None 00748 integer, intent(in) :: innode ! the node who's neighood is being built 00749 integer, intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00750 integer, intent(in) :: minasize,maxasize ! aggregate limits 00751 integer, intent(in out) :: nneigs 00752 integer, dimension(:), pointer :: stat 00753 integer, dimension(:), pointer :: rowstart,colnrs 00754 integer, dimension(:) :: nodes 00755 integer, dimension(:), pointer,optional :: aggrnum 00756 nneigs=0 00757 stat(innode)=1 ! mark it as if on layer 1 (although it's layer is 0) 00758 call colouring_neighs2(0,innode,neighood,minasize,maxasize,stat, & ! ok 00759 nneigs,nodes,rowstart,colnrs,aggrnum) 00760 end subroutine lets_colour2 00761 00762 recursive subroutine colouring_neighs2(layer,innode,neighood,minasize,maxasize,stat, & 00763 nneigs,nodes,rowstart,colnrs,aggrnum) 00764 use globals 00765 Implicit None 00766 integer,intent(in) :: layer ! the layer # of innode 00767 integer, intent(in) :: innode ! the node who's neighood is being built 00768 integer, intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00769 integer, intent(in) :: minasize,maxasize ! aggregate limits 00770 integer, dimension(:), pointer :: stat 00771 integer, dimension(:), pointer :: rowstart,colnrs 00772 integer, intent(in out) :: nneigs 00773 integer, dimension(:) :: nodes 00774 integer, dimension(:), pointer,optional :: aggrnum 00775 integer :: nneigs_in,nxtlayer ! nxtlayer -- the next layer to look at 00776 integer :: my_nneigs,radius2 00777 integer :: i,j,rs,re,cn 00778 00779 ! we assume, symmetric structure 00780 radius2=neighood+1 00781 nneigs_in=nneigs 00782 nxtlayer=layer+1 00783 rs=rowstart(innode) 00784 re=rowstart(innode+1)-1 00785 do i=rs,re ! put the neighs into the list if needed 00786 cn=colnrs(i) 00787 !print *,'cn=',cn 00788 if (stat(cn)==-1) then !node cn was in shade and not included yet 00789 stat(cn)=nxtlayer 00790 nneigs=nneigs+1 00791 nodes(nneigs)=-cn ! "-" marking previous shade 00792 else if (stat(cn)==0) then ! node cn to be included 00793 stat(cn)=nxtlayer 00794 nneigs=nneigs+1 00795 nodes(nneigs)=cn 00796 else if (stat(cn)==-3) then ! todo: comment it out as this should never 00797 ! happen! 00798 write(stream,*) 'warning: somehow jumped onto a isolated but too small aggr!' 00799 write(stream,*) 'innode,i,layer,nneigs_in,nneigs,cn,stat(cn):', & 00800 innode,i,layer,nneigs_in,nneigs,cn,stat(cn) 00801 write(stream,*) 'stat=' 00802 do j=1,size(stat) 00803 write(stream,'(i2"("i2") " )',ADVANCE='NO') stat(j),j 00804 enddo 00805 write(stream,*)' ' 00806 write(stream,*) 'aggrnum=' 00807 do j=1,size(aggrnum) 00808 write(stream,'(i2"("i2") " )',ADVANCE='NO') aggrnum(j),j 00809 enddo 00810 write(stream,*)' ' 00811 write(stream,*) ' min/max aggr size is:',minasize,maxasize 00812 stop 00813 endif 00814 enddo 00815 my_nneigs=nneigs 00816 !if (nxtlayer<2*neighood) then 00817 if (nxtlayer<radius2) then ! if there are more layers to process 00818 do i=nneigs_in+1,my_nneigs ! recursively check the neighbours 00819 if (i<=maxasize) then 00820 call colouring_neighs2(nxtlayer,abs(nodes(i)),neighood, & 00821 minasize,maxasize,stat,nneigs,nodes,rowstart,colnrs,aggrnum) 00822 endif 00823 enddo 00824 endif 00825 return 00826 end subroutine colouring_neighs2 00827 00828 recursive function aggregate_to_neighbour(innode,dist,aggrnum, & 00829 rowstart,colnrs) result(ok) 00830 Implicit None 00831 logical :: ok 00832 integer, intent(in) :: innode ! the node who's neighs's aggrnum is looked 00833 integer, intent(in out) :: dist ! distance to dig 00834 integer, dimension(:), intent(in out) :: aggrnum 00835 integer, dimension(:), pointer :: rowstart,colnrs 00836 integer :: i,j,rs,re,cn 00837 00838 ! we assume, symmetric structure 00839 rs=rowstart(innode) 00840 re=rowstart(innode+1)-1 00841 do i=rs,re ! look through the neigbours 00842 cn=colnrs(i) 00843 if (aggrnum(cn)/=0) then !we are done! 00844 aggrnum(innode)=aggrnum(cn) 00845 ok=.true. 00846 return 00847 else if (dist>0) then ! dig further around... 00848 dist=dist-1 00849 if (.not.aggregate_to_neighbour(innode,dist,aggrnum, & 00850 rowstart,colnrs)) then 00851 ok=.false. 00852 return 00853 endif 00854 endif 00855 enddo 00856 ! we have look everywhere deep enough, but there are no OK aggregates... 00857 ok=.false. 00858 return 00859 end function aggregate_to_neighbour 00860 00861 subroutine color_print_aggrs(n,aggrnum,coarse_aggrnum,overwrite,owner) 00862 use globals, only: stream 00863 ! print aggregates in case of small regular 2D problems: 00864 Implicit None 00865 integer,intent(in) :: n 00866 integer,dimension(:),pointer :: aggrnum 00867 integer,dimension(:),pointer,optional :: coarse_aggrnum 00868 logical,optional :: overwrite 00869 integer,dimension(:),pointer,optional :: owner 00870 integer :: i,j,k,kk 00871 integer, parameter :: isolcol1=6 00872 integer, parameter :: isolcol2=7 00873 k=sqrt(1.0*n) 00874 if (k*k==n.and.k<270) then 00875 if (present(overwrite).and.overwrite) then 00876 call cursor_up(k) 00877 endif 00878 if (k<=136) then 00879 kk=0 00880 do i=1,k 00881 do j=1,k 00882 kk=kk+1 00883 if (present(coarse_aggrnum)) then 00884 if (aggrnum(kk)>0) then 00885 if (coarse_aggrnum(aggrnum(kk))>0) then 00886 call cprint(char(modulo(aggrnum(kk)/10,10)+48), & 00887 coarse_aggrnum(aggrnum(kk))) 00888 call cprint(char(modulo(aggrnum(kk),10)+48), & 00889 coarse_aggrnum(aggrnum(kk))) 00890 else 00891 call cprint('#',isolcol2) 00892 call cprint('#',isolcol2) 00893 endif 00894 else 00895 call cprint('#',isolcol1) 00896 call cprint('#',isolcol1) 00897 endif 00898 elseif (present(owner)) then 00899 if (aggrnum(kk)>0) then 00900 call cprint(char(modulo(aggrnum(kk)/10,10)+48),owner(kk)) 00901 call cprint(char(modulo(aggrnum(kk),10)+48),owner(kk)) 00902 else 00903 call cprint('#',isolcol1) 00904 call cprint('#',isolcol1) 00905 endif 00906 else 00907 if (aggrnum(kk)>0) then 00908 call cprint(char(modulo(aggrnum(kk)/10,10)+48),aggrnum(kk)) 00909 call cprint(char(modulo(aggrnum(kk),10)+48),aggrnum(kk)) 00910 else 00911 call cprint('#',isolcol1) 00912 call cprint('#',isolcol1) 00913 endif 00914 endif 00915 enddo 00916 write(stream,*) 00917 enddo 00918 else 00919 kk=0 00920 do i=1,k 00921 do j=1,k 00922 kk=kk+1 00923 if (present(coarse_aggrnum)) then 00924 if (aggrnum(kk)>0) then 00925 if (coarse_aggrnum(aggrnum(kk))>0) then 00926 call cprint(char(modulo(aggrnum(kk),10)+48), & 00927 coarse_aggrnum(aggrnum(kk))) 00928 else 00929 call cprint('#',isolcol2) 00930 endif 00931 else 00932 call cprint('#',isolcol1) 00933 endif 00934 elseif (present(owner)) then 00935 if (aggrnum(kk)>0) then 00936 call cprint(char(modulo(aggrnum(kk),10)+48),owner(kk)) 00937 else 00938 call cprint('#',isolcol1) 00939 endif 00940 else 00941 if (aggrnum(kk)>0) then 00942 call cprint(char(modulo(aggrnum(kk),10)+48),aggrnum(kk)) 00943 else 00944 call cprint('#',isolcol1) 00945 endif 00946 endif 00947 enddo 00948 write(stream,*) 00949 enddo 00950 endif 00951 else 00952 write(stream,*) 'cannot colorprint aggrnums, k,n=',k,n 00953 endif 00954 end subroutine color_print_aggrs 00955 00956 subroutine cursor0() 00957 use globals, only: stream 00958 implicit none 00959 write(stream,'(2a,i1,a,i1,a1)',advance='no') & 00960 char(27),'[',0,';',0,'f' 00961 end subroutine cursor0 00962 00963 subroutine cursor_up(n) 00964 use globals, only: stream 00965 implicit none 00966 integer,intent(in) :: n 00967 if (n<10) then 00968 write(stream,'(2a,i1,a1)',advance='no') & 00969 char(27),'[',n,'A' 00970 elseif (n<100) then 00971 write(stream,'(2a,i2,a1)',advance='no') & 00972 char(27),'[',n,'A' 00973 else 00974 write(stream,'(2a,i3,a1)',advance='no') & 00975 char(27),'[',n,'A' 00976 endif 00977 end subroutine cursor_up 00978 00979 subroutine cursor_down(n) 00980 use globals, only: stream 00981 implicit none 00982 integer,intent(in) :: n 00983 if (n<10) then 00984 write(stream,'(2a,i1,a1)',advance='no') & 00985 char(27),'[',n,'B' 00986 elseif (n<100) then 00987 write(stream,'(2a,i2,a1)',advance='no') & 00988 char(27),'[',n,'B' 00989 else 00990 write(stream,'(2a,i3,a1)',advance='no') & 00991 char(27),'[',n,'B' 00992 endif 00993 end subroutine cursor_down 00994 00995 subroutine cprint(c,col) 00996 use globals, only: stream 00997 implicit none 00998 character,intent(in) :: c 00999 integer,intent(in) :: col 01000 integer :: b=0,i,j 01001 integer,dimension(20) :: t=(/ (i,i=41,46),(i,i=90,96),(i,i=100,106) /) 01002 j=modulo(col,size(t))+1 01003 if (t(j)<10) then 01004 write(stream,'(2a,i1,a,i1,a1,a1,a1,a1,i1,a1)',advance='no') & 01005 char(27),'[',b,';',t(j),'m',c,char(27),'[',0,'m' 01006 elseif (t(j)<100) then 01007 write(stream,'(2a,i1,a,i2,a1,a1,a1,a1,i1,a1)',advance='no') & 01008 char(27),'[',b,';',t(j),'m',c,char(27),'[',0,'m' 01009 else 01010 write(stream,'(2a,i1,a,i3,a1,a1,a1,a1,i1,a1)',advance='no') & 01011 char(27),'[',b,';',t(j),'m',c,char(27),'[',0,'m' 01012 endif 01013 end subroutine cprint 01014 01015 subroutine cprintall(c,col) ! to test what colours are there... 01016 !do i=1,512 01017 !call cprintall('#',i) 01018 !print *,i 01019 !enddo 01020 !stop 01021 use globals, only: stream 01022 implicit none 01023 character,intent(in) :: c 01024 integer,intent(in) :: col 01025 integer :: b=0,i,j 01026 if (col<10) then 01027 write(stream,'(2a,i1,a,i1,a1,a1,a1,a1,i1,a1)',advance='no') & 01028 char(27),'[',b,';',col,'m',c,char(27),'[',0,'m' 01029 elseif (col<100) then 01030 write(stream,'(2a,i1,a,i2,a1,a1,a1,a1,i1,a1)',advance='no') & 01031 char(27),'[',b,';',col,'m',c,char(27),'[',0,'m' 01032 else 01033 write(stream,'(2a,i1,a,i3,a1,a1,a1,a1,i1,a1)',advance='no') & 01034 char(27),'[',b,';',col,'m',c,char(27),'[',0,'m' 01035 endif 01036 end subroutine cprintall 01037 01038 !------------------------------------------------------ 01039 end Module Aggregate_mod 01040 !------------------------------------------------------
1.7.3-20110217