From e1c05d0e228383fcb140d91d1b6ad590816e23e1 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 30 Sep 2016 13:59:41 +0000 Subject: [PATCH] mld2p4-extaggr: Make.inc.in Makefile mlprec/impl/Makefile mlprec/impl/mld_d_bld_mlhier_array.f90 mlprec/impl/mld_d_dec_map_bld.f90 mlprec/impl/mld_d_hierarchy_bld.f90 mlprec/impl/mld_d_lev_aggrmap_bld.f90 mlprec/impl/mld_d_lev_aggrmat_asb.f90 mlprec/impl/mld_daggrmap_bld.f90 mlprec/impl/mld_daggrmat_asb.f90 mlprec/impl/mld_daggrmat_biz_asb.f90 mlprec/impl/mld_daggrmat_minnrg_asb.f90 mlprec/impl/mld_daggrmat_nosmth_asb.f90 mlprec/impl/mld_daggrmat_smth_asb.f90 mlprec/impl/mld_dcoarse_bld.f90 mlprec/impl/mld_dprecinit.F90 mlprec/mld_base_prec_type.F90 mlprec/mld_d_inner_mod.f90 tests/pdegen/mld_d_pde3d.f90 tests/pdegen/runs/ppde.inp Refactored map bld/mat asb routines. Now there is no need to call coarse_bld twice. To be cleaned up yet. --- Make.inc.in | 1 + Makefile | 4 + mlprec/impl/Makefile | 2 +- mlprec/impl/mld_d_bld_mlhier_array.f90 | 54 ++-- mlprec/impl/mld_d_dec_map_bld.f90 | 53 ++-- mlprec/impl/mld_d_hierarchy_bld.f90 | 353 +++++++++++++++++++++++- mlprec/impl/mld_d_lev_aggrmap_bld.f90 | 148 ++++++++++ mlprec/impl/mld_d_lev_aggrmat_asb.f90 | 252 +++++++++++++++++ mlprec/impl/mld_daggrmap_bld.f90 | 31 ++- mlprec/impl/mld_daggrmat_asb.f90 | 110 -------- mlprec/impl/mld_daggrmat_biz_asb.f90 | 25 +- mlprec/impl/mld_daggrmat_minnrg_asb.f90 | 23 +- mlprec/impl/mld_daggrmat_nosmth_asb.f90 | 31 +-- mlprec/impl/mld_daggrmat_smth_asb.f90 | 21 +- mlprec/impl/mld_dcoarse_bld.f90 | 322 +++++++++++---------- mlprec/impl/mld_dprecinit.F90 | 2 +- mlprec/mld_base_prec_type.F90 | 17 +- mlprec/mld_d_inner_mod.f90 | 49 ++-- tests/pdegen/mld_d_pde3d.f90 | 19 +- tests/pdegen/runs/ppde.inp | 22 +- 20 files changed, 1104 insertions(+), 435 deletions(-) create mode 100644 mlprec/impl/mld_d_lev_aggrmap_bld.f90 create mode 100644 mlprec/impl/mld_d_lev_aggrmat_asb.f90 diff --git a/Make.inc.in b/Make.inc.in index f03d5958..a89030d5 100644 --- a/Make.inc.in +++ b/Make.inc.in @@ -15,6 +15,7 @@ INSTALL_DIR=@INSTALL_DIR@ INSTALL_LIBDIR=@INSTALL_LIBDIR@ INSTALL_INCLUDEDIR=@INSTALL_INCLUDEDIR@ INSTALL_DOCSDIR=@INSTALL_DOCSDIR@ +INSTALL_SAMPLESDIR=@INSTALL_SAMPLESDIR@ ########################################################## diff --git a/Makefile b/Makefile index 5e66f185..413bb307 100644 --- a/Makefile +++ b/Makefile @@ -29,6 +29,10 @@ install: all /bin/cp -fr docs/*pdf docs/html $(INSTALL_DOCSDIR)) (./mkdir.sh $(INSTALL_DOCSDIR) && \ $(INSTALL_DATA) README LICENSE $(INSTALL_DOCSDIR)) + (./mkdir.sh $(INSTALL_SAMPLESDIR) && ./mkdir.sh $(INSTALL_SAMPLESDIR)/simple &&\ + ./mkdir.sh $(INSTALL_SAMPLESDIR)/advanced && \ + (cd examples; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/simple ) && \ + (cd tests; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/advanced )) cleanlib: (cd lib; /bin/rm -f *.a *$(.mod) *$(.fh)) (cd include; /bin/rm -f *.a *$(.mod) *$(.fh)) diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index c9c0e713..3db4d32f 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -26,7 +26,7 @@ DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o mld_d_bld_mlhier_aggsize.o mld mld_d_ml_prec_bld.o mld_d_hierarchy_bld.o \ mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \ mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \ - $(DMPFOBJS) mld_d_extprol_bld.o + $(DMPFOBJS) mld_d_extprol_bld.o mld_d_lev_aggrmap_bld.o mld_d_lev_aggrmat_asb.o SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o mld_s_bld_mlhier_aggsize.o mld_s_bld_mlhier_array.o \ mld_s_ml_prec_bld.o mld_s_hierarchy_bld.o \ diff --git a/mlprec/impl/mld_d_bld_mlhier_array.f90 b/mlprec/impl/mld_d_bld_mlhier_array.f90 index 3618de1b..bd0a0e77 100644 --- a/mlprec/impl/mld_d_bld_mlhier_array.f90 +++ b/mlprec/impl/mld_d_bld_mlhier_array.f90 @@ -37,23 +37,25 @@ !!$ !!$ -subroutine mld_d_bld_mlhier_array(nplevs,a,desc_a,precv,info) +subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) use psb_base_mod use mld_d_inner_mod, mld_protect_name => mld_d_bld_mlhier_array use mld_d_prec_mod implicit none - integer(psb_ipk_), intent(inout) :: nplevs + integer(psb_ipk_), intent(inout) :: nplevs, casize + real(psb_dpk_) :: mnaggratio type(psb_dspmat_type),intent(in), target :: a type(psb_desc_type), intent(inout), target :: desc_a type(mld_d_onelev_type),intent(inout), allocatable, target :: precv(:) integer(psb_ipk_), intent(out) :: info ! Local integer(psb_ipk_) :: ictxt, me,np - integer(psb_ipk_) :: err,i,k, err_act, newsz, iszv + integer(psb_ipk_) :: err,i,k, err_act, newsz, iszv, iaggsize integer(psb_ipk_) :: ipv(mld_ifpsz_), val class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm type(mld_dml_parms) :: baseparms, medparms, coarseparms type(mld_d_onelev_type), allocatable :: tprecv(:) + real(psb_dpk_) :: sizeratio integer(psb_ipk_) :: int_err(5) integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name, ch_err @@ -178,27 +180,45 @@ subroutine mld_d_bld_mlhier_array(nplevs,a,desc_a,precv,info) & write(debug_unit,*) me,' ',trim(name),& & 'Return from ',i,' call to mlprcbld ',info - if (i>2) then + iaggsize = sum(precv(i)%map%naggr) + if (iaggsize <= casize) then + newsz = i + end if + + if (i>2) then + sizeratio = iaggsize + sizeratio = sum(precv(i-1)%map%naggr)/sizeratio + if (sizeratio < mnaggratio) then + if (sizeratio > 1) then + newsz = i + else + ! + ! We are not gaining + ! + newsz = newsz-1 + end if + end if + if (all(precv(i)%map%naggr == precv(i-1)%map%naggr)) then newsz=i-1 + if (me == 0) then + write(debug_unit,*) trim(name),& + &': Warning: aggregates from level ',& + & newsz + write(debug_unit,*) trim(name),& + &': to level ',& + & iszv,' coincide.' + write(debug_unit,*) trim(name),& + &': Number of levels actually used :',newsz + write(debug_unit,*) + end if end if - call psb_bcast(ictxt,newsz) - if (newsz > 0) exit array_build_loop end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) exit array_build_loop end do array_build_loop if (newsz > 0) then - if (me == 0) then - write(debug_unit,*) trim(name),& - &': Warning: aggregates from level ',& - & newsz - write(debug_unit,*) trim(name),& - &': to level ',& - & iszv,' coincide.' - write(debug_unit,*) trim(name),& - &': Number of levels actually used :',newsz - write(debug_unit,*) - end if allocate(tprecv(newsz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& diff --git a/mlprec/impl/mld_d_dec_map_bld.f90 b/mlprec/impl/mld_d_dec_map_bld.f90 index fe263c2a..363aeec0 100644 --- a/mlprec/impl/mld_d_dec_map_bld.f90 +++ b/mlprec/impl/mld_d_dec_map_bld.f90 @@ -54,12 +54,12 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! Local variables integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& - & ideg(:), idxs(:) + & ideg(:), idxs(:), tmpaggr(:) real(psb_dpk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip, nisolate,nthird type(psb_d_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl - logical :: recovery, candidate + logical :: recovery, disjoint integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -114,6 +114,7 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! naggr = 0 icnt = 0 + nisolate = 0 step1: do ii=1, nr i = idxs(ii) @@ -138,20 +139,20 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if end if enddo - if (ip < 1) then - write(0,*) "Should at least contain the node itself ! " - cycle step1 - end if - candidate = .true. - do k=1, ip - candidate = candidate .and. (ilaggr(icol(k)) == -(nr+1)) - end do - if (candidate) then - ! - ! If the whole strongly coupled neighborhood of I is - ! as yet unconnected, turn it into the next aggregate - ! +!!$ if (ip <= 1) then +!!$ nisolate = nisolate + 1 +!!$ cycle step1 +!!$ end if + + ! + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate. + ! Same if ip==0 (in which case, neighborhood only + ! contains I even if it does not look from matrix) + ! + disjoint = all(ilaggr(icol(1:ip)) == -(nr+1)).or.(ip==0) + if (disjoint) then icnt = icnt + 1 naggr = naggr + 1 do k=1, ip @@ -161,6 +162,7 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if endif enddo step1 + if (debug_level >= psb_debug_outer_) then write(debug_unit,*) me,' ',trim(name),& & ' Check 1:',count(ilaggr == -(nr+1)) @@ -168,7 +170,8 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ! Phase two: join the neighbours - ! + ! + tmpaggr = ilaggr step2: do ii=1,nr i = idxs(ii) @@ -188,8 +191,9 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) do k=1, nz j = icol(k) if ((1<=j).and.(j<=nr)) then - if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& - & .and. (ilaggr(j) > 0).and. (abs(val(k)) > cpling)) then +!!$ if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& +!!$ & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then + if ((tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then ip = k cpling = abs(val(k)) end if @@ -205,10 +209,12 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) ! ! Phase three: sweep over leftovers, if any ! + nthird = 0 step3: do ii=1,nr i = idxs(ii) - if (ilaggr(i) == -(nr+1)) then + if (ilaggr(i) < 0) then + nthird = nthird + 1 call a%csget(i,i,nz,irow,icol,val,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -216,8 +222,8 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if ! - ! Find the most strongly connected neighbour that is - ! already aggregated, if any, and join its aggregate + ! Find its strongly connected neighbourhood not + ! already aggregated, and make it into a new aggregate. ! cpling = dzero ip = 0 @@ -282,7 +288,8 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) nlaggr(:) = 0 nlaggr(me+1) = naggr call psb_sum(ictxt,nlaggr(1:np)) - + + !write(*,*) me,'Info from dec_map_bld: ',naggr,nisolate,nthird call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/mld_d_hierarchy_bld.f90 b/mlprec/impl/mld_d_hierarchy_bld.f90 index c43238dd..7bf307bd 100644 --- a/mlprec/impl/mld_d_hierarchy_bld.f90 +++ b/mlprec/impl/mld_d_hierarchy_bld.f90 @@ -94,9 +94,13 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold) ! Local Variables integer(psb_ipk_) :: ictxt, me,np - integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs - real(psb_dpk_) :: mnaggratio - integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize + real(psb_dpk_) :: mnaggratio, sizeratio + class(mld_d_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_dml_parms) :: baseparms, medparms, coarseparms + integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_dspmat_type) :: op_prol + type(mld_d_onelev_type), allocatable :: tprecv(:) integer(psb_ipk_) :: int_err(5) character :: upd_ integer(psb_ipk_) :: debug_level, debug_unit @@ -191,10 +195,23 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold) goto 9999 endif + ! + ! The strategy: + ! 1. The maximum number of levels should be already encoded in the + ! size of the array; + ! 2. If the user did not specify anything, then a default coarse size + ! is generated, and the number of levels is set to the maximum; + ! 3. If the number of levels has been specified, make sure it's capped + ! at the maximum; + ! 4. If the size of the array is different from target number of levels, + ! reallocate; + ! 5. Build the matrix hierarchy, stopping early if either the target + ! coarse size is hit, or the gain fall below the mmin_aggr_ratio + ! threshold. + ! + + if (nplevs <= 0) then - ! - ! This should become the default strategy, we specify a target aggregation size. - ! if (casize <=0) then ! ! Default to the cubic root of the size at base level. @@ -204,15 +221,321 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold) casize = max(casize,ione) casize = casize*40_psb_ipk_ end if - call mld_bld_mlhier_aggsize(casize,mxplevs,mnaggratio,a,desc_a,p%precv,info) - else - ! - ! Oldstyle with fixed number of levels. - ! - nplevs = max(itwo,min(nplevs,mxplevs)) - call mld_bld_mlhier_array(nplevs,a,desc_a,p%precv,info) + nplevs = mxplevs end if - + + nplevs = max(itwo,min(nplevs,mxplevs)) + + coarseparms = p%precv(iszv)%parms + baseparms = p%precv(1)%parms + medparms = p%precv(2)%parms + + allocate(coarse_sm, source=p%precv(iszv)%sm,stat=info) + if (info == psb_success_) & + & allocate(med_sm, source=p%precv(2)%sm,stat=info) + if (info == psb_success_) & + & allocate(base_sm, source=p%precv(1)%sm,stat=info) + if (info /= psb_success_) then + write(0,*) 'Error in saving smoothers',info + call psb_errpush(psb_err_internal_error_,name,a_err='Base level precbuild.') + goto 9999 + end if + + ! + ! First set desired number of levels + ! + if (iszv /= nplevs) then + allocate(tprecv(nplevs),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + tprecv(1)%parms = baseparms + allocate(tprecv(1)%sm,source=base_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=2,nplevs-1 + tprecv(i)%parms = medparms + allocate(tprecv(i)%sm,source=med_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + end do + tprecv(nplevs)%parms = coarseparms + allocate(tprecv(nplevs)%sm,source=coarse_sm,stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,iszv + call p%precv(i)%free(info) + end do + call move_alloc(tprecv,p%precv) + iszv = size(p%precv) + end if + + ! + ! Finest level first; remember to fix base_a and base_desc + ! + p%precv(1)%base_a => a + p%precv(1)%base_desc => desc_a + newsz = 0 + if (.false.) then + old_array_build_loop: do i=2, iszv + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,p%precv(i)%parms) + + ! + ! Sanity checks on the parameters + ! + if (i= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + ! + ! Build the mapping between levels i-1 and i and the matrix + ! at level i + ! + if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& + & p%precv(i-1)%base_desc, p%precv(i),info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Init upper level preconditioner') + goto 9999 + endif + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Return from ',i,' call to mlprcbld ',info + + iaggsize = sum(p%precv(i)%map%naggr) + if (iaggsize <= casize) then + newsz = i + end if + + if (i>2) then + sizeratio = iaggsize + sizeratio = sum(p%precv(i-1)%map%naggr)/sizeratio + if (sizeratio < mnaggratio) then + if (sizeratio > 1) then + newsz = i + else + ! + ! We are not gaining + ! + newsz = newsz-1 + end if + end if + + if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then + newsz=i-1 + if (me == 0) then + write(debug_unit,*) trim(name),& + &': Warning: aggregates from level ',& + & newsz + write(debug_unit,*) trim(name),& + &': to level ',& + & iszv,' coincide.' + write(debug_unit,*) trim(name),& + &': Number of levels actually used :',newsz + write(debug_unit,*) + end if + end if + end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) exit old_array_build_loop + end do old_array_build_loop + + if (newsz > 0) then + ! + ! We exited early from the build loop, need to fix + ! the size. + ! + allocate(tprecv(newsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,newsz-1 + call p%precv(i)%move_alloc(tprecv(i),info) + end do + call p%precv(iszv)%move_alloc(tprecv(newsz),info) + do i=newsz+1, iszv + call p%precv(i)%free(info) + end do + call move_alloc(tprecv,p%precv) + ! Ignore errors from transfer + info = psb_success_ + ! + ! Restart + iszv = newsz + ! Fix the pointers, but the level 1 should + ! be already OK + do i=2, iszv - 1 + p%precv(i)%base_a => p%precv(i)%ac + p%precv(i)%base_desc => p%precv(i)%desc_ac + p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc + p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc + end do + + i = iszv + if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& + & p%precv(i-1)%base_desc, p%precv(i),info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') + goto 9999 + endif + end if + + else + + array_build_loop: do i=2, iszv + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,p%precv(i)%parms) + + ! + ! Sanity checks on the parameters + ! + if (i= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + ! + ! Build the mapping between levels i-1 and i and the matrix + ! at level i + ! + if (info == psb_success_) call mld_aggrmap_bld(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Map build') + goto 9999 + endif + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Return from ',i,' call to mlprcbld ',info + + + ! + ! Check for early termination of aggregation loop. + ! + iaggsize = sum(nlaggr) + if (iaggsize <= casize) then + newsz = i + end if + + if (i>2) then + sizeratio = iaggsize + sizeratio = sum(p%precv(i-1)%map%naggr)/sizeratio + if (sizeratio < mnaggratio) then + if (sizeratio > 1) then + newsz = i + else + ! + ! We are not gaining + ! + newsz = newsz-1 + end if + end if + + if (all(nlaggr == p%precv(i-1)%map%naggr)) then + newsz=i-1 + if (me == 0) then + write(debug_unit,*) trim(name),& + &': Warning: aggregates from level ',& + & newsz + write(debug_unit,*) trim(name),& + &': to level ',& + & iszv,' coincide.' + write(debug_unit,*) trim(name),& + &': Number of levels actually used :',newsz + write(debug_unit,*) + end if + end if + end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) & + & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) + + if (info == psb_success_) call mld_aggrmat_asb(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Map build') + goto 9999 + endif + + if (newsz > 0) exit array_build_loop + end do array_build_loop + + if (newsz > 0) then + ! + ! We exited early from the build loop, need to fix + ! the size. + ! + allocate(tprecv(newsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,& + & a_err='prec reallocation') + goto 9999 + endif + do i=1,newsz + call p%precv(i)%move_alloc(tprecv(i),info) + end do + do i=newsz+1, iszv + call p%precv(i)%free(info) + end do + call move_alloc(tprecv,p%precv) + ! Ignore errors from transfer + info = psb_success_ + ! + ! Restart + iszv = newsz + ! Fix the pointers, but the level 1 should + ! be already OK + do i=2, iszv + p%precv(i)%base_a => p%precv(i)%ac + p%precv(i)%base_desc => p%precv(i)%desc_ac + p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc + p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc + end do + end if + end if + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -220,7 +543,7 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold) endif iszv = size(p%precv) - + if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Exiting with',iszv,' levels' diff --git a/mlprec/impl/mld_d_lev_aggrmap_bld.f90 b/mlprec/impl/mld_d_lev_aggrmap_bld.f90 new file mode 100644 index 00000000..eec9a32c --- /dev/null +++ b/mlprec/impl/mld_d_lev_aggrmap_bld.f90 @@ -0,0 +1,148 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_dcoarse_bld.f90 +! +! Subroutine: mld_dcoarse_bld +! Version: real +! +! This routine builds the matrix associated to the current level of the +! multilevel preconditioner from the matrix associated to the previous level, +! by using a smoothed aggregation technique (therefore, it also builds the +! prolongation and restriction operators mapping the current level to the +! previous one and vice versa). Then the routine builds the base preconditioner +! at the current level. +! The current level is regarded as the coarse one, while the previous as +! the fine one. This is in agreement with the fact that the routine is called, +! by mld_mlprec_bld, only on levels >=2. +! +! +! Arguments: +! a - type(psb_dspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_d_onelev_type), input/output. +! The 'one-level' data structure containing the local part +! of the base preconditioner to be built as well as +! information concerning the prolongator and its transpose. +! info - integer, output. +! Error code. +! +subroutine mld_d_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_lev_aggrmap_bld + + implicit none + + ! Arguments + type(mld_d_onelev_type), intent(inout), target :: p + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: nzl, ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_d_lev_aggrmap_bld' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(p%parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(p%parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(p%parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + + select case(p%parms%aggr_alg) + case (mld_dec_aggr_, mld_sym_dec_aggr_) + + ! + ! Build a mapping between the row indices of the fine-level matrix + ! and the row indices of the coarse-level matrix, according to a decoupled + ! aggregation algorithm. This also defines a tentative prolongator from + ! the coarse to the fine level. + ! + call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld') + goto 9999 + end if + + case (mld_bcmatch_aggr_) + write(0,*) 'Matching is not implemented yet ' + info = -1111 + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/)) + goto 9999 + + case default + + info = -1 + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/)) + goto 9999 + + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_d_lev_aggrmap_bld diff --git a/mlprec/impl/mld_d_lev_aggrmat_asb.f90 b/mlprec/impl/mld_d_lev_aggrmat_asb.f90 new file mode 100644 index 00000000..ca47cc2e --- /dev/null +++ b/mlprec/impl/mld_d_lev_aggrmat_asb.f90 @@ -0,0 +1,252 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: mld_dcoarse_bld.f90 +! +! Subroutine: mld_dcoarse_bld +! Version: real +! +! This routine builds the matrix associated to the current level of the +! multilevel preconditioner from the matrix associated to the previous level, +! by using a smoothed aggregation technique (therefore, it also builds the +! prolongation and restriction operators mapping the current level to the +! previous one and vice versa). Then the routine builds the base preconditioner +! at the current level. +! The current level is regarded as the coarse one, while the previous as +! the fine one. This is in agreement with the fact that the routine is called, +! by mld_mlprec_bld, only on levels >=2. +! +! +! Arguments: +! a - type(psb_dspmat_type). +! The sparse matrix structure containing the local part of the +! fine-level matrix. +! desc_a - type(psb_desc_type), input. +! The communication descriptor of a. +! p - type(mld_d_onelev_type), input/output. +! The 'one-level' data structure containing the local part +! of the base preconditioner to be built as well as +! information concerning the prolongator and its transpose. +! info - integer, output. +! Error code. +! +subroutine mld_d_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_lev_aggrmat_asb + + implicit none + + ! Arguments + type(mld_d_onelev_type), intent(inout), target :: p + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(inout) :: op_prol + integer(psb_ipk_), intent(out) :: info + + + ! Local variables + character(len=20) :: name + integer(psb_mpik_) :: ictxt, np, me + integer(psb_ipk_) :: err_act + type(psb_dspmat_type) :: ac, op_restr + type(psb_d_coo_sparse_mat) :: acoo, bcoo + type(psb_d_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl, ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_dcoarse_bld' + if (psb_get_errstatus().ne.0) return + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + info = psb_success_ + ictxt = desc_a%get_context() + call psb_info(ictxt,me,np) + + call mld_check_def(p%parms%aggr_kind,'Smoother',& + & mld_smooth_prol_,is_legal_ml_aggr_kind) + call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + call mld_check_def(p%parms%aggr_filter,'Use filtered matrix',& + & mld_no_filter_mat_,is_legal_aggr_filter) + call mld_check_def(p%parms%smoother_pos,'smooth_pos',& + & mld_pre_smooth_,is_legal_ml_smooth_pos) + call mld_check_def(p%parms%aggr_omega_alg,'Omega Alg.',& + & mld_eig_est_,is_legal_ml_aggr_omega_alg) + call mld_check_def(p%parms%aggr_eig,'Eigenvalue estimate',& + & mld_max_norm_,is_legal_ml_aggr_eig) + call mld_check_def(p%parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega) + + + ! + ! Build the coarse-level matrix from the fine-level one, starting from + ! the mapping defined by mld_aggrmap_bld and applying the aggregation + ! algorithm specified by p%iprcparm(mld_aggr_kind_) + ! + call mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p%parms,ac,op_prol,op_restr,info) + + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') + goto 9999 + end if + + + ! Common code refactored here. + + ntaggr = sum(nlaggr) + + select case(p%parms%coarse_mat) + + case(mld_distr_mat_) + + call ac%mv_to(bcoo) + if (p%parms%clean_zeros) call bcoo%clean_zeros(info) + nzl = bcoo%get_nzeros() + + if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) + if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') + if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Assembld aux descr. distr.' + call p%ac%mv_from(bcoo) + + call p%ac%set_nrows(p%desc_ac%get_local_rows()) + call p%ac%set_ncols(p%desc_ac%get_local_cols()) + call p%ac%set_asb() + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') + goto 9999 + end if + + if (np>1) then + call op_prol%mv_to(acsr1) + nzl = acsr1%get_nzeros() + call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + goto 9999 + end if + call op_prol%mv_from(acsr1) + endif + call op_prol%set_ncols(p%desc_ac%get_local_cols()) + + if (np>1) then + call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) + call op_restr%mv_to(acoo) + nzl = acoo%get_nzeros() + if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') + call acoo%set_dupl(psb_dupl_add_) + if (info == psb_success_) call op_restr%mv_from(acoo) + if (info == psb_success_) call op_restr%cscnv(info,type='csr') + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Converting op_restr to local') + goto 9999 + end if + end if + ! + ! Clip to local rows. + ! + call op_restr%set_nrows(p%desc_ac%get_local_rows()) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' + + case(mld_repl_mat_) + ! + ! + call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) + if (info == psb_success_) & + & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) + + if (info /= psb_success_) goto 9999 + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select + + call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') + goto 9999 + end if + + ! + ! Copy the prolongation/restriction matrices into the descriptor map. + ! op_restr => PR^T i.e. restriction operator + ! op_prol => PR i.e. prolongation operator + ! + + p%map = psb_linmap(psb_map_aggr_,desc_a,& + & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) + if (info == psb_success_) call op_prol%free() + if (info == psb_success_) call op_restr%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + goto 9999 + end if + ! + ! Fix the base_a and base_desc pointers for handling of residuals. + ! This is correct because this routine is only called at levels >=2. + ! + p%base_a => p%ac + p%base_desc => p%desc_ac + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine mld_d_lev_aggrmat_asb diff --git a/mlprec/impl/mld_daggrmap_bld.f90 b/mlprec/impl/mld_daggrmap_bld.f90 index c652983a..de053aa6 100644 --- a/mlprec/impl/mld_daggrmap_bld.f90 +++ b/mlprec/impl/mld_daggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_d_inner_mod, mld_protect_name => mld_daggrmap_bld @@ -92,14 +92,15 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) real(psb_dpk_), intent(in) :: theta type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(out) :: op_prol integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer(psb_ipk_), intent(out) :: info ! Local variables integer(psb_ipk_), allocatable :: ils(:), neigh(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m,naggrm1, naggrp1, ntaggr type(psb_dspmat_type) :: atmp, atrans - logical :: recovery + type(psb_d_coo_sparse_mat) :: tmpcoo integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -151,6 +152,30 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) goto 9999 end if + naggr = nlaggr(me+1) + ntaggr = sum(nlaggr) + + naggrm1 = sum(nlaggr(1:me)) + naggrp1 = sum(nlaggr(1:me+1)) + ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 + call psb_halo(ilaggr,desc_a,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') + goto 9999 + end if + + call tmpcoo%allocate(ncol,ntaggr,ncol) + do i=1,ncol + tmpcoo%val(i) = done + tmpcoo%ia(i) = i + tmpcoo%ja(i) = ilaggr(i) + end do + call tmpcoo%set_nzeros(ncol) + call tmpcoo%set_dupl(psb_dupl_add_) + call tmpcoo%set_sorted() ! At this point this is in row-major + call op_prol%mv_from(tmpcoo) + call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/mld_daggrmat_asb.f90 b/mlprec/impl/mld_daggrmat_asb.f90 index 69965097..dd8dd9e4 100644 --- a/mlprec/impl/mld_daggrmat_asb.f90 +++ b/mlprec/impl/mld_daggrmat_asb.f90 @@ -166,116 +166,6 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,inf goto 9999 end if - -!!$ -!!$ ntaggr = sum(nlaggr) -!!$ -!!$ select case(p%parms%coarse_mat) -!!$ -!!$ case(mld_distr_mat_) -!!$ -!!$ call ac%mv_to(bcoo) -!!$ if (p%parms%clean_zeros) call bcoo%clean_zeros(info) -!!$ nzl = bcoo%get_nzeros() -!!$ -!!$ if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) -!!$ if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) -!!$ if (info == psb_success_) call psb_cdasb(p%desc_ac,info) -!!$ if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') -!!$ if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') -!!$ if (info /= psb_success_) then -!!$ call psb_errpush(psb_err_internal_error_,name,& -!!$ & a_err='Creating p%desc_ac and converting ac') -!!$ goto 9999 -!!$ end if -!!$ if (debug_level >= psb_debug_outer_) & -!!$ & write(debug_unit,*) me,' ',trim(name),& -!!$ & 'Assembld aux descr. distr.' -!!$ call p%ac%mv_from(bcoo) -!!$ -!!$ call p%ac%set_nrows(p%desc_ac%get_local_rows()) -!!$ call p%ac%set_ncols(p%desc_ac%get_local_cols()) -!!$ call p%ac%set_asb() -!!$ -!!$ if (info /= psb_success_) then -!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') -!!$ goto 9999 -!!$ end if -!!$ -!!$ if (np>1) then -!!$ call op_prol%mv_to(acsr1) -!!$ nzl = acsr1%get_nzeros() -!!$ call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') -!!$ if(info /= psb_success_) then -!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') -!!$ goto 9999 -!!$ end if -!!$ call op_prol%mv_from(acsr1) -!!$ endif -!!$ call op_prol%set_ncols(p%desc_ac%get_local_cols()) -!!$ -!!$ if (np>1) then -!!$ call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) -!!$ call op_restr%mv_to(acoo) -!!$ nzl = acoo%get_nzeros() -!!$ if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') -!!$ call acoo%set_dupl(psb_dupl_add_) -!!$ if (info == psb_success_) call op_restr%mv_from(acoo) -!!$ if (info == psb_success_) call op_restr%cscnv(info,type='csr') -!!$ if(info /= psb_success_) then -!!$ call psb_errpush(psb_err_internal_error_,name,& -!!$ & a_err='Converting op_restr to local') -!!$ goto 9999 -!!$ end if -!!$ end if -!!$ ! -!!$ ! Clip to local rows. -!!$ ! -!!$ call op_restr%set_nrows(p%desc_ac%get_local_rows()) -!!$ -!!$ if (debug_level >= psb_debug_outer_) & -!!$ & write(debug_unit,*) me,' ',trim(name),& -!!$ & 'Done ac ' -!!$ -!!$ case(mld_repl_mat_) -!!$ ! -!!$ ! -!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) -!!$ if (info == psb_success_) call psb_cdasb(p%desc_ac,info) -!!$ if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) -!!$ if (info == psb_success_) & -!!$ & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) -!!$ -!!$ if (info /= psb_success_) goto 9999 -!!$ -!!$ case default -!!$ info = psb_err_internal_error_ -!!$ call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') -!!$ goto 9999 -!!$ end select -!!$ -!!$ call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) -!!$ if(info /= psb_success_) then -!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') -!!$ goto 9999 -!!$ end if -!!$ -!!$ ! -!!$ ! Copy the prolongation/restriction matrices into the descriptor map. -!!$ ! op_restr => PR^T i.e. restriction operator -!!$ ! op_prol => PR i.e. prolongation operator -!!$ ! -!!$ -!!$ p%map = psb_linmap(psb_map_aggr_,desc_a,& -!!$ & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) -!!$ if (info == psb_success_) call op_prol%free() -!!$ if (info == psb_success_) call op_restr%free() -!!$ if(info /= psb_success_) then -!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') -!!$ goto 9999 -!!$ end if - - call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/mld_daggrmat_biz_asb.f90 b/mlprec/impl/mld_daggrmat_biz_asb.f90 index 90f95927..9d8dfa79 100644 --- a/mlprec/impl/mld_daggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_daggrmat_biz_asb.f90 @@ -89,7 +89,8 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) type(mld_dml_parms), intent(inout) :: parms - type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -128,19 +129,8 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr naggr = nlaggr(me+1) ntaggr = sum(nlaggr) - naggrm1 = sum(nlaggr(1:me)) - naggrp1 = sum(nlaggr(1:me+1)) - filter_mat = (parms%aggr_filter == mld_filter_mat_) - ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 - call psb_halo(ilaggr,desc_a,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') - goto 9999 - end if - ! naggr: number of local aggregates ! nrow: local rows. ! @@ -157,17 +147,10 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr end if ! 1. Allocate Ptilde in sparse matrix form - call tmpcoo%allocate(ncol,naggr,ncol) - do i=1,nrow - tmpcoo%val(i) = done - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) - end do - call tmpcoo%set_nzeros(nrow) - call tmpcoo%set_dupl(psb_dupl_add_) - + call op_prol%mv_to(tmpcoo) call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) + if (info /= psb_success_) goto 9999 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 index 09e7bad6..a08ab607 100644 --- a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 @@ -109,7 +109,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) type(mld_dml_parms), intent(inout) :: parms - type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -167,14 +168,6 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re filter_mat = (parms%aggr_filter == mld_filter_mat_) - ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 - call psb_halo(ilaggr,desc_a,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') - goto 9999 - end if - ! naggr: number of local aggregates ! nrow: local rows. ! @@ -209,20 +202,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! 1. Allocate Ptilde in sparse matrix form - call tmpcoo%allocate(ncol,ntaggr,ncol) - do i=1,ncol - tmpcoo%val(i) = done - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) - end do - call tmpcoo%set_nzeros(ncol) - call tmpcoo%set_dupl(psb_dupl_add_) - call tmpcoo%set_asb() + call op_prol%mv_to(tmpcoo) call ptilde%mv_from(tmpcoo) call ptilde%cscnv(info,type='csr') -!!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') - if (info == psb_success_) call a%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(da,info,type='csr',dupl=psb_dupl_add_) if (info /= psb_success_) then diff --git a/mlprec/impl/mld_daggrmat_nosmth_asb.f90 b/mlprec/impl/mld_daggrmat_nosmth_asb.f90 index 9fbec777..d4788037 100644 --- a/mlprec/impl/mld_daggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_daggrmat_nosmth_asb.f90 @@ -92,7 +92,8 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) type(mld_dml_parms), intent(inout) :: parms - type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -124,34 +125,12 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ntaggr = sum(nlaggr) naggrm1=sum(nlaggr(1:me)) - do i=1, nrow - ilaggr(i) = ilaggr(i) + naggrm1 - end do - call psb_halo(ilaggr,desc_a,info) - - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') - goto 9999 - end if - call acoo%allocate(ncol,ntaggr,ncol) - do i=1,nrow - acoo%val(i) = done - acoo%ia(i) = i - acoo%ja(i) = ilaggr(i) - end do - - call acoo%set_dupl(psb_dupl_add_) - call acoo%set_nzeros(nrow) - call acoo%set_asb() - call acoo%fix(info) - - - call op_prol%mv_from(acoo) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call op_prol%transp(op_restr) - + if (info /= psb_success_) goto 9999 + call op_prol%transp(op_restr) + call a%cp_to(ac_coo) nzt = ac_coo%get_nzeros() diff --git a/mlprec/impl/mld_daggrmat_smth_asb.f90 b/mlprec/impl/mld_daggrmat_smth_asb.f90 index db400747..46b34133 100644 --- a/mlprec/impl/mld_daggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_daggrmat_smth_asb.f90 @@ -104,7 +104,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) type(mld_dml_parms), intent(inout) :: parms - type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -147,13 +148,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest naggrp1 = sum(nlaggr(1:me+1)) filter_mat = (parms%aggr_filter == mld_filter_mat_) - ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 - call psb_halo(ilaggr,desc_a,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') - goto 9999 - end if ! naggr: number of local aggregates ! nrow: local rows. @@ -172,17 +166,10 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest end if ! 1. Allocate Ptilde in sparse matrix form - call tmpcoo%allocate(ncol,ntaggr,ncol) - do i=1,ncol - tmpcoo%val(i) = done - tmpcoo%ia(i) = i - tmpcoo%ja(i) = ilaggr(i) - end do - call tmpcoo%set_nzeros(ncol) - call tmpcoo%set_dupl(psb_dupl_add_) - call tmpcoo%set_sorted() ! At this point this is in row-major + call op_prol%mv_to(tmpcoo) call ptilde%mv_from_coo(tmpcoo,info) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) + if (info /= psb_success_) goto 9999 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/impl/mld_dcoarse_bld.f90 b/mlprec/impl/mld_dcoarse_bld.f90 index 11b368c6..b439d393 100644 --- a/mlprec/impl/mld_dcoarse_bld.f90 +++ b/mlprec/impl/mld_dcoarse_bld.f90 @@ -98,187 +98,213 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) ictxt = desc_a%get_context() call psb_info(ictxt,me,np) - call mld_check_def(p%parms%ml_type,'Multilevel type',& - & mld_mult_ml_,is_legal_ml_type) - call mld_check_def(p%parms%aggr_alg,'Aggregation',& - & mld_dec_aggr_,is_legal_ml_aggr_alg) - call mld_check_def(p%parms%aggr_ord,'Ordering',& - & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) - call mld_check_def(p%parms%aggr_kind,'Smoother',& - & mld_smooth_prol_,is_legal_ml_aggr_kind) - call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& - & mld_distr_mat_,is_legal_ml_coarse_mat) - call mld_check_def(p%parms%aggr_filter,'Use filtered matrix',& - & mld_no_filter_mat_,is_legal_aggr_filter) - call mld_check_def(p%parms%smoother_pos,'smooth_pos',& - & mld_pre_smooth_,is_legal_ml_smooth_pos) - call mld_check_def(p%parms%aggr_omega_alg,'Omega Alg.',& - & mld_eig_est_,is_legal_ml_aggr_omega_alg) - call mld_check_def(p%parms%aggr_eig,'Eigenvalue estimate',& - & mld_max_norm_,is_legal_ml_aggr_eig) - call mld_check_def(p%parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega) - call mld_check_def(p%parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) - - select case(p%parms%aggr_alg) - case (mld_dec_aggr_, mld_sym_dec_aggr_) - - ! - ! Build a mapping between the row indices of the fine-level matrix - ! and the row indices of the coarse-level matrix, according to a decoupled - ! aggregation algorithm. This also defines a tentative prolongator from - ! the coarse to the fine level. - ! - call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& - & a,desc_a,ilaggr,nlaggr,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld') + + if (.true.) then + call mld_aggrmap_bld(p,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) + else + call mld_check_def(p%parms%ml_type,'Multilevel type',& + & mld_mult_ml_,is_legal_ml_type) + call mld_check_def(p%parms%aggr_alg,'Aggregation',& + & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) + call mld_check_def(p%parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs) + + select case(p%parms%aggr_alg) + case (mld_dec_aggr_, mld_sym_dec_aggr_) + + ! + ! Build a mapping between the row indices of the fine-level matrix + ! and the row indices of the coarse-level matrix, according to a decoupled + ! aggregation algorithm. This also defines a tentative prolongator from + ! the coarse to the fine level. + ! + call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld') + goto 9999 + end if + + case (mld_bcmatch_aggr_) + write(0,*) 'Matching is not implemented yet ' + info = -1111 + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/)) goto 9999 - end if - + + case default + + info = -1 + call psb_errpush(psb_err_input_value_invalid_i_,name,& + & i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/)) + goto 9999 + + end select + end if + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmap_bld') + goto 9999 + end if + + + if (.true.) then + call mld_aggrmat_asb(p,& + & a,desc_a,ilaggr,nlaggr,op_prol,info) + else + + + call mld_check_def(p%parms%aggr_kind,'Smoother',& + & mld_smooth_prol_,is_legal_ml_aggr_kind) + call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& + & mld_distr_mat_,is_legal_ml_coarse_mat) + call mld_check_def(p%parms%aggr_filter,'Use filtered matrix',& + & mld_no_filter_mat_,is_legal_aggr_filter) + call mld_check_def(p%parms%smoother_pos,'smooth_pos',& + & mld_pre_smooth_,is_legal_ml_smooth_pos) + call mld_check_def(p%parms%aggr_omega_alg,'Omega Alg.',& + & mld_eig_est_,is_legal_ml_aggr_omega_alg) + call mld_check_def(p%parms%aggr_eig,'Eigenvalue estimate',& + & mld_max_norm_,is_legal_ml_aggr_eig) + call mld_check_def(p%parms%aggr_omega_val,'Omega',dzero,is_legal_d_omega) + + ! ! Build the coarse-level matrix from the fine-level one, starting from ! the mapping defined by mld_aggrmap_bld and applying the aggregation ! algorithm specified by p%iprcparm(mld_aggr_kind_) ! call mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p%parms,ac,op_prol,op_restr,info) - + if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') goto 9999 end if - case (mld_bcmatch_aggr_) - write(0,*) 'Matching is not implemented yet ' - info = -1111 - call psb_errpush(psb_err_input_value_invalid_i_,name,& - & i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/)) - goto 9999 - - case default - info = -1 - call psb_errpush(psb_err_input_value_invalid_i_,name,& - & i_err=(/ione,p%parms%aggr_alg,izero,izero,izero/)) - goto 9999 + ! Common code refactored here. - end select - + ntaggr = sum(nlaggr) - ! Common code refactored here. - - ntaggr = sum(nlaggr) + select case(p%parms%coarse_mat) - select case(p%parms%coarse_mat) + case(mld_distr_mat_) - case(mld_distr_mat_) - - call ac%mv_to(bcoo) - if (p%parms%clean_zeros) call bcoo%clean_zeros(info) - nzl = bcoo%get_nzeros() - - if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) - if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) - if (info == psb_success_) call psb_cdasb(p%desc_ac,info) - if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') - if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Creating p%desc_ac and converting ac') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Assembld aux descr. distr.' - call p%ac%mv_from(bcoo) + call ac%mv_to(bcoo) + if (p%parms%clean_zeros) call bcoo%clean_zeros(info) + nzl = bcoo%get_nzeros() - call p%ac%set_nrows(p%desc_ac%get_local_rows()) - call p%ac%set_ncols(p%desc_ac%get_local_cols()) - call p%ac%set_asb() + if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) + if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') + if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') + goto 9999 + end if + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Assembld aux descr. distr.' + call p%ac%mv_from(bcoo) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') - goto 9999 - end if + call p%ac%set_nrows(p%desc_ac%get_local_rows()) + call p%ac%set_ncols(p%desc_ac%get_local_cols()) + call p%ac%set_asb() - if (np>1) then - call op_prol%mv_to(acsr1) - nzl = acsr1%get_nzeros() - call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') goto 9999 end if - call op_prol%mv_from(acsr1) - endif - call op_prol%set_ncols(p%desc_ac%get_local_cols()) - - if (np>1) then - call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) - call op_restr%mv_to(acoo) - nzl = acoo%get_nzeros() - if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') - call acoo%set_dupl(psb_dupl_add_) - if (info == psb_success_) call op_restr%mv_from(acoo) - if (info == psb_success_) call op_restr%cscnv(info,type='csr') - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Converting op_restr to local') - goto 9999 + + if (np>1) then + call op_prol%mv_to(acsr1) + nzl = acsr1%get_nzeros() + call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') + goto 9999 + end if + call op_prol%mv_from(acsr1) + endif + call op_prol%set_ncols(p%desc_ac%get_local_cols()) + + if (np>1) then + call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_) + call op_restr%mv_to(acoo) + nzl = acoo%get_nzeros() + if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') + call acoo%set_dupl(psb_dupl_add_) + if (info == psb_success_) call op_restr%mv_from(acoo) + if (info == psb_success_) call op_restr%cscnv(info,type='csr') + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Converting op_restr to local') + goto 9999 + end if end if - end if - ! - ! Clip to local rows. - ! - call op_restr%set_nrows(p%desc_ac%get_local_rows()) + ! + ! Clip to local rows. + ! + call op_restr%set_nrows(p%desc_ac%get_local_rows()) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done ac ' + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Done ac ' - case(mld_repl_mat_) - ! - ! - call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) - if (info == psb_success_) call psb_cdasb(p%desc_ac,info) - if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) - if (info == psb_success_) & - & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) + case(mld_repl_mat_) + ! + ! + call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) + if (info == psb_success_) call psb_cdasb(p%desc_ac,info) + if ((info == psb_success_).and.p%parms%clean_zeros) call ac%clean_zeros(info) + if (info == psb_success_) & + & call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.) - if (info /= psb_success_) goto 9999 + if (info /= psb_success_) goto 9999 - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') - goto 9999 - end select + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') + goto 9999 + end select - call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') - goto 9999 - end if + call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') + goto 9999 + end if - ! - ! Copy the prolongation/restriction matrices into the descriptor map. - ! op_restr => PR^T i.e. restriction operator - ! op_prol => PR i.e. prolongation operator - ! + ! + ! Copy the prolongation/restriction matrices into the descriptor map. + ! op_restr => PR^T i.e. restriction operator + ! op_prol => PR i.e. prolongation operator + ! + + p%map = psb_linmap(psb_map_aggr_,desc_a,& + & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) + if (info == psb_success_) call op_prol%free() + if (info == psb_success_) call op_restr%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + goto 9999 + end if + ! + ! Fix the base_a and base_desc pointers for handling of residuals. + ! This is correct because this routine is only called at levels >=2. + ! + p%base_a => p%ac + p%base_desc => p%desc_ac + end if - p%map = psb_linmap(psb_map_aggr_,desc_a,& - & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr) - if (info == psb_success_) call op_prol%free() - if (info == psb_success_) call op_restr%free() if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') goto 9999 end if - ! - ! Fix the base_a and base_desc pointers for handling of residuals. - ! This is correct because this routine is only called at levels >=2. - ! - p%base_a => p%ac - p%base_desc => p%desc_ac + call psb_erractionrestore(err_act) return diff --git a/mlprec/impl/mld_dprecinit.F90 b/mlprec/impl/mld_dprecinit.F90 index 3276c9e7..c421d04d 100644 --- a/mlprec/impl/mld_dprecinit.F90 +++ b/mlprec/impl/mld_dprecinit.F90 @@ -171,7 +171,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev) nlev_ = max(1,nlev) p%n_prec_levs = nlev_ else - nlev_ = 3 + nlev_ = p%max_prec_levs p%n_prec_levs = -ione end if ilev_ = 1 diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index 30b328b6..2effe9c6 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -103,11 +103,12 @@ module mld_base_prec_type integer(psb_ipk_) :: coarse_mat, coarse_solve logical :: clean_zeros=.true. contains - procedure, pass(pm) :: clone => ml_parms_clone - procedure, pass(pm) :: descr => ml_parms_descr - procedure, pass(pm) :: mldescr => ml_parms_mldescr + procedure, pass(pm) :: get_coarse => ml_parms_get_coarse + procedure, pass(pm) :: clone => ml_parms_clone + procedure, pass(pm) :: descr => ml_parms_descr + procedure, pass(pm) :: mldescr => ml_parms_mldescr procedure, pass(pm) :: coarsedescr => ml_parms_coarsedescr - procedure, pass(pm) :: printout => ml_parms_printout + procedure, pass(pm) :: printout => ml_parms_printout end type mld_ml_parms @@ -506,6 +507,14 @@ contains end select end function mld_stringval + subroutine ml_parms_get_coarse(pm,pmin) + implicit none + class(mld_ml_parms), intent(inout) :: pm + class(mld_ml_parms), intent(in) :: pmin + pm%coarse_mat = pmin%coarse_mat + pm%coarse_solve = pmin%coarse_solve + end subroutine ml_parms_get_coarse + subroutine ml_parms_printout(pm,iout) diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 3c50afd7..f8012673 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -124,11 +124,12 @@ module mld_d_inner_mod end interface mld_bld_mlhier_aggsize interface mld_bld_mlhier_array - subroutine mld_d_bld_mlhier_array(nplevs,a,desc_a,precv,info) - use psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type + subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) + use psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type, psb_dpk_ use mld_d_prec_type, only : mld_d_onelev_type implicit none - integer(psb_ipk_), intent(inout) :: nplevs + integer(psb_ipk_), intent(inout) :: nplevs, casize + real(psb_dpk_) :: mnaggratio type(psb_dspmat_type),intent(in), target :: a type(psb_desc_type), intent(inout), target :: desc_a type(mld_d_onelev_type), allocatable, target, intent(inout) :: precv(:) @@ -137,7 +138,18 @@ module mld_d_inner_mod end interface mld_bld_mlhier_array interface mld_aggrmap_bld - subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) + subroutine mld_d_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + use mld_d_prec_type, only : mld_d_onelev_type + implicit none + type(mld_d_onelev_type), intent(inout), target :: p + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_lev_aggrmap_bld + subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none integer(psb_ipk_), intent(in) :: iorder @@ -146,6 +158,7 @@ module mld_d_inner_mod type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(out) :: op_prol integer(psb_ipk_), intent(out) :: info end subroutine mld_daggrmap_bld end interface mld_aggrmap_bld @@ -165,18 +178,19 @@ module mld_d_inner_mod end interface mld_dec_map_bld -!!$ interface mld_aggrmat_asb -!!$ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) -!!$ use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ -!!$ use mld_d_prec_type, only : mld_d_onelev_type -!!$ implicit none -!!$ type(psb_dspmat_type), intent(in) :: a -!!$ type(psb_desc_type), intent(in) :: desc_a -!!$ integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) -!!$ type(mld_d_onelev_type), intent(inout), target :: p -!!$ integer(psb_ipk_), intent(out) :: info -!!$ end subroutine mld_daggrmat_asb -!!$ end interface mld_aggrmat_asb + interface mld_aggrmat_asb + subroutine mld_d_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + use mld_d_prec_type, only : mld_d_onelev_type + implicit none + type(mld_d_onelev_type), intent(inout), target :: p + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_dspmat_type), intent(inout) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_d_lev_aggrmat_asb + end interface mld_aggrmat_asb @@ -189,7 +203,8 @@ module mld_d_inner_mod type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) type(mld_dml_parms), intent(inout) :: parms - type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol + type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info end subroutine mld_daggrmat_var_asb end interface diff --git a/tests/pdegen/mld_d_pde3d.f90 b/tests/pdegen/mld_d_pde3d.f90 index 07020a0c..d6a62892 100644 --- a/tests/pdegen/mld_d_pde3d.f90 +++ b/tests/pdegen/mld_d_pde3d.f90 @@ -186,6 +186,8 @@ program mld_d_pde3d type(precdata) :: prectype type(psb_d_coo_sparse_mat) :: acoo ! other variables + character(len=20) :: dump_prefix + logical :: dump_sol=.false., dump_prec=.false. integer(psb_ipk_) :: info, i character(len=20) :: name,ch_err @@ -214,7 +216,8 @@ program mld_d_pde3d ! ! get parameters ! - call get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps) + call get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps,& + &dump_prec,dump_prefix) ! ! allocate and fill in the coefficient matrix, rhs and initial guess @@ -380,6 +383,10 @@ program mld_d_pde3d write(psb_out_unit,'("Total memory occupation for PREC: ",i12)') precsize end if + if (dump_prec) call prec%dump(info,prefix=trim(dump_prefix),& + & ac=.true.,solver=.true.,smoother=.true.,rp=.true.,global_num=.true.) + + ! ! cleanup storage and exit ! @@ -404,13 +411,17 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps) + subroutine get_parms(ictxt,kmethd,prectype,afmt,idim,istopc,itmax,itrace,irst,eps,& + & dump_prec,dump_prefix) + integer(psb_ipk_) :: ictxt type(precdata) :: prectype character(len=*) :: kmethd, afmt integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst integer(psb_ipk_) :: np, iam, info real(psb_dpk_) :: eps + logical :: dump_prec + character(len=*) :: dump_prefix character(len=20) :: buffer call psb_info(ictxt, iam, np) @@ -424,6 +435,8 @@ contains call read_data(itrace,psb_inp_unit) call read_data(irst,psb_inp_unit) call read_data(eps,psb_inp_unit) + call read_data(dump_prec,psb_inp_unit) + call read_data(dump_prefix,psb_inp_unit) call read_data(prectype%descr,psb_inp_unit) ! verbose description of the prec call read_data(prectype%prec,psb_inp_unit) ! overall prectype call read_data(prectype%nlevs,psb_inp_unit) ! Prescribed number of levels @@ -462,6 +475,8 @@ contains call psb_bcast(ictxt,itrace) call psb_bcast(ictxt,irst) call psb_bcast(ictxt,eps) + call psb_bcast(ictxt,dump_prec) + call psb_bcast(ictxt,dump_prefix) call psb_bcast(ictxt,prectype%descr) ! verbose description of the prec call psb_bcast(ictxt,prectype%prec) ! overall prectype call psb_bcast(ictxt,prectype%nlevs) ! Prescribed number of levels diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index 66d7129d..f692ddef 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -6,30 +6,32 @@ CSR ! Storage format CSR COO JAD 10 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS +F ! Dump preconditioner on file T F +test-ml-unsm-our ! File prefix for preconditioner dump ML-MUL-RAS-ILU ! Descriptive name for preconditioner (up to 40 chars) ML ! Preconditioner NONE JACOBI BJAC AS ML --1 ! If ML: Prescribed number of levels; if <0, ignore it and use coarse size. --010 ! If ML: Target coarse size. If <0, then use library default +-4 ! If ML: Prescribed number of levels; if <0, ignore it and use coarse size. +-8000 ! If ML: Target coarse size. If <0, then use library default -1.5d0 ! If ML: Minimum aggregation ratio; if <0 use library default -0.10d0 ! If ML: Smoother Aggregation Threshold: >= 0.0 default if <0 -20 ! If ML: Maximum acceptable number of levels; if <0 use library default -SMOOTHED ! Type of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY -SYMDEC ! Type of aggregation: DEC SYMDEC +SMOOTHED ! Type of aggregation: SMOOTHED, UNSMOOTHED, MINENERGY +DEC ! Type of aggregation: DEC SYMDEC NATURAL ! Ordering of aggregation: NATURAL DEGREE MULT ! Type of multilevel correction: ADD MULT KCYCLE VCYCLE WCYCLE KCYCLESYM TWOSIDE ! Side of correction: PRE POST TWOSIDE (ignored for ADD) -4 ! Smoother sweeps +2 ! Smoother sweeps BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML 0 ! Number of overlap layers for AS preconditioner (at finest level) HALO ! AS Restriction operator NONE HALO NONE ! AS Prolongation operator NONE SUM AVG -GS ! Subdomain solver DSCALE ILU MILU ILUT FWGS BWGS MUMPS UMF SLU -4 ! Solver sweeps for GS +FWGS ! Subdomain solver DSCALE ILU MILU ILUT FWGS BWGS MUMPS UMF SLU +1 ! Solver sweeps for GS 0 ! Level-set N for ILU(N), and P for ILUT 1.d-4 ! Threshold T for ILU(T,P) DIST ! Coarse level: matrix distribution DIST REPL BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST MUMPS -ILU ! Coarse level: subsolver DSCALE GS BWGS ILU UMF SLU SLUDIST MUMPS -1 ! Coarse level: Level-set N for ILU(N) +FWGS ! Coarse level: subsolver DSCALE GS BWGS ILU UMF SLU SLUDIST MUMPS +0 ! Coarse level: Level-set N for ILU(N) 1.d-4 ! Coarse level: Threshold T for ILU(T,P) -4 ! Coarse level: Number of Jacobi sweeps +2 ! Coarse level: Number of Jacobi sweeps