source: src/main/aggr.F90 @ b378423

Revision b378423, 14.8 KB checked in by Oleg Batrashev <ogbash@…>, 2 years ago (diff)

Factor out Aggregates from sparse matrix structure.
  • Property mode set to 100644
Line 
1! DOUG - Domain decomposition On Unstructured Grids
2! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
3! Department of Mathematics, University of Bath
4!
5! This library is free software; you can redistribute it and/or
6! modify it under the terms of the GNU Lesser General Public
7! License as published by the Free Software Foundation; either
8! version 2.1 of the License, or (at your option) any later version.
9!
10! This library is distributed in the hope that it will be useful,
11! but WITHOUT ANY WARRANTY; without even the implied warranty of
12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13! Lesser General Public License for more details.
14!
15! You should have received a copy of the GNU Lesser General Public
16! License along with this library; if not, write to the Free Software
17! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
18! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
19! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
20! mailto:info(at)dougdevel.org)
21
22!> Main program for running DOUG with input files in assembled form.
23!> \section aggregation_running Running the aggregation-based DOUG code: (example)
24!>   <tt>mpirun -np 3 doug_aggr -f doug.ctl</tt>
25!>     where \c doug.ctl may contain the following fields
26!>   \subsection example Input-file example:
27!> \code
28!> solver 2
29!> solve_maxiters 300
30!> method 1
31!> coarse_method 1
32!> levels  2
33!> overlap -1
34!> smoothers 0
35!> input_type 2
36!> symmstruct T
37!> symmnumeric T
38!> # ###################
39!> # aggregate level 1:
40!> radius1 5
41!> strong1 0.67e0
42!> minasize1 2
43!> #maxasize1 19
44!> # aggregate level 2:
45!> radius2 35
46!> strong2 0.67e0
47!> minasize2 2
48!> #maxasize2 96
49!> # ###################
50!> matrix_type 1
51!> number_of_blocks 1
52!> initial_guess 2
53!> start_vec_file ./NOT.DEFINED.start_vec_file
54!> start_vec_type 2
55!> solve_tolerance 1.0e-12
56!> solution_format 2
57!> solution_file ./solution.file
58!> #debug -5
59!> debug 0
60!> verbose 0
61!> plotting 1
62!> assembled_mtx_file Hetero32.txt
63!> \endcode
64program main_aggr
65
66  use doug
67  use main_drivers
68  use Mesh_class
69  use SpMtx_mods
70  use Vect_mod
71  use DenseMtx_mod
72  use solvers_mod
73  use Aggregate_mod
74  use Aggregate_utils_mod
75  use CoarseMtx_mod
76  use CoarseAllgathers
77  use RobustCoarseMtx_mod
78  use pcgRobust_mod
79
80  use aggr_util_mod
81
82  implicit none
83
84#include<doug_config.h>
85
86#ifdef D_COMPLEX
87#define float complex
88#else
89#define float real
90#endif
91
92  type(Mesh)     :: M  !< Mesh
93
94  type(SpMtx)    :: A,A_ghost  !< System matrix (parallel sparse matrix)
95  type(SpMtx)    :: AC  !< coarse matrix
96  type(SpMtx)    :: LA  !< matrix without outer nodes
97  type(SpMtx)    :: Restrict !< Restriction matrix (for operation)
98  type(SumOfInversedSubMtx) :: B_RCS !< B matrix for the Robust Coarse Spaces
99  !type(SpMtx)    :: Rest_cmb !< Restriction matrix (for coarse matrix build)
100
101  float(kind=rk), dimension(:), pointer :: b  !< local RHS
102  float(kind=rk), dimension(:), pointer :: xl !< local solution vector
103  float(kind=rk), dimension(:), pointer :: x  !< global solution on master
104  float(kind=rk), dimension(:), pointer :: sol, rhs  !< for testing solver
105  real(kind=rk), pointer :: rhs_1(:), g(:)
106
107  ! Partitioning
108  integer               :: nparts !< number of partitons to partition a mesh
109  integer, dimension(6) :: part_opts = (/0,0,0,0,0,0/)
110
111  type(ConvInf) :: resStat
112
113  real(kind=rk) :: t,t1
114  integer :: i,j,k,n,it
115  float(kind=rk), dimension(:), pointer :: xchk, r, y
116  character :: str
117  character(len=40) :: frm
118  float(kind=rk) :: strong_conn1,strong_conn2,cond_num,nrm
119  integer :: aggr_radius1,aggr_radius2
120  integer :: min_asize1,min_asize2
121  integer :: max_asize1,max_asize2
122  integer :: aver_finesize,min_finesize,max_finesize
123  integer :: aver_subdsize,min_subdsize,max_subdsize
124  integer :: start_radius1,start_radius2
125  integer :: plotting, ol
126  integer,allocatable :: nodes(:), inds(:)
127
128  type(AggrInfo) :: fAggr !< fine aggregates
129  type(AggrInfo) :: cAggr !< coarse aggregates
130  type(Decomposition) :: DD !< domain decomposition
131  type(RobustPreconditionMtx) :: C
132  type(CoarseSpace) :: CS
133  ! Parallel coarse level
134  !type(CoarseData) :: cdat -- moved into the module itself...
135
136  ! Init DOUG
137  call DOUG_Init()
138
139  ! Master participates in calculations as well
140  nparts = numprocs
141
142  call parallelDistributeInput(sctls%input_type,M,A,b,nparts,part_opts,A_ghost)
143
144  ! overlap for subdomains
145  if (sctls%overlap<0) then ! autom. overlap from smoothing
146    ol = max(sctls%smoothers,0)
147  else
148    ol = sctls%overlap
149  endif
150
151  if (sctls%levels>1.or.(numprocs==1.and.sctls%levels==1)) then !todo remove
152    ! Testing aggregation: AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
153    if (sctls%strong1/=0.0_rk) then
154      strong_conn1=sctls%strong1
155    else
156      strong_conn1=0.67_rk
157    endif
158    if (sctls%radius1>0) then
159      aggr_radius1=sctls%radius1
160    else
161      aggr_radius1=2
162    endif
163    if (sctls%minasize1>0) then
164      min_asize1=sctls%minasize1
165    else
166       ! Changes R. Scheichl 21/06/05
167      ! min_asize1=2*aggr_radius1+1
168       min_asize1=0.5_rk*(2*aggr_radius1+1)**2
169    endif
170    if (sctls%maxasize1>0) then
171      max_asize1=sctls%maxasize1
172    else
173      max_asize1=(2*aggr_radius1+1)**2
174    endif
175    if (numprocs>1) then
176      plotting=0
177    else
178      plotting=sctls%plotting
179    endif
180
181    ! find fine aggregates
182    if (numprocs > 1) then
183      ! we need to create aggregates only on inner nodes, so use local matrix LA
184      !  instead of expanded (to overlap) local matrix A
185      LA = getLocal(A,M)
186      call SpMtx_find_strong(A=LA,alpha=strong_conn1)
187      call SpMtx_aggregate(LA,fAggr,aggr_radius1, &
188           minaggrsize=min_asize1,       &
189           maxaggrsize=max_asize1,       &
190           alpha=strong_conn1,           &
191           M=M,                          &
192           plotting=plotting)
193      call SpMtx_unscale(LA)
194
195      write(stream,*) "_----------", fAggr%inner%nagr, size(fAggr%inner%starts)
196  else
197      ! non-parallel case use the whole matrix
198      call SpMtx_find_strong(A=A,alpha=strong_conn1)
199      call SpMtx_aggregate(A,fAggr,aggr_radius1, &
200            minaggrsize=min_asize1,       &
201            maxaggrsize=max_asize1,       &
202            alpha=strong_conn1,           &
203            M=M,                          &
204            plotting=plotting)
205      call SpMtx_unscale(A)
206      !call Aggrs_readFile_fine(A%aggr, "aggregates.txt")
207    end if
208    ! profile info
209    if(pstream/=0) then
210      write(pstream, "(I0,':fine aggregates:',I0)") myrank, fAggr%inner%nagr
211    end if
212
213    !if (sctls%plotting>=2) then
214    !   call SpMtx_writeLogicalValues(A, A%strong, 'strong.txt')
215    !end if
216
217    call Mesh_printInfo(M)
218   
219    if (numprocs==1.and.sctls%plotting==2.and.M%nell>0) then
220      call Mesh_pl2D_plotAggregate(fAggr%inner,M,&
221                      A%strong_rowstart,A%strong_colnrs,&
222                      mctls%assembled_mtx_file, &
223                                 INIT_CONT_END=D_PLPLOT_INIT)
224                                 !D_PLPLOT_END)
225    endif
226   
227    ! .. Testing aggregationAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
228   
229    ! Testing coarse matrix and aggregation through it:
230    if (numprocs>1) then
231      call SpMtx_find_strong(A=A,alpha=strong_conn1,A_ghost=A_ghost,M=M)
232      call SpMtx_unscale(A)
233      call IntRestBuild(A,fAggr%inner,Restrict,A_ghost)
234      CS = CoarseSpace_Init(Restrict)
235      call CoarseData_Copy(cdat,cdat_vec)
236      call CoarseSpace_Expand(CS,Restrict,M,cdat)
237      call CoarseMtxBuild(A,cdat%LAC,Restrict,M%ninner,A_ghost)
238      call KeepGivenRowIndeces(Restrict, (/(i,i=1,fAggr%inner%nagr)/))
239
240      if (sctls%verbose>3.and.cdat%LAC%nnz<400) then
241        write(stream,*)'Restrict (local) is:=================='
242        call SpMtx_printRaw(A=Restrict)
243        write(stream,*)'A coarse (local) is:=================='
244        call SpMtx_printRaw(A=cdat%LAC)
245      endif
246
247    else
248      ! non-parallel
249      call IntRestBuild(A,fAggr%inner,Restrict)
250
251      if (sctls%coarse_method<=1) then ! if not specified or ==1
252         call CoarseMtxBuild(A,AC,Restrict,M%ninner)
253
254      else if (sctls%coarse_method==2) then
255         ! use the Robust Coarse Spaces algorithm
256         B_RCS = CoarseProjectionMtxsBuild(A,Restrict,fAggr%inner%nagr)
257         allocate(rhs_1(A%nrows))
258         allocate(g(A%ncols))
259
260         rhs_1 = 1.0
261         call pcg_forRCS(B_RCS,rhs_1,g)
262
263         if(sctls%verbose>5) then
264            write (stream, *) "Solution g = ", g
265         end if
266
267         call RobustRestrictMtxBuild(B_RCS,g,Restrict)
268         call CoarseMtxBuild(A,AC,Restrict,M%ninner)
269
270      else
271         write(stream,'(A," ",I2)') 'Wrong coarse method', sctls%coarse_method
272         call DOUG_abort('Error in aggr', -1)
273      endif
274           
275      if (sctls%verbose>3.and.AC%nnz<400) then
276        write(stream,*)'A coarse is:=================='
277        call SpMtx_printRaw(A=AC)
278      endif
279    endif
280   
281    if (numprocs==1) then
282      if (sctls%strong2>0) then
283        strong_conn2=sctls%strong2
284      else
285        strong_conn2=strong_conn1/2.0_rk
286      endif
287      call SpMtx_find_strong(AC,strong_conn2)
288   
289      if (sctls%radius2>0) then
290        aggr_radius2=sctls%radius2
291      else
292         n=sqrt(1.0_rk*A%nrows)
293         aggr_radius2=nint(3*sqrt(dble(n))/(2*aggr_radius1+1)-1)
294         write (stream,*) 'Coarse aggregation radius aggr_radius2 =',aggr_radius2
295      endif
296      if (sctls%minasize2>0) then
297        min_asize2=sctls%minasize2
298      elseif (sctls%radius2>0) then
299        min_asize2=2*sctls%radius2+1
300      else
301        min_asize2=0.5_rk*(2*aggr_radius2+1)**2
302      endif
303      if (sctls%maxasize2>0) then
304        max_asize2=sctls%maxasize2
305      else
306        !max_asize2=max_asize1
307        max_asize2=(2*aggr_radius2+1)**2
308      endif
309      call SpMtx_aggregate(AC,cAggr,aggr_radius2, &
310            minaggrsize=min_asize2,          &
311            maxaggrsize=max_asize2,          &
312            alpha=strong_conn2,              &
313            aggr_fine=fAggr)
314      call SpMtx_unscale(AC)
315      !call Aggrs_readFile_coarse(cAggr, "aggregates.txt")
316
317      ! profile info
318      if(pstream/=0) then
319        write(pstream, "(I0,':coarse aggregates:',I0)") myrank, cAggr%inner%nagr
320      end if
321
322      if (sctls%plotting==2) then
323         call Aggr_writeFile(fAggr%inner, 'aggr2.txt', cAggr%inner)
324      end if
325      if (sctls%plotting==2.and.M%nell>0) then
326        !print *,'press Key<Enter>'
327        !read *,str
328        call Mesh_pl2D_plotAggregate(fAggr%inner,M,&
329                        A%strong_rowstart,A%strong_colnrs,&
330                        mctls%assembled_mtx_file, &
331                        caggrnum=cAggr%inner%num, &
332                      INIT_CONT_END=D_PLPLOT_END)!, &
333                      !INIT_CONT_END=D_PLPLOT_CONT)!, &
334                                  ! D_PLPLOT_END)
335      endif
336      write(stream,*)'# coarse aggregates:',cAggr%inner%nagr
337
338    endif
339  endif
340
341  if (numprocs==1) then
342    DD = Decomposition_from_aggrs(A, cAggr%full, fAggr%full, ol)
343    call AggrInfo_Destroy(cAggr)
344    call AggrInfo_Destroy(fAggr)
345  else
346    DD = Decomposition_full(A,A_ghost,M%ninner,ol)
347    ! call Aggrs_writeFile(M, fAggr, cdat, "aggregates.txt")
348    if (sctls%levels>1) call AggrInfo_Destroy(fAggr)
349  end if
350
351  ! Testing UMFPACK:
352  allocate(sol(A%nrows))
353  allocate(rhs(A%ncols))
354  ! rhs=1.0_rk
355
356  ! Solve the system
357! allocate(xl(A%nrows))
358! allocate(b(A%nrows))
359  allocate(xl(M%nlf))
360  xl = 0.0_rk
361  if (.FALSE..and.sctls%input_type==DCTL_INPUT_TYPE_ASSEMBLED.and.len_trim(mctls%assembled_rhs_file)<=0) then ! just test the solver
362    ! Set solution to random vector and calcluate RHS via b = A*x
363    write(stream,'(a,a)')' ##### (testing the solver: random RHS with known answer)'
364    allocate(xchk(M%nlf))
365    call random_number(xchk(1:M%ninner))
366    xchk(1:M%ninner) = 0.5_rk - xchk(1:M%ninner)
367    call update_outer_ol(xchk,M)
368    !call Print_Glob_Vect(xchk,M,'global xchk===')
369    call SpMtx_pmvm(b,A,xchk,M)
370    !call Print_Glob_Vect(b,M,'global b===')
371    !call MPI_BARRIER(MPI_COMM_WORLD,i)
372    ! and also normalise the rhs:
373    nrm=Vect_dot_product(b,b)
374    b=b/dsqrt(nrm)
375    xchk=xchk/dsqrt(nrm)
376    ! Another good test
377    !write(stream,'(a,a)') ' ##### (using unit vector as RHS) ##### '
378    !b=1.0_rk
379  endif
380
381  select case(sctls%solver)
382  case (DCTL_SOLVE_CG)
383     ! Conjugate gradient
384     call cg(A, b, xl, M, solinf=resStat, resvects_=.true.)
385     !call cg(A, b, xl, M, solinf=resStat)
386  case (DCTL_SOLVE_PCG)
387     ! Preconditioned conjugate gradient
388     !call pcg(A, b, xl, M, solinf=resStat, resvects_in=.true.)
389     t1 = MPI_WTIME()
390     if (numprocs==1) then
391       write(stream,*)'calling pcg_weigs /1/...'
392       call pcg_weigs(A=A,b=b,x=xl,Msh=M,DomDec=DD,it=it,cond_num=cond_num, &
393          CoarseMtx_=AC,Restrict=Restrict,refactor_=.true.)
394     else
395       !if (max(sctls%overlap,sctls%smoothers)>0) then
396       !  write(stream,*)'calling pcg_weigs /2/...'
397       !  call pcg_weigs(A=A,b=b,x=xl,Msh=M,it=it,cond_num=cond_num, &
398       !    A_interf_=A_ghost,refactor_=.true.)
399       !else
400         if (sctls%levels==2) then
401           write(stream,*)'calling pcg_weigs /3/...'
402           call pcg_weigs(A=A,b=b,x=xl,Msh=M,DomDec=DD,it=it,cond_num=cond_num, &
403                A_interf_=A_ghost, &
404                  CoarseMtx_=AC,Restrict=Restrict, &
405                  refactor_=.true.)
406         else
407           write(stream,*)'calling pcg_weigs /4/...'
408           call pcg_weigs(A, b, xl, M,DD,it,cond_num, &
409                A_interf_=A_ghost, refactor_=.true.)
410         endif
411       !endif
412     endif
413     t=MPI_WTIME()-t1
414     write(stream,*) 'time spent in pcg():',t
415     t1=total_setup_time()
416     write(stream,*) '    ...of which setup:',t1
417     write(stream,*) '       ...of which factorisation:', &
418        total_factorisation_time()
419     write(stream,*) 'solve time without setup:',t-t1
420  case default
421     call DOUG_abort('[DOUG main] : Wrong solution method specified', -1)
422  end select
423  if (associated(xchk)) then
424    ! Check the error
425    if (numprocs==1) then
426      allocate(r(A%nrows))
427    else
428      allocate(r(M%nlf))
429    endif
430    r = xl - xchk
431    write(stream,*) 'CHECK: The norm of the error is ', &
432           sqrt(Vect_dot_product(r,r)/Vect_dot_product(xchk,xchk))
433    deallocate(xchk)
434    deallocate(r)
435  endif
436
437  if (numprocs>1) then
438    ! Assemble result on master
439    if (ismaster()) then
440       print *, "freedoms", M%ngf
441      allocate(x(M%ngf)); x = 0.0_rk
442    end if
443    call Vect_Gather(xl, x, M)
444    if (ismaster().and.sctls%verbose>2.and.(size(x) <= 100)) &
445      call Vect_Print(x, 'sol > ')
446    if (ismaster()) then
447       call WriteSolutionToFile(x)
448    end if
449
450  else
451    if (sctls%verbose>2.and.size(xl)<=100) then
452      call Vect_Print(xl, 'sol > ')
453    endif
454    call WriteSolutionToFile(xl)
455  endif
456
457
458  ! Destroy objects
459  call Mesh_Destroy(M)
460  call SpMtx_Destroy(A)
461  call ConvInf_Destroy(resStat)
462  if (associated(b)) deallocate(b)
463  if (associated(xl)) deallocate(xl)
464  if (associated(x)) deallocate(x)
465  if (associated(sol)) deallocate(sol)
466
467  call DOUG_Finalize()
468
469
470end program main_aggr
Note: See TracBrowser for help on using the repository browser.