diff --git a/mlprec/impl/mld_d_dec_map_bld.f90 b/mlprec/impl/mld_d_dec_map_bld.f90 index 363aeec0..c4128a19 100644 --- a/mlprec/impl/mld_d_dec_map_bld.f90 +++ b/mlprec/impl/mld_d_dec_map_bld.f90 @@ -191,9 +191,9 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) do k=1, nz j = icol(k) if ((1<=j).and.(j<=nr)) then -!!$ if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& -!!$ & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then - if ((tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then +!!$ if ((tmpaggr(j) > 0).and. (abs(val(k)) > cpling)) then ip = k cpling = abs(val(k)) end if diff --git a/mlprec/impl/mld_d_hierarchy_bld.f90 b/mlprec/impl/mld_d_hierarchy_bld.f90 index 7bf307bd..42a33e94 100644 --- a/mlprec/impl/mld_d_hierarchy_bld.f90 +++ b/mlprec/impl/mld_d_hierarchy_bld.f90 @@ -287,255 +287,134 @@ subroutine mld_d_hierarchy_bld(a,desc_a,p,info,amold,vmold,imold) p%precv(1)%base_a => a p%precv(1)%base_desc => desc_a newsz = 0 - if (.false.) then - old_array_build_loop: do i=2, iszv + array_build_loop: do i=2, iszv + ! + ! Check on the iprcparm contents: they should be the same + ! on all processes. + ! + call psb_bcast(ictxt,p%precv(i)%parms) + + ! + ! Sanity checks on the parameters + ! + if (i= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Calling mlprcbld at level ',i + ! + ! Build the mapping between levels i-1 and i and the matrix + ! at level i + ! + if (info == psb_success_) call mld_aggrmap_bld(p%precv(i),& + & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& + & ilaggr,nlaggr,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i - ! - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Map build') + goto 9999 + endif - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Init upper level preconditioner') - goto 9999 - endif + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Return from ',i,' call to mlprcbld ',info - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info - iaggsize = sum(p%precv(i)%map%naggr) - if (iaggsize <= casize) then - newsz = i - end if + ! + ! 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 + if (i>2) then + sizeratio = iaggsize + sizeratio = sum(p%precv(i-1)%map%naggr)/sizeratio + if (sizeratio < mnaggratio) then + if (sizeratio > 1) then + newsz = i + else + ! + ! We are not gaining + ! + newsz = newsz-1 end if + end if - if (all(p%precv(i)%map%naggr == p%precv(i-1)%map%naggr)) then - newsz=i-1 - if (me == 0) then - write(debug_unit,*) trim(name),& - &': Warning: aggregates from level ',& - & newsz - write(debug_unit,*) trim(name),& - &': to level ',& - & iszv,' coincide.' - write(debug_unit,*) trim(name),& - &': Number of levels actually used :',newsz - write(debug_unit,*) - end if + 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 - call psb_bcast(ictxt,newsz) - if (newsz > 0) exit old_array_build_loop - end do old_array_build_loop - - if (newsz > 0) then - ! - ! We exited early from the build loop, need to fix - ! the size. - ! - allocate(tprecv(newsz),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,& - & a_err='prec reallocation') - goto 9999 - endif - do i=1,newsz-1 - call p%precv(i)%move_alloc(tprecv(i),info) - end do - call p%precv(iszv)%move_alloc(tprecv(newsz),info) - do i=newsz+1, iszv - call p%precv(i)%free(info) - end do - call move_alloc(tprecv,p%precv) - ! Ignore errors from transfer - info = psb_success_ - ! - ! Restart - iszv = newsz - ! Fix the pointers, but the level 1 should - ! be already OK - do i=2, iszv - 1 - p%precv(i)%base_a => p%precv(i)%ac - p%precv(i)%base_desc => p%precv(i)%desc_ac - p%precv(i)%map%p_desc_X => p%precv(i-1)%base_desc - p%precv(i)%map%p_desc_Y => p%precv(i)%base_desc - end do - - i = iszv - if (info == psb_success_) call mld_coarse_bld(p%precv(i-1)%base_a,& - & p%precv(i-1)%base_desc, p%precv(i),info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='coarse rebuild') - goto 9999 - endif end if + call psb_bcast(ictxt,newsz) + if (newsz > 0) & + & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) - else - - array_build_loop: do i=2, iszv - ! - ! Check on the iprcparm contents: they should be the same - ! on all processes. - ! - call psb_bcast(ictxt,p%precv(i)%parms) - - ! - ! Sanity checks on the parameters - ! - if (i= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Calling mlprcbld at level ',i - ! - ! Build the mapping between levels i-1 and i and the matrix - ! at level i - ! - if (info == psb_success_) call mld_aggrmap_bld(p%precv(i),& - & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Map build') - goto 9999 - endif - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Return from ',i,' call to mlprcbld ',info + 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 - ! - ! Check for early termination of aggregation loop. - ! - iaggsize = sum(nlaggr) - if (iaggsize <= casize) then - newsz = i - end if - - if (i>2) then - sizeratio = iaggsize - sizeratio = sum(p%precv(i-1)%map%naggr)/sizeratio - if (sizeratio < mnaggratio) then - if (sizeratio > 1) then - newsz = i - else - ! - ! We are not gaining - ! - newsz = newsz-1 - end if - end if - - if (all(nlaggr == p%precv(i-1)%map%naggr)) then - newsz=i-1 - if (me == 0) then - write(debug_unit,*) trim(name),& - &': Warning: aggregates from level ',& - & newsz - write(debug_unit,*) trim(name),& - &': to level ',& - & iszv,' coincide.' - write(debug_unit,*) trim(name),& - &': Number of levels actually used :',newsz - write(debug_unit,*) - end if - end if - end if - call psb_bcast(ictxt,newsz) - if (newsz > 0) & - & call p%precv(i)%parms%get_coarse(p%precv(iszv)%parms) - - if (info == psb_success_) call mld_aggrmat_asb(p%precv(i),& - & p%precv(i-1)%base_a,p%precv(i-1)%base_desc,& - & ilaggr,nlaggr,op_prol,info) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Map build') - goto 9999 - endif - - if (newsz > 0) exit array_build_loop - end do array_build_loop + if (newsz > 0) 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 (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' ) diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index f8012673..27599ffe 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -178,7 +178,7 @@ module mld_d_inner_mod end interface mld_dec_map_bld - interface mld_aggrmat_asb + interface mld_lev_mat_asb subroutine mld_d_lev_aggrmat_asb(p,a,desc_a,ilaggr,nlaggr,op_prol,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ use mld_d_prec_type, only : mld_d_onelev_type @@ -190,9 +190,8 @@ module mld_d_inner_mod type(psb_dspmat_type), intent(inout) :: op_prol integer(psb_ipk_), intent(out) :: info end subroutine mld_d_lev_aggrmat_asb - end interface mld_aggrmat_asb + end interface mld_lev_mat_asb - abstract interface subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info) diff --git a/tests/pdegen/mld_d_pde3d.f90 b/tests/pdegen/mld_d_pde3d.f90 index d6a62892..6a29ed84 100644 --- a/tests/pdegen/mld_d_pde3d.f90 +++ b/tests/pdegen/mld_d_pde3d.f90 @@ -261,6 +261,7 @@ program mld_d_pde3d call mld_precset(prec,'aggr_kind', prectype%aggrkind,info) call mld_precset(prec,'aggr_alg', prectype%aggr_alg,info) call mld_precset(prec,'aggr_ord', prectype%aggr_ord,info) + call mld_precset(prec,'aggr_filter', mld_filter_mat_, info) call psb_barrier(ictxt) t1 = psb_wtime()