From 3c6fc5d14d862f38d0790c5217c81e8befd8adc3 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 20 Jul 2016 11:04:38 +0000 Subject: [PATCH] mld2p4-2: mlprec/impl/mld_caggrmat_biz_asb.f90 mlprec/impl/mld_caggrmat_minnrg_asb.f90 mlprec/impl/mld_cmlprec_aply.f90 mlprec/impl/mld_daggrmat_biz_asb.f90 mlprec/impl/mld_daggrmat_minnrg_asb.f90 mlprec/impl/mld_dmlprec_aply.f90 mlprec/impl/mld_saggrmat_biz_asb.f90 mlprec/impl/mld_saggrmat_minnrg_asb.f90 mlprec/impl/mld_smlprec_aply.f90 mlprec/impl/mld_zaggrmat_biz_asb.f90 mlprec/impl/mld_zaggrmat_minnrg_asb.f90 mlprec/impl/mld_zmlprec_aply.f90 Fixed use of op_prol vs tmp_prol in coarse matrix build. Fixed implementation of KCycle/KCycleSym. Seems to be working now. --- mlprec/impl/mld_caggrmat_biz_asb.f90 | 14 ++--- mlprec/impl/mld_caggrmat_minnrg_asb.f90 | 11 ++-- mlprec/impl/mld_cmlprec_aply.f90 | 83 ++++++++++++++----------- mlprec/impl/mld_daggrmat_biz_asb.f90 | 14 ++--- mlprec/impl/mld_daggrmat_minnrg_asb.f90 | 11 ++-- mlprec/impl/mld_dmlprec_aply.f90 | 83 ++++++++++++++----------- mlprec/impl/mld_saggrmat_biz_asb.f90 | 14 ++--- mlprec/impl/mld_saggrmat_minnrg_asb.f90 | 11 ++-- mlprec/impl/mld_smlprec_aply.f90 | 83 ++++++++++++++----------- mlprec/impl/mld_zaggrmat_biz_asb.f90 | 14 ++--- mlprec/impl/mld_zaggrmat_minnrg_asb.f90 | 11 ++-- mlprec/impl/mld_zmlprec_aply.f90 | 83 ++++++++++++++----------- 12 files changed, 228 insertions(+), 204 deletions(-) diff --git a/mlprec/impl/mld_caggrmat_biz_asb.f90 b/mlprec/impl/mld_caggrmat_biz_asb.f90 index 3b9b7dd4..307ca827 100644 --- a/mlprec/impl/mld_caggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_caggrmat_biz_asb.f90 @@ -97,7 +97,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act integer(psb_ipk_) ::ictxt, np, me character(len=20) :: name - type(psb_cspmat_type) :: am3, am4,tmp_prol + type(psb_cspmat_type) :: am3, am4 type(psb_c_coo_sparse_mat) :: tmpcoo type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_spk_), allocatable :: adiag(:) @@ -350,29 +350,29 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr call acsr1%set_dupl(psb_dupl_add_) call op_prol%mv_from(acsr1) - call op_prol%clone(tmp_prol,info) - call psb_rwextd(ncol,tmp_prol,info) + + call psb_rwextd(ncol,op_prol,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') goto 9999 end if - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ - call tmp_prol%transp(op_restr) + call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call tmp_prol%free() + call psb_rwextd(ncol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') diff --git a/mlprec/impl/mld_caggrmat_minnrg_asb.f90 b/mlprec/impl/mld_caggrmat_minnrg_asb.f90 index 74e15e28..7ebcb689 100644 --- a/mlprec/impl/mld_caggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_caggrmat_minnrg_asb.f90 @@ -120,7 +120,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re character(len=20) :: name type(psb_cspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da - type(psb_cspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol + type(psb_cspmat_type) :: dat, datp, datdatp, atmp3 type(psb_c_coo_sparse_mat) :: tmpcoo type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc @@ -512,10 +512,9 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Now we have to gather the halo of op_prol, and add it to itself ! to multiply it by A, ! - call op_prol%clone(tmp_prol,info) - if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,& + call psb_sphalo(op_prol,desc_a,am4,info,& & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) + if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4) if (info == psb_success_) call am4%free() if(info /= psb_success_) then @@ -548,13 +547,13 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2' diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index f0c66596..79c6c3ff 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -474,7 +474,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,' inner_ml_aply at level ',level + write(debug_unit,*) me,' Start inner_ml_aply at level ',level end if select case(p%precv(level)%parms%ml_type) @@ -538,6 +538,9 @@ contains goto 9999 end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if call psb_erractionrestore(err_act) return @@ -603,7 +606,7 @@ contains call p%precv(level)%sm%apply(cone,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + & sweeps,work,info,init='Z') if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during ADD smoother_apply') @@ -615,7 +618,6 @@ contains call psb_map_X2Y(cone,mlprec_wrk(level)%vx2l,& & czero,mlprec_wrk(level+1)%vx2l,& & p%precv(level+1)%map,info,work=work) - call mlprec_wrk(level+1)%vy2l%zero() if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -917,7 +919,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,name,' at level ',level + write(debug_unit,*) me,name,' start at level ',level end if if ((level<1).or.(level>nlev)) then @@ -937,7 +939,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else if (level < nlev) then @@ -946,13 +948,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(cone,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_post if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then @@ -983,7 +985,7 @@ contains call psb_map_X2Y(cone,mlprec_wrk(level)%vty,& & czero,mlprec_wrk(level + 1)%vx2l,& & p%precv(level + 1)%map,info,work=work) - call mlprec_wrk(level + 1)%vy2l%zero() + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -992,11 +994,11 @@ contains !Set the preconditioner - if (level < nlev ) then + if (level <= nlev - 2 ) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then - call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') + call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR') else call psb_errpush(psb_err_internal_error_,name,& & a_err='Bad value for ml_type') @@ -1091,7 +1093,8 @@ contains type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) integer(psb_ipk_), intent(in) :: level - character, intent(in) :: trans, innersolv + character, intent(in) :: trans + character(len=*), intent(in) :: innersolv complex(psb_spk_),target :: work(:) !Other variables @@ -1099,10 +1102,11 @@ contains type(psb_c_vect_type), dimension(0:1) :: d complex(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta - real(psb_spk_) :: l2_norm, delta, rtol=0.25 + real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm complex(psb_spk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx - + character(len=20) :: name = 'innerit_k_cycle' + !Assemble rhs, w, v, v1, x call psb_geasb(rhs,& @@ -1120,6 +1124,14 @@ contains call psb_geasb(x,& & p%precv(level)%base_desc,info,& & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call psb_geasb(d(1),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call x%zero() @@ -1133,30 +1145,20 @@ contains nc2l = p%precv(level)%base_desc%get_local_cols() info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& - & a_err='TYPE@(psb_spk_)') + & a_err='complex(psb_spk_)') goto 9999 end if - delta = psb_gedot(w, w, p%precv(level)%base_desc, info) + delta0 = psb_genrm2(w, p%precv(level)%base_desc, info) !Apply the preconditioner - - call mlprec_wrk(level)%vy2l%set(czero) + call mlprec_wrk(level)%vy2l%zero() idx=0 call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) - !Assemble d(0) and d(1) - call psb_geasb(d(0),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geaxpby(cone,mlprec_wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info) - call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v,p%precv(level)%base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -1165,24 +1167,27 @@ contains end if !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info) - !CGR - else + !GCR + else if (psb_toupper(trim(innersolv)) == 'GCR') then delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) - tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif alpha = delta_old/tau !Update residual w call psb_geaxpby(-alpha, v, cone, w, p%precv(level)%base_desc, info) - l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info) + l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) iter = 0 - - if (l2_norm <= rtol*delta) then + if (l2_norm <= rtol*delta0) then !Update solution x call psb_geaxpby(alpha, d(idx), cone, x, p%precv(level)%base_desc, info) else @@ -1204,18 +1209,20 @@ contains end if !tau1, tau2, tau3, tau4 - !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau - !CGR - else + else if (psb_toupper(trim(innersolv)) == 'GCR') then tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif !Update solution @@ -1225,8 +1232,8 @@ contains call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info) endif - !Free vectors call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) + !Free vectors call psb_gefree(v, p%precv(level)%base_desc, info) call psb_gefree(v1, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info) diff --git a/mlprec/impl/mld_daggrmat_biz_asb.f90 b/mlprec/impl/mld_daggrmat_biz_asb.f90 index ff631904..90f95927 100644 --- a/mlprec/impl/mld_daggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_daggrmat_biz_asb.f90 @@ -97,7 +97,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act integer(psb_ipk_) ::ictxt, np, me character(len=20) :: name - type(psb_dspmat_type) :: am3, am4,tmp_prol + type(psb_dspmat_type) :: am3, am4 type(psb_d_coo_sparse_mat) :: tmpcoo type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_dpk_), allocatable :: adiag(:) @@ -350,29 +350,29 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr call acsr1%set_dupl(psb_dupl_add_) call op_prol%mv_from(acsr1) - call op_prol%clone(tmp_prol,info) - call psb_rwextd(ncol,tmp_prol,info) + + call psb_rwextd(ncol,op_prol,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') goto 9999 end if - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ - call tmp_prol%transp(op_restr) + call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call tmp_prol%free() + call psb_rwextd(ncol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 index 46989b65..09e7bad6 100644 --- a/mlprec/impl/mld_daggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_daggrmat_minnrg_asb.f90 @@ -120,7 +120,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re character(len=20) :: name type(psb_dspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da - type(psb_dspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol + type(psb_dspmat_type) :: dat, datp, datdatp, atmp3 type(psb_d_coo_sparse_mat) :: tmpcoo type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc @@ -512,10 +512,9 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Now we have to gather the halo of op_prol, and add it to itself ! to multiply it by A, ! - call op_prol%clone(tmp_prol,info) - if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,& + call psb_sphalo(op_prol,desc_a,am4,info,& & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) + if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4) if (info == psb_success_) call am4%free() if(info /= psb_success_) then @@ -548,13 +547,13 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2' diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index b2ba4ad5..ef860cb4 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -474,7 +474,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,' inner_ml_aply at level ',level + write(debug_unit,*) me,' Start inner_ml_aply at level ',level end if select case(p%precv(level)%parms%ml_type) @@ -538,6 +538,9 @@ contains goto 9999 end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if call psb_erractionrestore(err_act) return @@ -603,7 +606,7 @@ contains call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + & sweeps,work,info,init='Z') if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during ADD smoother_apply') @@ -615,7 +618,6 @@ contains call psb_map_X2Y(done,mlprec_wrk(level)%vx2l,& & dzero,mlprec_wrk(level+1)%vx2l,& & p%precv(level+1)%map,info,work=work) - call mlprec_wrk(level+1)%vy2l%zero() if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -917,7 +919,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,name,' at level ',level + write(debug_unit,*) me,name,' start at level ',level end if if ((level<1).or.(level>nlev)) then @@ -937,7 +939,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else if (level < nlev) then @@ -946,13 +948,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_post if (info == psb_success_) call p%precv(level)%sm2%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then @@ -983,7 +985,7 @@ contains call psb_map_X2Y(done,mlprec_wrk(level)%vty,& & dzero,mlprec_wrk(level + 1)%vx2l,& & p%precv(level + 1)%map,info,work=work) - call mlprec_wrk(level + 1)%vy2l%zero() + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -992,11 +994,11 @@ contains !Set the preconditioner - if (level < nlev ) then + if (level <= nlev - 2 ) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then - call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') + call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR') else call psb_errpush(psb_err_internal_error_,name,& & a_err='Bad value for ml_type') @@ -1091,7 +1093,8 @@ contains type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) integer(psb_ipk_), intent(in) :: level - character, intent(in) :: trans, innersolv + character, intent(in) :: trans + character(len=*), intent(in) :: innersolv real(psb_dpk_),target :: work(:) !Other variables @@ -1099,10 +1102,11 @@ contains type(psb_d_vect_type), dimension(0:1) :: d real(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta - real(psb_dpk_) :: l2_norm, delta, rtol=0.25 + real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm real(psb_dpk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx - + character(len=20) :: name = 'innerit_k_cycle' + !Assemble rhs, w, v, v1, x call psb_geasb(rhs,& @@ -1120,6 +1124,14 @@ contains call psb_geasb(x,& & p%precv(level)%base_desc,info,& & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call psb_geasb(d(1),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call x%zero() @@ -1133,30 +1145,20 @@ contains nc2l = p%precv(level)%base_desc%get_local_cols() info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& - & a_err='TYPE@(psb_dpk_)') + & a_err='real(psb_dpk_)') goto 9999 end if - delta = psb_gedot(w, w, p%precv(level)%base_desc, info) + delta0 = psb_genrm2(w, p%precv(level)%base_desc, info) !Apply the preconditioner - - call mlprec_wrk(level)%vy2l%set(dzero) + call mlprec_wrk(level)%vy2l%zero() idx=0 call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) - !Assemble d(0) and d(1) - call psb_geasb(d(0),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info) - call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -1165,24 +1167,27 @@ contains end if !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info) - !CGR - else + !GCR + else if (psb_toupper(trim(innersolv)) == 'GCR') then delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) - tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif alpha = delta_old/tau !Update residual w call psb_geaxpby(-alpha, v, done, w, p%precv(level)%base_desc, info) - l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info) + l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) iter = 0 - - if (l2_norm <= rtol*delta) then + if (l2_norm <= rtol*delta0) then !Update solution x call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info) else @@ -1204,18 +1209,20 @@ contains end if !tau1, tau2, tau3, tau4 - !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau - !CGR - else + else if (psb_toupper(trim(innersolv)) == 'GCR') then tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif !Update solution @@ -1225,8 +1232,8 @@ contains call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info) endif - !Free vectors call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) + !Free vectors call psb_gefree(v, p%precv(level)%base_desc, info) call psb_gefree(v1, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info) diff --git a/mlprec/impl/mld_saggrmat_biz_asb.f90 b/mlprec/impl/mld_saggrmat_biz_asb.f90 index ffc1c65e..3e715cd6 100644 --- a/mlprec/impl/mld_saggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_saggrmat_biz_asb.f90 @@ -97,7 +97,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act integer(psb_ipk_) ::ictxt, np, me character(len=20) :: name - type(psb_sspmat_type) :: am3, am4,tmp_prol + type(psb_sspmat_type) :: am3, am4 type(psb_s_coo_sparse_mat) :: tmpcoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_spk_), allocatable :: adiag(:) @@ -350,29 +350,29 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr call acsr1%set_dupl(psb_dupl_add_) call op_prol%mv_from(acsr1) - call op_prol%clone(tmp_prol,info) - call psb_rwextd(ncol,tmp_prol,info) + + call psb_rwextd(ncol,op_prol,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') goto 9999 end if - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ - call tmp_prol%transp(op_restr) + call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call tmp_prol%free() + call psb_rwextd(ncol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') diff --git a/mlprec/impl/mld_saggrmat_minnrg_asb.f90 b/mlprec/impl/mld_saggrmat_minnrg_asb.f90 index 974645b2..e9e15e4a 100644 --- a/mlprec/impl/mld_saggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_saggrmat_minnrg_asb.f90 @@ -120,7 +120,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re character(len=20) :: name type(psb_sspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da - type(psb_sspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol + type(psb_sspmat_type) :: dat, datp, datdatp, atmp3 type(psb_s_coo_sparse_mat) :: tmpcoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc @@ -512,10 +512,9 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Now we have to gather the halo of op_prol, and add it to itself ! to multiply it by A, ! - call op_prol%clone(tmp_prol,info) - if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,& + call psb_sphalo(op_prol,desc_a,am4,info,& & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) + if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4) if (info == psb_success_) call am4%free() if(info /= psb_success_) then @@ -548,13 +547,13 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2' diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index d5caa998..4a24c246 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -474,7 +474,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,' inner_ml_aply at level ',level + write(debug_unit,*) me,' Start inner_ml_aply at level ',level end if select case(p%precv(level)%parms%ml_type) @@ -538,6 +538,9 @@ contains goto 9999 end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if call psb_erractionrestore(err_act) return @@ -603,7 +606,7 @@ contains call p%precv(level)%sm%apply(sone,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + & sweeps,work,info,init='Z') if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during ADD smoother_apply') @@ -615,7 +618,6 @@ contains call psb_map_X2Y(sone,mlprec_wrk(level)%vx2l,& & szero,mlprec_wrk(level+1)%vx2l,& & p%precv(level+1)%map,info,work=work) - call mlprec_wrk(level+1)%vy2l%zero() if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -917,7 +919,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,name,' at level ',level + write(debug_unit,*) me,name,' start at level ',level end if if ((level<1).or.(level>nlev)) then @@ -937,7 +939,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else if (level < nlev) then @@ -946,13 +948,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(sone,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_post if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then @@ -983,7 +985,7 @@ contains call psb_map_X2Y(sone,mlprec_wrk(level)%vty,& & szero,mlprec_wrk(level + 1)%vx2l,& & p%precv(level + 1)%map,info,work=work) - call mlprec_wrk(level + 1)%vy2l%zero() + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -992,11 +994,11 @@ contains !Set the preconditioner - if (level < nlev ) then + if (level <= nlev - 2 ) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then - call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') + call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR') else call psb_errpush(psb_err_internal_error_,name,& & a_err='Bad value for ml_type') @@ -1091,7 +1093,8 @@ contains type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) integer(psb_ipk_), intent(in) :: level - character, intent(in) :: trans, innersolv + character, intent(in) :: trans + character(len=*), intent(in) :: innersolv real(psb_spk_),target :: work(:) !Other variables @@ -1099,10 +1102,11 @@ contains type(psb_s_vect_type), dimension(0:1) :: d real(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta - real(psb_spk_) :: l2_norm, delta, rtol=0.25 + real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm real(psb_spk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx - + character(len=20) :: name = 'innerit_k_cycle' + !Assemble rhs, w, v, v1, x call psb_geasb(rhs,& @@ -1120,6 +1124,14 @@ contains call psb_geasb(x,& & p%precv(level)%base_desc,info,& & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call psb_geasb(d(1),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call x%zero() @@ -1133,30 +1145,20 @@ contains nc2l = p%precv(level)%base_desc%get_local_cols() info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& - & a_err='TYPE@(psb_spk_)') + & a_err='real(psb_spk_)') goto 9999 end if - delta = psb_gedot(w, w, p%precv(level)%base_desc, info) + delta0 = psb_genrm2(w, p%precv(level)%base_desc, info) !Apply the preconditioner - - call mlprec_wrk(level)%vy2l%set(szero) + call mlprec_wrk(level)%vy2l%zero() idx=0 call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) - !Assemble d(0) and d(1) - call psb_geasb(d(0),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geaxpby(sone,mlprec_wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info) - call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v,p%precv(level)%base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -1165,24 +1167,27 @@ contains end if !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info) - !CGR - else + !GCR + else if (psb_toupper(trim(innersolv)) == 'GCR') then delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) - tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif alpha = delta_old/tau !Update residual w call psb_geaxpby(-alpha, v, sone, w, p%precv(level)%base_desc, info) - l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info) + l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) iter = 0 - - if (l2_norm <= rtol*delta) then + if (l2_norm <= rtol*delta0) then !Update solution x call psb_geaxpby(alpha, d(idx), sone, x, p%precv(level)%base_desc, info) else @@ -1204,18 +1209,20 @@ contains end if !tau1, tau2, tau3, tau4 - !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau - !CGR - else + else if (psb_toupper(trim(innersolv)) == 'GCR') then tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif !Update solution @@ -1225,8 +1232,8 @@ contains call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info) endif - !Free vectors call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) + !Free vectors call psb_gefree(v, p%precv(level)%base_desc, info) call psb_gefree(v1, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info) diff --git a/mlprec/impl/mld_zaggrmat_biz_asb.f90 b/mlprec/impl/mld_zaggrmat_biz_asb.f90 index 7ef0cffd..64bcf643 100644 --- a/mlprec/impl/mld_zaggrmat_biz_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_biz_asb.f90 @@ -97,7 +97,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act integer(psb_ipk_) ::ictxt, np, me character(len=20) :: name - type(psb_zspmat_type) :: am3, am4,tmp_prol + type(psb_zspmat_type) :: am3, am4 type(psb_z_coo_sparse_mat) :: tmpcoo type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_dpk_), allocatable :: adiag(:) @@ -350,29 +350,29 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr call acsr1%set_dupl(psb_dupl_add_) call op_prol%mv_from(acsr1) - call op_prol%clone(tmp_prol,info) - call psb_rwextd(ncol,tmp_prol,info) + + call psb_rwextd(ncol,op_prol,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') goto 9999 end if - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ - call tmp_prol%transp(op_restr) + call op_prol%transp(op_restr) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call tmp_prol%free() + call psb_rwextd(ncol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') diff --git a/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 b/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 index fc5c7c03..ba0238e0 100644 --- a/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_minnrg_asb.f90 @@ -120,7 +120,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re character(len=20) :: name type(psb_zspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da - type(psb_zspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol + type(psb_zspmat_type) :: dat, datp, datdatp, atmp3 type(psb_z_coo_sparse_mat) :: tmpcoo type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc @@ -512,10 +512,9 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re ! Now we have to gather the halo of op_prol, and add it to itself ! to multiply it by A, ! - call op_prol%clone(tmp_prol,info) - if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,& + call psb_sphalo(op_prol,desc_a,am4,info,& & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) + if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4) if (info == psb_success_) call am4%free() if(info /= psb_success_) then @@ -548,13 +547,13 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re & write(debug_unit,*) me,' ',trim(name),& & 'starting sphalo/ rwxtd' - call psb_symbmm(a,tmp_prol,am3,info) + call psb_symbmm(a,op_prol,am3,info) if(info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='symbmm 2') goto 9999 end if - call psb_numbmm(a,tmp_prol,am3) + call psb_numbmm(a,op_prol,am3) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Done NUMBMM 2' diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index f99ba738..1ab04619 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -474,7 +474,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,' inner_ml_aply at level ',level + write(debug_unit,*) me,' Start inner_ml_aply at level ',level end if select case(p%precv(level)%parms%ml_type) @@ -538,6 +538,9 @@ contains goto 9999 end select + if(debug_level > 1) then + write(debug_unit,*) me,' End inner_ml_aply at level ',level + end if call psb_erractionrestore(err_act) return @@ -603,7 +606,7 @@ contains call p%precv(level)%sm%apply(zone,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info) + & sweeps,work,info,init='Z') if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during ADD smoother_apply') @@ -615,7 +618,6 @@ contains call psb_map_X2Y(zone,mlprec_wrk(level)%vx2l,& & zzero,mlprec_wrk(level+1)%vx2l,& & p%precv(level+1)%map,info,work=work) - call mlprec_wrk(level+1)%vy2l%zero() if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -917,7 +919,7 @@ contains call psb_info(ictxt, me, np) if(debug_level > 1) then - write(debug_unit,*) me,name,' at level ',level + write(debug_unit,*) me,name,' start at level ',level end if if ((level<1).or.(level>nlev)) then @@ -937,7 +939,7 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else if (level < nlev) then @@ -946,13 +948,13 @@ contains if (info == psb_success_) call p%precv(level)%sm%apply(zone,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') else sweeps = p%precv(level)%parms%sweeps_post if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& - & sweeps,work,info,init='Y') + & sweeps,work,info,init='Z') end if if (info /= psb_success_) then @@ -983,7 +985,7 @@ contains call psb_map_X2Y(zone,mlprec_wrk(level)%vty,& & zzero,mlprec_wrk(level + 1)%vx2l,& & p%precv(level + 1)%map,info,work=work) - call mlprec_wrk(level + 1)%vy2l%zero() + if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='Error during restriction') @@ -992,11 +994,11 @@ contains !Set the preconditioner - if (level < nlev ) then + if (level <= nlev - 2 ) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then - call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') + call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR') else call psb_errpush(psb_err_internal_error_,name,& & a_err='Bad value for ml_type') @@ -1091,7 +1093,8 @@ contains type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) integer(psb_ipk_), intent(in) :: level - character, intent(in) :: trans, innersolv + character, intent(in) :: trans + character(len=*), intent(in) :: innersolv complex(psb_dpk_),target :: work(:) !Other variables @@ -1099,10 +1102,11 @@ contains type(psb_z_vect_type), dimension(0:1) :: d complex(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta - real(psb_dpk_) :: l2_norm, delta, rtol=0.25 + real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm complex(psb_dpk_), allocatable :: temp_v(:) integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx - + character(len=20) :: name = 'innerit_k_cycle' + !Assemble rhs, w, v, v1, x call psb_geasb(rhs,& @@ -1120,6 +1124,14 @@ contains call psb_geasb(x,& & p%precv(level)%base_desc,info,& & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) + !Assemble d(0) and d(1) + call psb_geasb(d(0),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call psb_geasb(d(1),& + & p%precv(level)%base_desc,info,& + & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) + call x%zero() @@ -1133,30 +1145,20 @@ contains nc2l = p%precv(level)%base_desc%get_local_cols() info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& - & a_err='TYPE@(psb_dpk_)') + & a_err='complex(psb_dpk_)') goto 9999 end if - delta = psb_gedot(w, w, p%precv(level)%base_desc, info) + delta0 = psb_genrm2(w, p%precv(level)%base_desc, info) !Apply the preconditioner - - call mlprec_wrk(level)%vy2l%set(zzero) + call mlprec_wrk(level)%vy2l%zero() idx=0 call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) - !Assemble d(0) and d(1) - call psb_geasb(d(0),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geasb(d(1),& - & p%precv(level)%base_desc,info,& - & scratch=.true.,mold=mlprec_wrk(level)%vy2l%v) - call psb_geaxpby(zone,mlprec_wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info) - call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v,p%precv(level)%base_desc,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -1165,24 +1167,27 @@ contains end if !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info) - !CGR - else + !GCR + else if (psb_toupper(trim(innersolv)) == 'GCR') then delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) - tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + tau = psb_gedot(v, v, p%precv(level)%base_desc, info) + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif alpha = delta_old/tau !Update residual w call psb_geaxpby(-alpha, v, zone, w, p%precv(level)%base_desc, info) - l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info) + l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info) iter = 0 - - if (l2_norm <= rtol*delta) then + if (l2_norm <= rtol*delta0) then !Update solution x call psb_geaxpby(alpha, d(idx), zone, x, p%precv(level)%base_desc, info) else @@ -1204,18 +1209,20 @@ contains end if !tau1, tau2, tau3, tau4 - !FCG - if (innersolv == 'FCG') then + if (psb_toupper(trim(innersolv)) == 'FCG') then tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau - !CGR - else + else if (psb_toupper(trim(innersolv)) == 'GCR') then tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info) tau4= tau2 - (tau1*tau1)/tau + else + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Invalid inner solver') + goto 9999 endif !Update solution @@ -1225,8 +1232,8 @@ contains call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info) endif - !Free vectors call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) + !Free vectors call psb_gefree(v, p%precv(level)%base_desc, info) call psb_gefree(v1, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info)