DOUG 0.2

CoarsePreconditioner_smooth.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 CoarsePreconditioner_smooth_mod
00024   use Preconditioner_base_mod
00025   use SpMtx_distribution_mod
00026   use SpMtx_operation
00027   use SpMtx_aggregation
00028   use CoarseMtx_mod
00029 
00030   implicit none
00031 
00033   type CoarseSpace
00034      integer :: nsupports !< number of nodes in the coarse space
00035      integer,pointer :: support_nodes(:) !< fine node indices for the coarse supports
00036      integer,pointer :: support_bounds(:) !< bounds for support_inds
00037      integer,pointer :: esupport_nodes(:) !< expanded to the overlap: neighbour coarse nodes
00038      integer,pointer :: esupport_bounds(:) !< bounds for esupport_nodes
00039   end type CoarseSpace
00040   
00041   private
00042   public :: CoarsePreconditioner_smooth_Init
00043 
00044 contains
00045 
00047   subroutine CoarsePreconditioner_smooth_Init(CP, D, P)
00048     type(CoarsePreconditioner),intent(inout) :: CP
00049     type(Distribution),intent(inout) :: D
00050     type(Partitionings),intent(in) :: P
00051 
00052     type(CoarseSpace) :: CS
00053     integer,pointer :: aggrnum(:)
00054     integer :: nagr, i
00055 
00056     CP%type = COARSE_PRECONDITIONER_TYPE_SMOOTH
00057 
00058     if (numprocs>1) then
00059       allocate(aggrnum(D%mesh%nlf))
00060       nagr = P%fAggr%inner%nagr
00061       aggrnum = 0
00062       aggrnum(1:D%mesh%ninner) = P%fAggr%inner%num
00063       call setup_aggr_cdat(CP%cdat, CP%cdat_vec, nagr, D%mesh%ninner,aggrnum,D%mesh)
00064 
00065       call SpMtx_find_strong(A=D%A,alpha=P%strong_conn1,A_ghost=D%A_ghost)
00066       call SpMtx_exchange_strong(D%A,D%A_ghost,D%mesh)
00067       call SpMtx_symm_strong(D%A,D%A_ghost,.false.)
00068       call SpMtx_unscale(D%A)
00069     end if
00070 
00071     call IntRestBuild(D%A,P%fAggr%inner,CP%R,D%A_ghost)
00072 
00073     if (numprocs>1) then    
00074       CS = CoarseSpace_Init(CP%R)
00075       call CoarseData_Copy(CP%cdat,CP%cdat_vec)
00076       call CoarseSpace_Expand(CS,CP%R,D%mesh,CP%cdat)
00077     end if
00078     
00079     if (numprocs>1) then
00080       call CoarseMtxBuild(D%A,CP%cdat%LAC,CP%R,D%mesh%ninner,D%A_ghost)
00081     else
00082       call CoarseMtxBuild(D%A,CP%AC,CP%R,D%mesh%ninner)      
00083     end if
00084 
00085     if (numprocs>1) then
00086       call KeepGivenRowIndeces(CP%R, (/(i,i=1,P%fAggr%inner%nagr)/),.false.)
00087     end if
00088 
00089   end subroutine CoarsePreconditioner_smooth_Init
00090 
00092   function CoarseSpace_Init(Restrict) result(CS)
00093     type(SpMtx), intent(in) :: Restrict
00094     type(CoarseSpace) :: CS
00095 
00096     integer, allocatable :: nnodes(:), cnodes(:)
00097     integer :: i, isupport
00098 
00099     CS%nsupports = Restrict%nrows
00100 
00101     ! count the number of fine mesh nodes in each coarse node support
00102     allocate(nnodes(CS%nsupports))
00103     nnodes = 0
00104     do i = 1,Restrict%nnz
00105       isupport = Restrict%indi(i)
00106       if (isupport<=CS%nsupports) then
00107         nnodes(isupport) = nnodes(isupport)+1
00108       end if
00109     end do
00110 
00111     ! scan-reduce to get bounds
00112     allocate(CS%support_bounds(CS%nsupports+1))
00113     CS%support_bounds(1) = 1
00114     do i = 1,CS%nsupports
00115       CS%support_bounds(i+1) = CS%support_bounds(i)+nnodes(i)
00116     end do
00117     !write(stream,*) "CS%support_bounds", CS%support_bounds
00118 
00119     ! store indices
00120     allocate(CS%support_nodes(CS%support_bounds(CS%nsupports+1)-1))
00121     allocate(cnodes(CS%nsupports))
00122     cnodes = CS%support_bounds(1:CS%nsupports)
00123     do i = 1,Restrict%nnz
00124       isupport = Restrict%indi(i)
00125       if (isupport<=CS%nsupports) then
00126         CS%support_nodes(cnodes(isupport)) = Restrict%indj(i)
00127         cnodes(isupport) = cnodes(isupport)+1
00128       end if
00129     end do
00130 
00131     !write(stream,*) "CS%support_nodes", CS%support_nodes
00132   end function CoarseSpace_Init
00133 
00137   subroutine CoarseSpace_Expand(CS,R,M,cdat)
00138     use CoarseAllgathers, only: CoarseData
00139     type(CoarseSpace), intent(inout) :: CS
00140     type(SpMtx), intent(inout) :: R !< Restriction matrix
00141     type(Mesh), intent(in) :: M
00142     type(CoarseData), intent(inout) :: cdat
00143 
00144     real(kind=rk), pointer :: val(:)
00145     type(SpMtx) :: iR, i2ecR, eR
00146     integer :: i
00147     
00148     ! ---- collect restrict values
00149     call collectRestrictValues(R, i2ecR)
00150 
00151     ! extend cdat with the neighbour coarse nodes
00152     i2ecR%indj = M%gl_fmap(i2ecR%indj)
00153     call add_indices(cdat, i2ecR%indi)
00154 
00155     ! localize i2ecR and add new elements to the restriction matrix
00156     i2ecR%indi = cdat%gl_cfmap(i2ecR%indi)
00157     i2ecR%nrows = cdat%nlfc
00158     i2ecR%ncols = M%nlf
00159     iR = SpMtx_add(R, i2ecR, 1.0_rk, 1.0_rk)
00160     call KeepGivenColumnIndeces(iR,(/(i,i=1,M%ninner)/), keepShape=.TRUE.)
00161 
00162     ! ---- distribute restrict values
00163     call distributeRestrictValues(iR, eR)
00164 
00165     ! extend cdat once more, this may contain third partition coarse nodes
00166     eR%indj = M%gl_fmap(eR%indj)
00167     call add_indices(cdat, eR%indi)
00168 
00169     ! localize eR and add new elements
00170     eR%indi = cdat%gl_cfmap(eR%indi)
00171     eR%nrows = cdat%nlfc
00172     eR%ncols = M%nlf
00173 
00174     ! add internal and external elements
00175     call SpMtx_Destroy(R)
00176     R = SpMtx_add(iR, eR, 1.0_rk, 1.0_rk)
00177 
00178     call SpMtx_Destroy(eR)
00179     call SpMtx_Destroy(iR)
00180 
00181   contains
00182     subroutine add_indices(cdata, g_cinds)
00183       type(CoarseData), intent(inout) :: cdata
00184       integer, intent(in) :: g_cinds(:) !< coarse nodes to add (with duplicates)
00185       
00186       integer :: gci,lci,k,nf ! global coarse index, local coarse index
00187       integer,pointer :: lg_cfmap(:)
00188 
00189       ! first mark new indices with negative indices
00190       nf = cdata%nlfc
00191       do k=1,size(g_cinds)
00192         gci = g_cinds(k)
00193         if (cdata%gl_cfmap(gci)==0) then
00194           nf = nf+1
00195           cdata%gl_cfmap(gci)=-nf
00196         end if
00197       end do
00198 
00199       ! now extend lg_cfmap
00200       allocate(lg_cfmap(nf))
00201       lg_cfmap(1:cdata%nlfc) = cdata%lg_cfmap
00202       nf = cdata%nlfc
00203       do k=1,size(g_cinds)
00204         gci = g_cinds(k)
00205         if (cdata%gl_cfmap(gci)<0) then
00206           nf = nf+1
00207           lci = -cdata%gl_cfmap(gci)
00208           cdata%gl_cfmap(gci) = lci
00209           lg_cfmap(lci) = gci
00210         end if
00211       end do
00212       
00213       cdata%nlfc = nf
00214       deallocate(cdata%lg_cfmap)
00215       cdata%lg_cfmap => lg_cfmap
00216       
00217     end subroutine add_indices
00218 
00219     subroutine collectRestrictValues(R, eR)
00220       type(SpMtx),intent(in) :: R
00221       type(SpMtx),intent(out) :: eR
00222 
00223       type(indlist), allocatable :: smooth_sends(:)
00224       integer :: k, i, j, isupport, ptn, ngh
00225 
00226       ! first count
00227       allocate(smooth_sends(M%nnghbrs))
00228       smooth_sends%ninds=0
00229       do k=1,R%nnz
00230          isupport = R%indi(k)
00231          j = R%indj(k)
00232          ptn = M%eptnmap(M%lg_fmap(j))
00233          ! if local coarse node value to be prolonged to external fine node value
00234          if (isupport<=CS%nsupports.and.ptn/=myrank+1) then
00235            ! find neighbour number
00236            do i=1,M%nnghbrs
00237              if (M%nghbrs(i)+1==ptn) then
00238                ngh = i
00239                exit
00240              end if
00241            end do
00242            smooth_sends(ngh)%ninds = smooth_sends(ngh)%ninds+1
00243          end if
00244       end do
00245 
00246       ! then collect
00247       do i=1,M%nnghbrs
00248         allocate(smooth_sends(i)%inds(smooth_sends(i)%ninds))
00249       end do
00250       smooth_sends%ninds = 0
00251       do k=1,R%nnz
00252          isupport = R%indi(k)
00253          j = R%indj(k)
00254          ptn = M%eptnmap(M%lg_fmap(j))
00255          ! if local coarse node value to be prolonged to external fine node value
00256          if (isupport<=CS%nsupports.and.ptn/=myrank+1) then
00257             do i=1,M%nnghbrs
00258                if (M%nghbrs(i)+1==ptn) then
00259                   ngh = i
00260                   exit
00261                end if
00262             end do
00263             i = smooth_sends(ngh)%ninds
00264             smooth_sends(ngh)%ninds = i+1
00265             smooth_sends(ngh)%inds(i+1) = k
00266          end if
00267       end do
00268 
00269       eR = SpMtx_exchange(R, smooth_sends, M, cdat%lg_cfmap, M%lg_fmap)
00270 
00271     end subroutine collectRestrictValues
00272     
00273     subroutine distributeRestrictValues(iR, eR)
00274       type(SpMtx),intent(in) :: iR
00275       type(SpMtx),intent(out) :: eR
00276       type(indlist), allocatable :: sends(:)
00277 
00278       integer :: i
00279 
00280       allocate(sends(M%nnghbrs))
00281       do i=1,M%nnghbrs
00282         call SpMtx_findColumnElems(iR, M%ol_inner(i)%inds, sends(i)%inds)
00283         sends(i)%ninds = size(sends(i)%inds)
00284       end do
00285       eR = SpMtx_exchange(iR, sends, M, cdat%lg_cfmap, M%lg_fmap)
00286       
00287       ! deallocate
00288       do i=1,M%nnghbrs
00289         deallocate(sends(i)%inds)
00290       end do
00291       deallocate(sends)
00292     end subroutine distributeRestrictValues
00293 
00294   end subroutine CoarseSpace_Expand
00295 
00296 end module CoarsePreconditioner_smooth_mod