From 74ab1d540f2a62339078e52db12fa00c7334e313 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 30 Sep 2016 19:19:23 +0000 Subject: [PATCH] mld2p4-extaggr: mlprec/impl/Makefile mlprec/impl/mld_c_dec_map_bld.f90 mlprec/impl/mld_c_hierarchy_bld.f90 mlprec/impl/mld_c_lev_aggrmap_bld.f90 mlprec/impl/mld_c_lev_aggrmat_asb.f90 mlprec/impl/mld_caggrmap_bld.f90 mlprec/impl/mld_caggrmat_asb.f90 mlprec/impl/mld_caggrmat_biz_asb.f90 mlprec/impl/mld_caggrmat_minnrg_asb.f90 mlprec/impl/mld_caggrmat_nosmth_asb.f90 mlprec/impl/mld_caggrmat_smth_asb.f90 mlprec/impl/mld_d_bld_mlhier_array.f90 mlprec/impl/mld_d_dec_map_bld.f90.new mlprec/impl/mld_d_dec_map_bld.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_smth_asb.f90 mlprec/impl/mld_dcoarse_bld.f90 mlprec/impl/mld_dprecinit.F90 mlprec/impl/mld_s_dec_map_bld.f90 mlprec/impl/mld_s_hierarchy_bld.f90 mlprec/impl/mld_s_lev_aggrmap_bld.f90 mlprec/impl/mld_s_lev_aggrmat_asb.f90 mlprec/impl/mld_saggrmap_bld.f90 mlprec/impl/mld_saggrmat_asb.f90 mlprec/impl/mld_saggrmat_biz_asb.f90 mlprec/impl/mld_saggrmat_minnrg_asb.f90 mlprec/impl/mld_saggrmat_nosmth_asb.f90 mlprec/impl/mld_saggrmat_smth_asb.f90 mlprec/impl/mld_z_dec_map_bld.f90 mlprec/impl/mld_z_hierarchy_bld.f90 mlprec/impl/mld_z_lev_aggrmap_bld.f90 mlprec/impl/mld_z_lev_aggrmat_asb.f90 mlprec/impl/mld_zaggrmap_bld.f90 mlprec/impl/mld_zaggrmat_asb.f90 mlprec/impl/mld_zaggrmat_biz_asb.f90 mlprec/impl/mld_zaggrmat_minnrg_asb.f90 mlprec/impl/mld_zaggrmat_nosmth_asb.f90 mlprec/impl/mld_zaggrmat_smth_asb.f90 mlprec/mld_c_inner_mod.f90 mlprec/mld_d_inner_mod.f90 mlprec/mld_s_inner_mod.f90 mlprec/mld_z_inner_mod.f90 New organization of aggregation routines. --- mlprec/impl/Makefile | 6 +- mlprec/impl/mld_c_dec_map_bld.f90 | 39 ++- mlprec/impl/mld_c_hierarchy_bld.f90 | 228 +++++++++++++++- mlprec/impl/mld_c_lev_aggrmap_bld.f90 | 148 ++++++++++ mlprec/impl/mld_c_lev_aggrmat_asb.f90 | 252 +++++++++++++++++ mlprec/impl/mld_caggrmap_bld.f90 | 29 +- mlprec/impl/mld_caggrmat_asb.f90 | 115 +------- mlprec/impl/mld_caggrmat_biz_asb.f90 | 25 +- mlprec/impl/mld_caggrmat_minnrg_asb.f90 | 29 +- mlprec/impl/mld_caggrmat_nosmth_asb.f90 | 31 +-- mlprec/impl/mld_caggrmat_smth_asb.f90 | 23 +- mlprec/impl/mld_d_bld_mlhier_array.f90 | 54 ++-- mlprec/impl/mld_d_dec_map_bld.f90 | 16 +- mlprec/impl/mld_d_dec_map_bld.f90.new | 343 ++++++++++++++++++++++++ mlprec/impl/mld_daggrmap_bld.f90 | 8 +- mlprec/impl/mld_daggrmat_asb.f90 | 14 +- mlprec/impl/mld_daggrmat_biz_asb.f90 | 2 +- mlprec/impl/mld_daggrmat_minnrg_asb.f90 | 6 +- mlprec/impl/mld_daggrmat_smth_asb.f90 | 2 +- mlprec/impl/mld_dcoarse_bld.f90 | 252 ++++------------- mlprec/impl/mld_dprecinit.F90 | 2 +- mlprec/impl/mld_s_dec_map_bld.f90 | 39 ++- mlprec/impl/mld_s_hierarchy_bld.f90 | 228 +++++++++++++++- mlprec/impl/mld_s_lev_aggrmap_bld.f90 | 148 ++++++++++ mlprec/impl/mld_s_lev_aggrmat_asb.f90 | 252 +++++++++++++++++ mlprec/impl/mld_saggrmap_bld.f90 | 29 +- mlprec/impl/mld_saggrmat_asb.f90 | 115 +------- mlprec/impl/mld_saggrmat_biz_asb.f90 | 25 +- mlprec/impl/mld_saggrmat_minnrg_asb.f90 | 29 +- mlprec/impl/mld_saggrmat_nosmth_asb.f90 | 31 +-- mlprec/impl/mld_saggrmat_smth_asb.f90 | 23 +- mlprec/impl/mld_z_dec_map_bld.f90 | 37 ++- mlprec/impl/mld_z_hierarchy_bld.f90 | 228 +++++++++++++++- mlprec/impl/mld_z_lev_aggrmap_bld.f90 | 148 ++++++++++ mlprec/impl/mld_z_lev_aggrmat_asb.f90 | 252 +++++++++++++++++ mlprec/impl/mld_zaggrmap_bld.f90 | 29 +- mlprec/impl/mld_zaggrmat_asb.f90 | 115 +------- mlprec/impl/mld_zaggrmat_biz_asb.f90 | 25 +- mlprec/impl/mld_zaggrmat_minnrg_asb.f90 | 29 +- mlprec/impl/mld_zaggrmat_nosmth_asb.f90 | 31 +-- mlprec/impl/mld_zaggrmat_smth_asb.f90 | 23 +- mlprec/mld_c_inner_mod.f90 | 30 ++- mlprec/mld_d_inner_mod.f90 | 30 ++- mlprec/mld_s_inner_mod.f90 | 30 ++- mlprec/mld_z_inner_mod.f90 | 30 ++- 45 files changed, 2578 insertions(+), 1002 deletions(-) create mode 100644 mlprec/impl/mld_c_lev_aggrmap_bld.f90 create mode 100644 mlprec/impl/mld_c_lev_aggrmat_asb.f90 create mode 100644 mlprec/impl/mld_d_dec_map_bld.f90.new create mode 100644 mlprec/impl/mld_s_lev_aggrmap_bld.f90 create mode 100644 mlprec/impl/mld_s_lev_aggrmat_asb.f90 create mode 100644 mlprec/impl/mld_z_lev_aggrmap_bld.f90 create mode 100644 mlprec/impl/mld_z_lev_aggrmat_asb.f90 diff --git a/mlprec/impl/Makefile b/mlprec/impl/Makefile index 3db4d32f..6b7c2bf5 100644 --- a/mlprec/impl/Makefile +++ b/mlprec/impl/Makefile @@ -32,19 +32,19 @@ SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o mld_s_bld_mlhier_aggsize.o mld mld_s_ml_prec_bld.o mld_s_hierarchy_bld.o \ mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \ mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \ - $(SMPFOBJS) mld_s_extprol_bld.o + $(SMPFOBJS) mld_s_extprol_bld.o mld_s_lev_aggrmap_bld.o mld_s_lev_aggrmat_asb.o ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o mld_z_bld_mlhier_aggsize.o mld_z_bld_mlhier_array.o \ mld_z_ml_prec_bld.o mld_z_hierarchy_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \ - $(ZMPFOBJS) mld_z_extprol_bld.o + $(ZMPFOBJS) mld_z_extprol_bld.o mld_z_lev_aggrmap_bld.o mld_z_lev_aggrmat_asb.o CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o mld_c_bld_mlhier_aggsize.o mld_c_bld_mlhier_array.o \ mld_c_ml_prec_bld.o mld_c_hierarchy_bld.o \ mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \ mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \ - $(CMPFOBJS) mld_c_extprol_bld.o + $(CMPFOBJS) mld_c_extprol_bld.o mld_c_lev_aggrmap_bld.o mld_c_lev_aggrmat_asb.o INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS) diff --git a/mlprec/impl/mld_c_dec_map_bld.f90 b/mlprec/impl/mld_c_dec_map_bld.f90 index 35452b80..0084692d 100644 --- a/mlprec/impl/mld_c_dec_map_bld.f90 +++ b/mlprec/impl/mld_c_dec_map_bld.f90 @@ -54,12 +54,12 @@ subroutine mld_c_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(:) complex(psb_spk_), allocatable :: val(:), diag(:) integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_c_csr_sparse_mat) :: acsr real(psb_spk_) :: cpling, tcl - logical :: recovery, candidate + logical :: disjoint integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -138,20 +138,15 @@ subroutine mld_c_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 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 +156,7 @@ subroutine mld_c_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 +164,8 @@ subroutine mld_c_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) @@ -189,7 +186,7 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then ip = k cpling = abs(val(k)) end if @@ -208,7 +205,7 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) step3: do ii=1,nr i = idxs(ii) - if (ilaggr(i) == -(nr+1)) then + if (ilaggr(i) < 0) then call a%csget(i,i,nz,irow,icol,val,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -216,10 +213,10 @@ subroutine mld_c_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 = szero + cpling = dzero ip = 0 do k=1, nz j = icol(k) diff --git a/mlprec/impl/mld_c_hierarchy_bld.f90 b/mlprec/impl/mld_c_hierarchy_bld.f90 index 46977d3c..04d8548b 100644 --- a/mlprec/impl/mld_c_hierarchy_bld.f90 +++ b/mlprec/impl/mld_c_hierarchy_bld.f90 @@ -94,9 +94,13 @@ subroutine mld_c_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_spk_) :: mnaggratio - integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize + real(psb_spk_) :: mnaggratio, sizeratio + class(mld_c_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_sml_parms) :: baseparms, medparms, coarseparms + integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_cspmat_type) :: op_prol + type(mld_c_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_c_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,200 @@ subroutine mld_c_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 + 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 + 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 ! - ! Oldstyle with fixed number of levels. + 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 + + ! - nplevs = max(itwo,min(nplevs,mxplevs)) - call mld_bld_mlhier_array(nplevs,a,desc_a,p%precv,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_lev_mat_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 - + + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -220,7 +422,7 @@ subroutine mld_c_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_c_lev_aggrmap_bld.f90 b/mlprec/impl/mld_c_lev_aggrmap_bld.f90 new file mode 100644 index 00000000..0e2034bf --- /dev/null +++ b/mlprec/impl/mld_c_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_ccoarse_bld.f90 +! +! Subroutine: mld_ccoarse_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_cspmat_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_c_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_c_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_lev_aggrmap_bld + + implicit none + + ! Arguments + type(mld_c_onelev_type), intent(inout), target :: p + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_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_c_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',szero,is_legal_s_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_c_lev_aggrmap_bld diff --git a/mlprec/impl/mld_c_lev_aggrmat_asb.f90 b/mlprec/impl/mld_c_lev_aggrmat_asb.f90 new file mode 100644 index 00000000..8dbb2370 --- /dev/null +++ b/mlprec/impl/mld_c_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_ccoarse_bld.f90 +! +! Subroutine: mld_ccoarse_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_cspmat_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_c_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_c_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_lev_aggrmat_asb + + implicit none + + ! Arguments + type(mld_c_onelev_type), intent(inout), target :: p + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_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_cspmat_type) :: ac, op_restr + type(psb_c_coo_sparse_mat) :: acoo, bcoo + type(psb_c_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl, ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_ccoarse_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',szero,is_legal_s_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_caggrmat_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_c_lev_aggrmat_asb diff --git a/mlprec/impl/mld_caggrmap_bld.f90 b/mlprec/impl/mld_caggrmap_bld.f90 index e2ed909b..365c6a87 100644 --- a/mlprec/impl/mld_caggrmap_bld.f90 +++ b/mlprec/impl/mld_caggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmap_bld @@ -93,13 +93,14 @@ subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(out) :: op_prol 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_cspmat_type) :: atmp, atrans - logical :: recovery + type(psb_c_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,28 @@ subroutine mld_caggrmap_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) = cone + 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_caggrmat_asb.f90 b/mlprec/impl/mld_caggrmat_asb.f90 index 82d654e0..1cef24ee 100644 --- a/mlprec/impl/mld_caggrmat_asb.f90 +++ b/mlprec/impl/mld_caggrmat_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmat_asb @@ -109,11 +109,12 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(mld_sml_parms), intent(inout) :: parms type(mld_c_onelev_type), intent(inout), target :: p + type(psb_cspmat_type), intent(inout) :: ac, op_prol,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables - type(psb_cspmat_type) :: ac, op_prol,op_restr type(psb_c_coo_sparse_mat) :: acoo, bcoo type(psb_c_csr_sparse_mat) :: acsr1 integer(psb_ipk_) :: nzl,ntaggr, err_act @@ -165,116 +166,6 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) 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_caggrmat_biz_asb.f90 b/mlprec/impl/mld_caggrmat_biz_asb.f90 index 307ca827..794bb9f1 100644 --- a/mlprec/impl/mld_caggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_caggrmat_biz_asb.f90 @@ -89,7 +89,8 @@ subroutine mld_caggrmat_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_sml_parms), intent(inout) :: parms - type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_cspmat_type), intent(inout) :: op_prol + type(psb_cspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -128,19 +129,8 @@ subroutine mld_caggrmat_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_caggrmat_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) = cone - 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_caggrmat_minnrg_asb.f90 b/mlprec/impl/mld_caggrmat_minnrg_asb.f90 index 7ebcb689..4884252a 100644 --- a/mlprec/impl/mld_caggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_caggrmat_minnrg_asb.f90 @@ -108,8 +108,9 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_sml_parms), intent(inout) :: parms - type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr + type(mld_sml_parms), intent(inout) :: parms + type(psb_cspmat_type), intent(inout) :: op_prol + type(psb_cspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -167,14 +168,6 @@ subroutine mld_caggrmat_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_caggrmat_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) = cone - 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 @@ -276,7 +259,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re call am3%mv_to(acsr3) ! Compute omega_int - ommx = cmplx(szero,szero) + ommx = czero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) @@ -454,7 +437,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = cmplx(szero,szero) + ommx = czero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) diff --git a/mlprec/impl/mld_caggrmat_nosmth_asb.f90 b/mlprec/impl/mld_caggrmat_nosmth_asb.f90 index 91b01614..b47c1dd7 100644 --- a/mlprec/impl/mld_caggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_caggrmat_nosmth_asb.f90 @@ -92,7 +92,8 @@ subroutine mld_caggrmat_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_sml_parms), intent(inout) :: parms - type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_cspmat_type), intent(inout) :: op_prol + type(psb_cspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -124,34 +125,12 @@ subroutine mld_caggrmat_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) = cone - 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_caggrmat_smth_asb.f90 b/mlprec/impl/mld_caggrmat_smth_asb.f90 index 0f24e6d6..83ad51ad 100644 --- a/mlprec/impl/mld_caggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_caggrmat_smth_asb.f90 @@ -104,7 +104,8 @@ subroutine mld_caggrmat_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_sml_parms), intent(inout) :: parms - type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_cspmat_type), intent(inout) :: op_prol + type(psb_cspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -147,14 +148,7 @@ subroutine mld_caggrmat_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_caggrmat_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) = cone - 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_d_bld_mlhier_array.f90 b/mlprec/impl/mld_d_bld_mlhier_array.f90 index bd0a0e77..3618de1b 100644 --- a/mlprec/impl/mld_d_bld_mlhier_array.f90 +++ b/mlprec/impl/mld_d_bld_mlhier_array.f90 @@ -37,25 +37,23 @@ !!$ !!$ -subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) +subroutine mld_d_bld_mlhier_array(nplevs,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, casize - real(psb_dpk_) :: mnaggratio + integer(psb_ipk_), intent(inout) :: nplevs 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, iaggsize + integer(psb_ipk_) :: err,i,k, err_act, newsz, iszv 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 @@ -180,45 +178,27 @@ subroutine mld_d_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) & write(debug_unit,*) me,' ',trim(name),& & 'Return from ',i,' call to mlprcbld ',info - 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 (i>2) then 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 c4128a19..ecbf5a8b 100644 --- a/mlprec/impl/mld_d_dec_map_bld.f90 +++ b/mlprec/impl/mld_d_dec_map_bld.f90 @@ -56,10 +56,10 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& & 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, nisolate,nthird + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_d_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl - logical :: recovery, disjoint + logical :: disjoint integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -114,7 +114,6 @@ 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) @@ -140,11 +139,6 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) end if enddo -!!$ 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. @@ -193,7 +187,6 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) if ((1<=j).and.(j<=nr)) 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 @@ -209,12 +202,10 @@ 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) < 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_ @@ -288,8 +279,7 @@ 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_dec_map_bld.f90.new b/mlprec/impl/mld_d_dec_map_bld.f90.new new file mode 100644 index 00000000..67e4437b --- /dev/null +++ b/mlprec/impl/mld_d_dec_map_bld.f90.new @@ -0,0 +1,343 @@ +!!$ +!!$ +!!$ 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. +!!$ +!!$ + +subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_dec_map_bld + + implicit none + + ! Arguments + integer(psb_ipk_), intent(in) :: iorder + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer(psb_ipk_), intent(out) :: info + + ! Local variables + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:) + real(psb_dpk_), allocatable :: val(:), diag(:) + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + type(psb_d_csr_sparse_mat) :: acsr + real(psb_dpk_) :: cpling, tcl + logical :: recovery + integer(psb_ipk_) :: debug_level, debug_unit,err_act + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if (psb_get_errstatus() /= 0) return + info=psb_success_ + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ictxt=desc_a%get_context() + call psb_info(ictxt,me,np) + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + + nr = a%get_nrows() + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + diag = a%get_diag(info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + call a%cp_to(acsr) + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = acsr%irp(i+1) - acsr%irp(i) + end do + call acsr%free() + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if + ! Note: -(nr+1) Untouched as yet + ! -i 1<=i<=nr Adjacent to aggregate i + ! i 1<=i<=nr Belonging to aggregate i + ! + ! Phase one: group nodes together. + ! Very simple minded strategy. + ! + naggr = 0 + nlp = 0 + do + icnt = 0 + do ii=1, nr + i = idxs(ii) + if (ilaggr(i) == -(nr+1)) then + ! + ! 1. Untouched nodes are marked >0 together + ! with their neighbours + ! + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='csget') + goto 9999 + end if + + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + ilg = ilaggr(j) + if ((ilg<0).and.(i /= j)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ilaggr(j) = naggr + else + ilaggr(j) = -naggr + endif + end if + end if + enddo + + ! + ! 2. Untouched neighbours of these nodes are marked <0. + ! + call a%get_neigh(i,neigh,n_ne,info,lev=2_psb_ipk_) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_neigh') + goto 9999 + end if + + do n = 1, n_ne + m = neigh(n) + if ((1<=m).and.(m<=nr)) then + if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr + endif + enddo + endif + enddo + nlp = nlp + 1 + if (icnt == 0) exit + enddo + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: sweep over leftovers. + ! + allocate(ils(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/naggr+10,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + do i=1, size(ils) + ils(i) = 0 + end do + do i=1, nr + n = ilaggr(i) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') + goto 9999 + else + ils(n) = ils(n) + 1 + end if + + end if + end do + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: number of aggregates ',naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: nodes aggregated ',sum(ils) + end if + + recovery=.false. + do i=1, nr + if (ilaggr(i) < 0) then + ! + ! Now some silly rule to break ties: + ! Group with adjacent aggregate. + ! + isz = nr+1 + ia = -1 + cpling = dzero + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getrow') + goto 9999 + end if + + do j=1, nz + k = icol(j) + if ((1<=k).and.(k<=nr).and.(k /= i)) then + tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k))) + if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then +!!$ if (tcl > theta) then + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') + goto 9999 + end if + + if ((abs(val(j))>cpling) .or. & + & ((abs(val(j)) == cpling).and. (ils(n) < isz))) then +!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then + ia = n + isz = ils(n) + cpling = abs(val(j)) +!!$ cpling = tcl + endif + endif + endif + end if + enddo + + if (ia == -1) then + ! At this point, the easiest thing is to start a new aggregate + naggr = naggr + 1 + ilaggr(i) = naggr + ils(ilaggr(i)) = 1 + + else + + ilaggr(i) = ia + + if (ia>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') + goto 9999 + end if + ils(ia) = ils(ia) + 1 + endif + + end if + end do + if (debug_level >= psb_debug_outer_) then + if (recovery) then + write(debug_unit,*) me,' ',trim(name),& + & 'Had to recover from strange situation in loc_aggregate.' + write(debug_unit,*) me,' ',trim(name),& + & 'Perhaps an unsymmetric pattern?' + endif + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 2: number of aggregates ',naggr,sum(ils) + do i=1, naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) + enddo + write(debug_unit,*) me,' ',trim(name),& + & maxval(ils(1:naggr)) + write(debug_unit,*) me,' ',trim(name),& + & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' + end if + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + deallocate(ils,neigh,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,izero,izero,izero,izero/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_d_dec_map_bld + diff --git a/mlprec/impl/mld_daggrmap_bld.f90 b/mlprec/impl/mld_daggrmap_bld.f90 index de053aa6..f8e6d7cc 100644 --- a/mlprec/impl/mld_daggrmap_bld.f90 +++ b/mlprec/impl/mld_daggrmap_bld.f90 @@ -92,8 +92,8 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_pro 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(:) + type(psb_dspmat_type), intent(out) :: op_prol integer(psb_ipk_), intent(out) :: info ! Local variables @@ -152,14 +152,12 @@ subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_pro goto 9999 end if - naggr = nlaggr(me+1) - ntaggr = sum(nlaggr) - + 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 diff --git a/mlprec/impl/mld_daggrmat_asb.f90 b/mlprec/impl/mld_daggrmat_asb.f90 index dd8dd9e4..2c7ab573 100644 --- a/mlprec/impl/mld_daggrmat_asb.f90 +++ b/mlprec/impl/mld_daggrmat_asb.f90 @@ -110,9 +110,9 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,inf type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) type(mld_dml_parms), intent(inout) :: parms -!!$ type(mld_d_onelev_type), intent(inout), target :: p - integer(psb_ipk_), intent(out) :: info + type(mld_d_onelev_type), intent(inout), target :: p type(psb_dspmat_type), intent(inout) :: ac, op_prol,op_restr + integer(psb_ipk_), intent(out) :: info ! Local variables type(psb_d_coo_sparse_mat) :: acoo, bcoo @@ -134,26 +134,26 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,inf call psb_info(ictxt, me, np) - select case (parms%aggr_kind) + select case (p%parms%aggr_kind) case (mld_no_smooth_) call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,& - & parms,ac,op_prol,op_restr,info) + & p%parms,ac,op_prol,op_restr,info) case(mld_smooth_prol_) call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, & - & parms,ac,op_prol,op_restr,info) + & p%parms,ac,op_prol,op_restr,info) case(mld_biz_prol_) call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, & - & parms,ac,op_prol,op_restr,info) + & p%parms,ac,op_prol,op_restr,info) case(mld_min_energy_) call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, & - & parms,ac,op_prol,op_restr,info) + & p%parms,ac,op_prol,op_restr,info) case default info = psb_err_internal_error_ diff --git a/mlprec/impl/mld_daggrmat_biz_asb.f90 b/mlprec/impl/mld_daggrmat_biz_asb.f90 index 9d8dfa79..3e7d2a8a 100644 --- a/mlprec/impl/mld_daggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_daggrmat_biz_asb.f90 @@ -89,7 +89,7 @@ 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(inout) :: op_prol + type(psb_dspmat_type), intent(inout) :: op_prol type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 index a08ab607..965443ca 100644 --- a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 @@ -108,7 +108,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_dml_parms), intent(inout) :: parms + type(mld_dml_parms), intent(inout) :: parms type(psb_dspmat_type), intent(inout) :: op_prol type(psb_dspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info @@ -259,7 +259,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re call am3%mv_to(acsr3) ! Compute omega_int - ommx = cmplx(dzero,dzero) + ommx = dzero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) @@ -437,7 +437,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = cmplx(dzero,dzero) + ommx = dzero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) diff --git a/mlprec/impl/mld_daggrmat_smth_asb.f90 b/mlprec/impl/mld_daggrmat_smth_asb.f90 index 46b34133..455a7f2a 100644 --- a/mlprec/impl/mld_daggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_daggrmat_smth_asb.f90 @@ -148,7 +148,7 @@ 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_) - + ! ! naggr: number of local aggregates ! nrow: local rows. ! diff --git a/mlprec/impl/mld_dcoarse_bld.f90 b/mlprec/impl/mld_dcoarse_bld.f90 index b439d393..b52c4215 100644 --- a/mlprec/impl/mld_dcoarse_bld.f90 +++ b/mlprec/impl/mld_dcoarse_bld.f90 @@ -83,229 +83,69 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) integer(psb_mpik_) :: ictxt, np, me integer(psb_ipk_) :: err_act integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) - type(psb_dspmat_type) :: ac, op_prol,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_ + info = psb_success_ ictxt = desc_a%get_context() call psb_info(ictxt,me,np) - - 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 - - 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 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) + + + ! + ! 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') 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 - - - ! 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 - end if + ! + ! 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_aggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) if(info /= psb_success_) then 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 c421d04d..3276c9e7 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_ = p%max_prec_levs + nlev_ = 3 p%n_prec_levs = -ione end if ilev_ = 1 diff --git a/mlprec/impl/mld_s_dec_map_bld.f90 b/mlprec/impl/mld_s_dec_map_bld.f90 index 10be26e0..afd107fa 100644 --- a/mlprec/impl/mld_s_dec_map_bld.f90 +++ b/mlprec/impl/mld_s_dec_map_bld.f90 @@ -54,12 +54,12 @@ subroutine mld_s_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_spk_), allocatable :: val(:), diag(:) integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_s_csr_sparse_mat) :: acsr real(psb_spk_) :: cpling, tcl - logical :: recovery, candidate + logical :: disjoint integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -138,20 +138,15 @@ subroutine mld_s_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 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 +156,7 @@ subroutine mld_s_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 +164,8 @@ subroutine mld_s_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) @@ -189,7 +186,7 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then ip = k cpling = abs(val(k)) end if @@ -208,7 +205,7 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) step3: do ii=1,nr i = idxs(ii) - if (ilaggr(i) == -(nr+1)) then + if (ilaggr(i) < 0) then call a%csget(i,i,nz,irow,icol,val,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -216,10 +213,10 @@ subroutine mld_s_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 = szero + cpling = dzero ip = 0 do k=1, nz j = icol(k) diff --git a/mlprec/impl/mld_s_hierarchy_bld.f90 b/mlprec/impl/mld_s_hierarchy_bld.f90 index 5216b314..3faf3f57 100644 --- a/mlprec/impl/mld_s_hierarchy_bld.f90 +++ b/mlprec/impl/mld_s_hierarchy_bld.f90 @@ -94,9 +94,13 @@ subroutine mld_s_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_spk_) :: mnaggratio - integer(psb_ipk_) :: ipv(mld_ifpsz_), val + integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs, iaggsize + real(psb_spk_) :: mnaggratio, sizeratio + class(mld_s_base_smoother_type), allocatable :: coarse_sm, base_sm, med_sm + type(mld_sml_parms) :: baseparms, medparms, coarseparms + integer(psb_ipk_), allocatable :: ilaggr(:), nlaggr(:) + type(psb_sspmat_type) :: op_prol + type(mld_s_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_s_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,200 @@ subroutine mld_s_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 + 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 + 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 ! - ! Oldstyle with fixed number of levels. + 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 + + ! - nplevs = max(itwo,min(nplevs,mxplevs)) - call mld_bld_mlhier_array(nplevs,a,desc_a,p%precv,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_lev_mat_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 - + + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -220,7 +422,7 @@ subroutine mld_s_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_s_lev_aggrmap_bld.f90 b/mlprec/impl/mld_s_lev_aggrmap_bld.f90 new file mode 100644 index 00000000..fa4b78aa --- /dev/null +++ b/mlprec/impl/mld_s_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_scoarse_bld.f90 +! +! Subroutine: mld_scoarse_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_sspmat_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_s_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_s_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_lev_aggrmap_bld + + implicit none + + ! Arguments + type(mld_s_onelev_type), intent(inout), target :: p + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_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_s_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',szero,is_legal_s_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_s_lev_aggrmap_bld diff --git a/mlprec/impl/mld_s_lev_aggrmat_asb.f90 b/mlprec/impl/mld_s_lev_aggrmat_asb.f90 new file mode 100644 index 00000000..0d746a94 --- /dev/null +++ b/mlprec/impl/mld_s_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_scoarse_bld.f90 +! +! Subroutine: mld_scoarse_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_sspmat_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_s_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_s_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_lev_aggrmat_asb + + implicit none + + ! Arguments + type(mld_s_onelev_type), intent(inout), target :: p + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_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_sspmat_type) :: ac, op_restr + type(psb_s_coo_sparse_mat) :: acoo, bcoo + type(psb_s_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl, ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_scoarse_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',szero,is_legal_s_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_saggrmat_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_s_lev_aggrmat_asb diff --git a/mlprec/impl/mld_saggrmap_bld.f90 b/mlprec/impl/mld_saggrmap_bld.f90 index b46bb98a..5dfd1af5 100644 --- a/mlprec/impl/mld_saggrmap_bld.f90 +++ b/mlprec/impl/mld_saggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmap_bld @@ -93,13 +93,14 @@ subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(out) :: op_prol 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_sspmat_type) :: atmp, atrans - logical :: recovery + type(psb_s_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,28 @@ subroutine mld_saggrmap_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) = sone + 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_saggrmat_asb.f90 b/mlprec/impl/mld_saggrmat_asb.f90 index c0fe4b22..36bcb4db 100644 --- a/mlprec/impl/mld_saggrmat_asb.f90 +++ b/mlprec/impl/mld_saggrmat_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmat_asb @@ -109,11 +109,12 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(mld_sml_parms), intent(inout) :: parms type(mld_s_onelev_type), intent(inout), target :: p + type(psb_sspmat_type), intent(inout) :: ac, op_prol,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables - type(psb_sspmat_type) :: ac, op_prol,op_restr type(psb_s_coo_sparse_mat) :: acoo, bcoo type(psb_s_csr_sparse_mat) :: acsr1 integer(psb_ipk_) :: nzl,ntaggr, err_act @@ -165,116 +166,6 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) 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_saggrmat_biz_asb.f90 b/mlprec/impl/mld_saggrmat_biz_asb.f90 index 3e715cd6..f91c8102 100644 --- a/mlprec/impl/mld_saggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_saggrmat_biz_asb.f90 @@ -89,7 +89,8 @@ subroutine mld_saggrmat_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_sml_parms), intent(inout) :: parms - type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_sspmat_type), intent(inout) :: op_prol + type(psb_sspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -128,19 +129,8 @@ subroutine mld_saggrmat_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_saggrmat_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) = sone - 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_saggrmat_minnrg_asb.f90 b/mlprec/impl/mld_saggrmat_minnrg_asb.f90 index e9e15e4a..3868d9c6 100644 --- a/mlprec/impl/mld_saggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_saggrmat_minnrg_asb.f90 @@ -108,8 +108,9 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_sml_parms), intent(inout) :: parms - type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr + type(mld_sml_parms), intent(inout) :: parms + type(psb_sspmat_type), intent(inout) :: op_prol + type(psb_sspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -167,14 +168,6 @@ subroutine mld_saggrmat_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_saggrmat_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) = sone - 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 @@ -276,7 +259,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re call am3%mv_to(acsr3) ! Compute omega_int - ommx = cmplx(szero,szero) + ommx = szero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) @@ -454,7 +437,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = cmplx(szero,szero) + ommx = szero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) diff --git a/mlprec/impl/mld_saggrmat_nosmth_asb.f90 b/mlprec/impl/mld_saggrmat_nosmth_asb.f90 index 78608ea0..bbb8246f 100644 --- a/mlprec/impl/mld_saggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_saggrmat_nosmth_asb.f90 @@ -92,7 +92,8 @@ subroutine mld_saggrmat_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_sml_parms), intent(inout) :: parms - type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_sspmat_type), intent(inout) :: op_prol + type(psb_sspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -124,34 +125,12 @@ subroutine mld_saggrmat_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) = sone - 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_saggrmat_smth_asb.f90 b/mlprec/impl/mld_saggrmat_smth_asb.f90 index 1d00dda1..03eb2155 100644 --- a/mlprec/impl/mld_saggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_saggrmat_smth_asb.f90 @@ -104,7 +104,8 @@ subroutine mld_saggrmat_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_sml_parms), intent(inout) :: parms - type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_sspmat_type), intent(inout) :: op_prol + type(psb_sspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -147,14 +148,7 @@ subroutine mld_saggrmat_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_saggrmat_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) = sone - 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_z_dec_map_bld.f90 b/mlprec/impl/mld_z_dec_map_bld.f90 index 876b8a7e..dffb3402 100644 --- a/mlprec/impl/mld_z_dec_map_bld.f90 +++ b/mlprec/impl/mld_z_dec_map_bld.f90 @@ -54,12 +54,12 @@ subroutine mld_z_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(:) complex(psb_dpk_), allocatable :: val(:), diag(:) integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_z_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl - logical :: recovery, candidate + logical :: disjoint integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -138,20 +138,15 @@ subroutine mld_z_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 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 +156,7 @@ subroutine mld_z_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 +164,8 @@ subroutine mld_z_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) @@ -189,7 +186,7 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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 + & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then ip = k cpling = abs(val(k)) end if @@ -208,7 +205,7 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) step3: do ii=1,nr i = idxs(ii) - if (ilaggr(i) == -(nr+1)) then + if (ilaggr(i) < 0) then call a%csget(i,i,nz,irow,icol,val,info) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -216,8 +213,8 @@ subroutine mld_z_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 diff --git a/mlprec/impl/mld_z_hierarchy_bld.f90 b/mlprec/impl/mld_z_hierarchy_bld.f90 index ab14a40a..47689ee5 100644 --- a/mlprec/impl/mld_z_hierarchy_bld.f90 +++ b/mlprec/impl/mld_z_hierarchy_bld.f90 @@ -94,9 +94,13 @@ subroutine mld_z_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_z_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_zspmat_type) :: op_prol + type(mld_z_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_z_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,200 @@ subroutine mld_z_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 + 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 + 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 ! - ! Oldstyle with fixed number of levels. + 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 + + ! - nplevs = max(itwo,min(nplevs,mxplevs)) - call mld_bld_mlhier_array(nplevs,a,desc_a,p%precv,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_lev_mat_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 - + + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Internal hierarchy build' ) @@ -220,7 +422,7 @@ subroutine mld_z_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_z_lev_aggrmap_bld.f90 b/mlprec/impl/mld_z_lev_aggrmap_bld.f90 new file mode 100644 index 00000000..52892fde --- /dev/null +++ b/mlprec/impl/mld_z_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_zcoarse_bld.f90 +! +! Subroutine: mld_zcoarse_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_zspmat_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_z_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_z_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_lev_aggrmap_bld + + implicit none + + ! Arguments + type(mld_z_onelev_type), intent(inout), target :: p + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_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_z_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_z_lev_aggrmap_bld diff --git a/mlprec/impl/mld_z_lev_aggrmat_asb.f90 b/mlprec/impl/mld_z_lev_aggrmat_asb.f90 new file mode 100644 index 00000000..6101fbed --- /dev/null +++ b/mlprec/impl/mld_z_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_zcoarse_bld.f90 +! +! Subroutine: mld_zcoarse_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_zspmat_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_z_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_z_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_lev_aggrmat_asb + + implicit none + + ! Arguments + type(mld_z_onelev_type), intent(inout), target :: p + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_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_zspmat_type) :: ac, op_restr + type(psb_z_coo_sparse_mat) :: acoo, bcoo + type(psb_z_csr_sparse_mat) :: acsr1 + integer(psb_ipk_) :: nzl, ntaggr + integer(psb_ipk_) :: debug_level, debug_unit + + name='mld_zcoarse_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_zaggrmat_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_z_lev_aggrmat_asb diff --git a/mlprec/impl/mld_zaggrmap_bld.f90 b/mlprec/impl/mld_zaggrmap_bld.f90 index 2058fd66..9af55340 100644 --- a/mlprec/impl/mld_zaggrmap_bld.f90 +++ b/mlprec/impl/mld_zaggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmap_bld @@ -93,13 +93,14 @@ subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(out) :: op_prol 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_zspmat_type) :: atmp, atrans - logical :: recovery + type(psb_z_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,28 @@ subroutine mld_zaggrmap_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) = zone + 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_zaggrmat_asb.f90 b/mlprec/impl/mld_zaggrmat_asb.f90 index c6ebe19a..f7a4dbe7 100644 --- a/mlprec/impl/mld_zaggrmat_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_asb.f90 @@ -98,7 +98,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) +subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_asb @@ -109,11 +109,12 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) + type(mld_dml_parms), intent(inout) :: parms type(mld_z_onelev_type), intent(inout), target :: p + type(psb_zspmat_type), intent(inout) :: ac, op_prol,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables - type(psb_zspmat_type) :: ac, op_prol,op_restr type(psb_z_coo_sparse_mat) :: acoo, bcoo type(psb_z_csr_sparse_mat) :: acsr1 integer(psb_ipk_) :: nzl,ntaggr, err_act @@ -165,116 +166,6 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) 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_zaggrmat_biz_asb.f90 b/mlprec/impl/mld_zaggrmat_biz_asb.f90 index 64bcf643..96b3f49f 100644 --- a/mlprec/impl/mld_zaggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_biz_asb.f90 @@ -89,7 +89,8 @@ subroutine mld_zaggrmat_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_zspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_zspmat_type), intent(inout) :: op_prol + type(psb_zspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -128,19 +129,8 @@ subroutine mld_zaggrmat_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_zaggrmat_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) = zone - 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_zaggrmat_minnrg_asb.f90 b/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 index ba0238e0..80143f4d 100644 --- a/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 @@ -108,8 +108,9 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer(psb_ipk_), intent(inout) :: ilaggr(:), nlaggr(:) - type(mld_dml_parms), intent(inout) :: parms - type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr + type(mld_dml_parms), intent(inout) :: parms + type(psb_zspmat_type), intent(inout) :: op_prol + type(psb_zspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -167,14 +168,6 @@ subroutine mld_zaggrmat_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_zaggrmat_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) = zone - 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 @@ -276,7 +259,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re call am3%mv_to(acsr3) ! Compute omega_int - ommx = cmplx(dzero,dzero) + ommx = zzero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) @@ -454,7 +437,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = cmplx(dzero,dzero) + ommx = zzero do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) diff --git a/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 b/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 index 580b4e73..4d3960bb 100644 --- a/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_nosmth_asb.f90 @@ -92,7 +92,8 @@ subroutine mld_zaggrmat_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_zspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_zspmat_type), intent(inout) :: op_prol + type(psb_zspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -124,34 +125,12 @@ subroutine mld_zaggrmat_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) = zone - 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_zaggrmat_smth_asb.f90 b/mlprec/impl/mld_zaggrmat_smth_asb.f90 index 794210b5..8e85e14e 100644 --- a/mlprec/impl/mld_zaggrmat_smth_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_smth_asb.f90 @@ -104,7 +104,8 @@ subroutine mld_zaggrmat_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_zspmat_type), intent(out) :: ac,op_prol,op_restr + type(psb_zspmat_type), intent(inout) :: op_prol + type(psb_zspmat_type), intent(out) :: ac,op_restr integer(psb_ipk_), intent(out) :: info ! Local variables @@ -147,14 +148,7 @@ subroutine mld_zaggrmat_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_zaggrmat_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) = zone - 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/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index 9da1ce64..6b349f97 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -124,11 +124,12 @@ module mld_c_inner_mod end interface mld_bld_mlhier_aggsize interface mld_bld_mlhier_array - subroutine mld_c_bld_mlhier_array(nplevs,a,desc_a,precv,info) + subroutine mld_c_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) use psb_base_mod, only : psb_ipk_, psb_cspmat_type, psb_desc_type use mld_c_prec_type, only : mld_c_onelev_type implicit none - integer(psb_ipk_), intent(inout) :: nplevs + integer(psb_ipk_), intent(inout) :: nplevs, casize + real(psb_spk_) :: mnaggratio type(psb_cspmat_type),intent(in), target :: a type(psb_desc_type), intent(inout), target :: desc_a type(mld_c_onelev_type), allocatable, target, intent(inout) :: precv(:) @@ -137,6 +138,17 @@ module mld_c_inner_mod end interface mld_bld_mlhier_array interface mld_aggrmap_bld + subroutine mld_c_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + use mld_c_prec_type, only : mld_c_onelev_type + implicit none + type(mld_c_onelev_type), intent(inout), target :: p + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_lev_aggrmap_bld subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ implicit none @@ -165,6 +177,20 @@ module mld_c_inner_mod end interface mld_dec_map_bld + interface mld_lev_mat_asb + subroutine mld_c_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + use mld_c_prec_type, only : mld_c_onelev_type + implicit none + type(mld_c_onelev_type), intent(inout), target :: p + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_cspmat_type), intent(inout) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_lev_aggrmat_asb + end interface mld_lev_mat_asb + interface mld_aggrmat_asb subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 27599ffe..309119ce 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -125,11 +125,11 @@ module mld_d_inner_mod interface mld_bld_mlhier_array 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 psb_base_mod, only : psb_ipk_, psb_dspmat_type, psb_desc_type use mld_d_prec_type, only : mld_d_onelev_type implicit none - integer(psb_ipk_), intent(inout) :: nplevs, casize - real(psb_dpk_) :: mnaggratio + 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(:) @@ -145,11 +145,11 @@ module mld_d_inner_mod 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(:) + integer(psb_ipk_), 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) + subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none integer(psb_ipk_), intent(in) :: iorder @@ -158,7 +158,6 @@ 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 @@ -192,6 +191,20 @@ module mld_d_inner_mod end subroutine mld_d_lev_aggrmat_asb end interface mld_lev_mat_asb + 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 + + abstract interface subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) @@ -202,8 +215,7 @@ 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(inout) :: op_prol - type(psb_dspmat_type), intent(out) :: ac,op_restr + type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr integer(psb_ipk_), intent(out) :: info end subroutine mld_daggrmat_var_asb end interface @@ -211,7 +223,7 @@ module mld_d_inner_mod procedure(mld_daggrmat_var_asb) :: mld_daggrmat_nosmth_asb, & & mld_daggrmat_smth_asb, mld_daggrmat_minnrg_asb, & - & mld_daggrmat_biz_asb, mld_daggrmat_asb + & mld_daggrmat_biz_asb end module mld_d_inner_mod diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index ed6dfa77..3adf013d 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -124,11 +124,12 @@ module mld_s_inner_mod end interface mld_bld_mlhier_aggsize interface mld_bld_mlhier_array - subroutine mld_s_bld_mlhier_array(nplevs,a,desc_a,precv,info) + subroutine mld_s_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) use psb_base_mod, only : psb_ipk_, psb_sspmat_type, psb_desc_type use mld_s_prec_type, only : mld_s_onelev_type implicit none - integer(psb_ipk_), intent(inout) :: nplevs + integer(psb_ipk_), intent(inout) :: nplevs, casize + real(psb_spk_) :: mnaggratio type(psb_sspmat_type),intent(in), target :: a type(psb_desc_type), intent(inout), target :: desc_a type(mld_s_onelev_type), allocatable, target, intent(inout) :: precv(:) @@ -137,6 +138,17 @@ module mld_s_inner_mod end interface mld_bld_mlhier_array interface mld_aggrmap_bld + subroutine mld_s_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + use mld_s_prec_type, only : mld_s_onelev_type + implicit none + type(mld_s_onelev_type), intent(inout), target :: p + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_lev_aggrmap_bld subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ implicit none @@ -165,6 +177,20 @@ module mld_s_inner_mod end interface mld_dec_map_bld + interface mld_lev_mat_asb + subroutine mld_s_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ + use mld_s_prec_type, only : mld_s_onelev_type + implicit none + type(mld_s_onelev_type), intent(inout), target :: p + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_sspmat_type), intent(inout) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_lev_aggrmat_asb + end interface mld_lev_mat_asb + interface mld_aggrmat_asb subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index 2af263d8..b1f811ce 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -124,11 +124,12 @@ module mld_z_inner_mod end interface mld_bld_mlhier_aggsize interface mld_bld_mlhier_array - subroutine mld_z_bld_mlhier_array(nplevs,a,desc_a,precv,info) + subroutine mld_z_bld_mlhier_array(nplevs,casize,mnaggratio,a,desc_a,precv,info) use psb_base_mod, only : psb_ipk_, psb_zspmat_type, psb_desc_type use mld_z_prec_type, only : mld_z_onelev_type implicit none - integer(psb_ipk_), intent(inout) :: nplevs + integer(psb_ipk_), intent(inout) :: nplevs, casize + real(psb_dpk_) :: mnaggratio type(psb_zspmat_type),intent(in), target :: a type(psb_desc_type), intent(inout), target :: desc_a type(mld_z_onelev_type), allocatable, target, intent(inout) :: precv(:) @@ -137,6 +138,17 @@ module mld_z_inner_mod end interface mld_bld_mlhier_array interface mld_aggrmap_bld + subroutine mld_z_lev_aggrmap_bld(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + use mld_z_prec_type, only : mld_z_onelev_type + implicit none + type(mld_z_onelev_type), intent(inout), target :: p + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(out) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_lev_aggrmap_bld subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none @@ -165,6 +177,20 @@ module mld_z_inner_mod end interface mld_dec_map_bld + interface mld_lev_mat_asb + subroutine mld_z_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) + use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ + use mld_z_prec_type, only : mld_z_onelev_type + implicit none + type(mld_z_onelev_type), intent(inout), target :: p + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(inout) :: ilaggr(:),nlaggr(:) + type(psb_zspmat_type), intent(inout) :: op_prol + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_lev_aggrmat_asb + end interface mld_lev_mat_asb + interface mld_aggrmat_asb subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_