DOUG 0.2

aggr.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 
00067 program main_aggr
00068 
00069   use Distribution_mod
00070   use Partitioning_mod
00071   use Partitioning_aggr_mod
00072   use Partitioning_full_mod
00073   use Partitioning_metis_mod
00074   use Preconditioner_mod
00075   use FinePreconditioner_complete_mod
00076   use FinePreconditioner_sgs_mod
00077   use CoarsePreconditioner_smooth_mod
00078   use CoarsePreconditioner_robust_mod
00079   use Mesh_class
00080   use Mesh_plot_mod
00081   use SpMtx_mods
00082   use Vect_mod
00083   use DenseMtx_mod
00084   use solvers_mod
00085   use Aggregate_mod
00086   use Aggregate_utils_mod
00087   use CoarseMtx_mod
00088   use CoarseAllgathers
00089 
00090   implicit none
00091 
00092 #include<doug_config.h>
00093 
00094 #ifdef D_COMPLEX
00095 #define float complex
00096 #else
00097 #define float real
00098 #endif
00099 
00100   float(kind=rk), dimension(:), pointer :: xl !< local solution vector
00101   float(kind=rk), dimension(:), pointer :: x  !< global solution on master
00102   float(kind=rk), dimension(:), pointer :: sol, rhs  !< for testing solver
00103 
00104   ! Partitioning
00105   integer               :: nparts !< number of partitons to partition a mesh
00106   integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/)
00107 
00108   type(ConvInf) :: resStat
00109 
00110   real(kind=rk) :: t,t1
00111   integer :: i,j,k,it,ol
00112   float(kind=rk), dimension(:), pointer :: r, y
00113   character :: str
00114   character(len=40) :: frm
00115   float(kind=rk) :: cond_num,nrm
00116   integer,allocatable :: nodes(:), inds(:)
00117 
00118   type(Distribution) :: D !< mesh and matrix distribution
00119   type(Partitionings) :: P !< fine and coarse aggregates
00120   type(FinePreconditioner) :: FP
00121   type(CoarsePreconditioner) :: CP
00122 
00123   ! Init DOUG
00124   call DOUG_Init()
00125 
00126   ! Master participates in calculations as well
00127   nparts = numprocs
00128 
00129   D = Distribution_NewInit(sctls%input_type,nparts,part_opts)
00130 
00131   if (sctls%levels>1.or.(numprocs==1.and.sctls%levels==1)) then !todo remove
00132     P = Partitionings_New()
00133     call Partitionings_aggr_InitFine(P,D)
00134     ! profile info
00135     if(pstream/=0) then
00136       write(pstream, "(I0,':fine aggregates:',I0)") myrank, P%fAggr%inner%nagr
00137     end if
00138 
00139     !if (sctls%plotting>=2) then
00140     !   call SpMtx_writeLogicalValues(A, D%A%strong, 'strong.txt')
00141     !end if
00142 
00143     call Mesh_printInfo(D%mesh)
00144     
00145     if (numprocs==1.and.sctls%plotting==2.and.D%mesh%nell>0) then
00146       call Mesh_pl2D_plotAggregate(P%fAggr%inner,D%mesh,&
00147                       D%A%strong_rowstart,D%A%strong_colnrs,&
00148                       mctls%assembled_mtx_file, &
00149                                  INIT_CONT_END=D_PLPLOT_INIT)
00150                                  !D_PLPLOT_END)
00151     endif
00152     
00153     CP = CoarsePreconditioner_New()
00154     ! Testing coarse matrix and aggregation through it:
00155     if (sctls%coarse_method<=1) then ! if not specified or ==1
00156       call CoarsePreconditioner_smooth_Init(CP, D, P)
00157 
00158     else if (sctls%coarse_method==2) then
00159       ! use the Robust Coarse Spaces algorithm
00160       call CoarsePreconditioner_robust_Init(CP, D, P)
00161 
00162     else
00163       write(stream,'(A," ",I2)') 'Wrong coarse method', sctls%coarse_method
00164       call DOUG_abort('Error in aggr', -1)
00165     endif
00166               
00167     ! coarse aggregates
00168     if (numprocs==1) then
00169       call Partitionings_aggr_InitCoarse(P,D,CP%AC)
00170       !call Aggrs_readFile_coarse(P%cAggr, "aggregates.txt")
00171 
00172       ! profile info
00173       if(pstream/=0) then
00174         write(pstream, "(I0,':coarse aggregates:',I0)") myrank, P%cAggr%inner%nagr
00175       end if
00176 
00177       if (sctls%plotting==2) then
00178          call Aggr_writeFile(P%fAggr%inner, 'aggr2.txt', P%cAggr%inner)
00179       end if
00180       if (sctls%plotting==2.and.D%mesh%nell>0) then
00181         !print *,'press Key<Enter>'
00182         !read *,str
00183         call Mesh_pl2D_plotAggregate(P%fAggr%inner,D%mesh,&
00184                         D%A%strong_rowstart,D%A%strong_colnrs,&
00185                         mctls%assembled_mtx_file, &
00186                         caggrnum=P%cAggr%inner%num, &
00187                       INIT_CONT_END=D_PLPLOT_END)!, &
00188                       !INIT_CONT_END=D_PLPLOT_CONT)!, &
00189                                   ! D_PLPLOT_END)
00190       endif
00191       write(stream,*)'# coarse aggregates:',P%cAggr%inner%nagr
00192 
00193     endif 
00194 
00195   else ! 1 level, several procs
00196     ! required for metis coarse subdomains
00197     P = Partitionings_New()
00198     call Partitionings_aggr_InitFine(P,D)
00199   endif
00200 
00201   if (numprocs>1) then
00202     !call Partitionings_full_InitCoarse(P,D)
00203     if (sctls%num_subdomains<=0) sctls%num_subdomains=1
00204     call Partitionings_metis_InitCoarse(P,D,sctls%num_subdomains)
00205   end if
00206 
00207   ! overlap for subdomains
00208   if (sctls%overlap<0) then ! autom. overlap from smoothing
00209     ol = max(sctls%smoothers,0)
00210   else
00211     ol = sctls%overlap
00212   endif
00213 
00214   FP = FinePreconditioner_New(D)
00215   call FinePreconditioner_Init(FP, D, P, ol)
00216   if (sctls%fine_method==FINE_PRECONDITIONER_TYPE_NONE) then
00217     ! do nothing
00218   else if (sctls%fine_method==FINE_PRECONDITIONER_TYPE_COMPLETE) then
00219     call FinePreconditioner_complete_Init(FP)
00220   else if (sctls%fine_method==FINE_PRECONDITIONER_TYPE_SGS) then
00221     if (sctls%num_iters<=0) sctls%num_iters=3
00222     call FinePreconditioner_sgs_Init(FP,sctls%num_iters)
00223   else
00224     write(stream,'(A," ",I2)') 'Wrong fine method', sctls%fine_method
00225     call DOUG_abort('Error in aggr', -1)
00226   end if
00227 
00228   if (numprocs==1) then
00229     call AggrInfo_Destroy(P%cAggr)
00230     call AggrInfo_Destroy(P%fAggr)
00231   else
00232     ! call Aggrs_writeFile(M, P%fAggr, CP%cdat, "aggregates.txt")
00233     if (sctls%levels>1) call AggrInfo_Destroy(P%fAggr)
00234   end if
00235 
00236   ! Testing UMFPACK:
00237   allocate(sol(D%A%nrows))
00238   allocate(rhs(D%A%ncols))
00239 
00240   ! Solve the system
00241   allocate(xl(D%mesh%nlf))
00242   xl = 0.0_rk
00243 
00244   select case(sctls%solver)
00245   case (DCTL_SOLVE_PCG)
00246      ! Preconditioned conjugate gradient
00247      t1 = MPI_WTIME()
00248      write(stream,*)'calling pcg_weigs'
00249      call pcg_weigs(D, x=xl,&
00250           finePrec=FP,coarsePrec=CP,&
00251           it=it,cond_num=cond_num)
00252      t=MPI_WTIME()-t1
00253      write(stream,*) 'time spent in pcg():',t
00254      t1=total_setup_time()
00255      write(stream,*) '    ...of which setup:',t1
00256      write(stream,*) '       ...of which factorisation:', &
00257         total_factorisation_time()
00258      write(stream,*) 'solve time without setup:',t-t1
00259   case default
00260      call DOUG_abort('[DOUG main] : Wrong solution method specified', -1)
00261   end select
00262 
00263   if (numprocs>1) then
00264     ! Assemble result on master
00265     if (ismaster()) then
00266        print *, "freedoms", D%mesh%ngf
00267       allocate(x(D%mesh%ngf)); x = 0.0_rk
00268     end if
00269     call Vect_Gather(xl, x, D%mesh)
00270     if (ismaster().and.sctls%verbose>2.and.(size(x) <= 100)) &
00271       call Vect_Print(x, 'sol > ')
00272     if (ismaster()) then
00273        call WriteSolutionToFile(x)
00274     end if
00275 
00276   else
00277     if (sctls%verbose>2.and.size(xl)<=100) then
00278       call Vect_Print(xl, 'sol > ')
00279     endif
00280     call WriteSolutionToFile(xl)
00281   endif
00282 
00283 
00284   ! Destroy objects
00285   call Mesh_Destroy(D%mesh)
00286   call SpMtx_Destroy(D%A)
00287   call ConvInf_Destroy(resStat)
00288   if (associated(xl)) deallocate(xl)
00289   if (associated(x)) deallocate(x)
00290   if (associated(sol)) deallocate(sol)
00291 
00292   call DOUG_Finalize()
00293 
00294 
00295 end program main_aggr