DOUG 0.2

Aggregate_mod.F90

Go to the documentation of this file.
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 !------------------------------------------------------