DOUG 0.2

SpMtx_aggregation.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 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