diff --git a/amgprec/amg_c_base_aggregator_mod.f90 b/amgprec/amg_c_base_aggregator_mod.f90 index 42d7d140..93250ba5 100644 --- a/amgprec/amg_c_base_aggregator_mod.f90 +++ b/amgprec/amg_c_base_aggregator_mod.f90 @@ -126,7 +126,7 @@ module amg_c_base_aggregator_mod & psb_c_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ implicit none type(psb_c_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_c_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr @@ -144,7 +144,7 @@ module amg_c_base_aggregator_mod & psb_c_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ implicit none type(psb_c_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_c_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr diff --git a/amgprec/amg_d_base_aggregator_mod.f90 b/amgprec/amg_d_base_aggregator_mod.f90 index ef87f42e..14e2cd64 100644 --- a/amgprec/amg_d_base_aggregator_mod.f90 +++ b/amgprec/amg_d_base_aggregator_mod.f90 @@ -126,7 +126,7 @@ module amg_d_base_aggregator_mod & psb_d_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ implicit none type(psb_d_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_d_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr @@ -144,7 +144,7 @@ module amg_d_base_aggregator_mod & psb_d_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ implicit none type(psb_d_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_d_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr diff --git a/amgprec/amg_d_parmatch_aggregator_mod.F90 b/amgprec/amg_d_parmatch_aggregator_mod.F90 index fdbf961a..5451e687 100644 --- a/amgprec/amg_d_parmatch_aggregator_mod.F90 +++ b/amgprec/amg_d_parmatch_aggregator_mod.F90 @@ -168,7 +168,7 @@ module amg_d_parmatch_aggregator_mod type(amg_dml_parms), intent(inout) :: parms type(amg_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(inout) :: a - type(psb_desc_type), intent(inout) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_ldspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info @@ -235,7 +235,7 @@ module amg_d_parmatch_aggregator_mod & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data implicit none type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol @@ -257,7 +257,7 @@ module amg_d_parmatch_aggregator_mod integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol - type(psb_dspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_dspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_d_parmatch_unsmth_bld @@ -275,7 +275,7 @@ module amg_d_parmatch_aggregator_mod integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol - type(psb_dspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_dspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_d_parmatch_smth_bld @@ -288,11 +288,11 @@ module amg_d_parmatch_aggregator_mod & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data implicit none type(psb_dspmat_type), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol - type(psb_dspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_dspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(out) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_d_parmatch_spmm_bld_ov @@ -306,11 +306,11 @@ module amg_d_parmatch_aggregator_mod & psb_d_csr_sparse_mat, psb_ld_csr_sparse_mat implicit none type(psb_d_csr_sparse_mat), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol - type(psb_dspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_dspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(out) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_d_parmatch_spmm_bld_inner diff --git a/amgprec/amg_s_base_aggregator_mod.f90 b/amgprec/amg_s_base_aggregator_mod.f90 index f1039639..2c07fc4a 100644 --- a/amgprec/amg_s_base_aggregator_mod.f90 +++ b/amgprec/amg_s_base_aggregator_mod.f90 @@ -126,7 +126,7 @@ module amg_s_base_aggregator_mod & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ implicit none type(psb_s_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr @@ -144,7 +144,7 @@ module amg_s_base_aggregator_mod & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_ implicit none type(psb_s_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr diff --git a/amgprec/amg_s_parmatch_aggregator_mod.F90 b/amgprec/amg_s_parmatch_aggregator_mod.F90 index 32a3b711..fd1defa0 100644 --- a/amgprec/amg_s_parmatch_aggregator_mod.F90 +++ b/amgprec/amg_s_parmatch_aggregator_mod.F90 @@ -168,7 +168,7 @@ module amg_s_parmatch_aggregator_mod type(amg_sml_parms), intent(inout) :: parms type(amg_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(inout) :: a - type(psb_desc_type), intent(inout) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lsspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info @@ -235,7 +235,7 @@ module amg_s_parmatch_aggregator_mod & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data implicit none type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol @@ -257,7 +257,7 @@ module amg_s_parmatch_aggregator_mod integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol - type(psb_sspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_sspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_s_parmatch_unsmth_bld @@ -275,7 +275,7 @@ module amg_s_parmatch_aggregator_mod integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol - type(psb_sspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_sspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_s_parmatch_smth_bld @@ -288,11 +288,11 @@ module amg_s_parmatch_aggregator_mod & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data implicit none type(psb_sspmat_type), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol - type(psb_sspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_sspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(out) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_s_parmatch_spmm_bld_ov @@ -306,11 +306,11 @@ module amg_s_parmatch_aggregator_mod & psb_s_csr_sparse_mat, psb_ls_csr_sparse_mat implicit none type(psb_s_csr_sparse_mat), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol - type(psb_sspmat_type), intent(out) :: op_prol,ac, op_restr + type(psb_sspmat_type), intent(inout) :: op_prol,ac, op_restr type(psb_desc_type), intent(out) :: desc_ac integer(psb_ipk_), intent(out) :: info end subroutine amg_s_parmatch_spmm_bld_inner diff --git a/amgprec/amg_z_base_aggregator_mod.f90 b/amgprec/amg_z_base_aggregator_mod.f90 index 84c52a9d..81858fb7 100644 --- a/amgprec/amg_z_base_aggregator_mod.f90 +++ b/amgprec/amg_z_base_aggregator_mod.f90 @@ -126,7 +126,7 @@ module amg_z_base_aggregator_mod & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ implicit none type(psb_z_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr @@ -144,7 +144,7 @@ module amg_z_base_aggregator_mod & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_ implicit none type(psb_z_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr diff --git a/amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90 index cb8fb6a7..4efaf61d 100644 --- a/amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90 @@ -83,8 +83,8 @@ subroutine amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_c_dec_aggregator_type), target, intent(inout) :: ag type(amg_sml_parms), intent(inout) :: parms type(amg_saggr_data), intent(in) :: ag_data - type(psb_cspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_cspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lcspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 index 9974a233..4a71212e 100644 --- a/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90 @@ -86,8 +86,8 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_c_symdec_aggregator_type), target, intent(inout) :: ag type(amg_sml_parms), intent(inout) :: parms type(amg_saggr_data), intent(in) :: ag_data - type(psb_cspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_cspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lcspmat_type), intent(out) :: op_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 index e14001c6..612c2953 100644 --- a/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90 @@ -105,7 +105,7 @@ ! ! subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,info) + & ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_c_inner_mod, amg_protect_name => amg_caggrmat_minnrg_bld @@ -117,8 +117,8 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms - type(psb_lcspmat_type), intent(inout) :: op_prol - type(psb_lcspmat_type), intent(out) :: ac,op_restr + type(psb_lcspmat_type), intent(inout) :: t_prol + type(psb_cspmat_type), intent(inout) :: op_prol, ac,op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info @@ -171,6 +171,8 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + !NEEDS TO BE REWORKED !! + ! naggr: number of local aggregates ! nrow: local rows. ! @@ -183,361 +185,361 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& goto 9999 end if - ! Get the diagonal D - adiag = a%get_diag(info) - if (info == psb_success_) & - & call psb_realloc(ncol,adiag,info) - if (info == psb_success_) & - & call psb_halo(adiag,desc_a,info) - if (info == psb_success_) call a%cp_to_l(la) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') - goto 9999 - end if - - do i=1,size(adiag) - if (adiag(i) /= czero) then - adinv(i) = cone / adiag(i) - else - adinv(i) = cone - end if - end do - - - - ! 1. Allocate Ptilde in sparse matrix form - call op_prol%mv_to(tmpcoo) - call ptilde%mv_from(tmpcoo) - call ptilde%cscnv(info,type='csr') - - if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call la%cscnv(da,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 - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies done.' - - call da%scal(adinv,info) - - call psb_spspmm(da,ptilde,dap,info) - - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') - goto 9999 - end if - - call dap%clone(atmp,info) - - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) - if (info == psb_success_) call am4%free() - - call psb_spspmm(da,atmp,dadap,info) - call atmp%free() - - ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) - ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) - call dap%mv_to(csc_dap) - call dadap%mv_to(csc_dadap) - - call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) - call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - ! !$ write(0,*) trim(name),' OMP :',omp - ! !$ write(0,*) trim(name),' ODEN:',oden - - omp = omp/oden - - ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - call am3%mv_to(acsr3) - ! Compute omega_int - ommx = czero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = czero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - do i=1, nrow - omf(i) = ommx - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero - if(psb_minreal(omf(i)) < szero) omf(i) = czero - end do - - omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) - - if (filter_mat) then - ! - ! Build the filtered matrix Af from A - ! - call la%cscnv(acsrf,info,dupl=psb_dupl_add_) - - do i=1,nrow - tmp = czero - jd = -1 - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) jd = j - if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then - tmp=tmp+acsrf%val(j) - acsrf%val(j)=czero - endif - enddo - if (jd == -1) then - write(0,*) 'Wrong input: we need the diagonal!!!!', i - else - acsrf%val(jd)=acsrf%val(jd)-tmp - end if - enddo - ! Take out zeroed terms - call acsrf%clean_zeros(info) - - ! - ! Build the smoothed prolongator using the filtered matrix - ! - do i=1,acsrf%get_nrows() - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) then - acsrf%val(j) = cone - omf(i)*acsrf%val(j) - else - acsrf%val(j) = - omf(i)*acsrf%val(j) - end if - end do - end do - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - - call af%mv_from(acsrf) - ! - ! op_prol = (I-w*D*Af)Ptilde - ! Doing it this way means to consider diag(Af_i) - ! - ! - call psb_spspmm(af,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 1' - else - ! - ! Build the smoothed prolongator using the original matrix - ! - do i=1,acsr3%get_nrows() - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if (acsr3%ja(j) == i) then - acsr3%val(j) = cone - omf(i)*acsr3%val(j) - else - acsr3%val(j) = - omf(i)*acsr3%val(j) - end if - end do - end do - - call am3%mv_from(acsr3) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - ! - ! - ! op_prol = (I-w*D*A)Ptilde - ! - ! - call psb_spspmm(am3,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - end if - - - ! - ! Ok, let's start over with the restrictor - ! - call ptilde%transc(rtilde) - call la%cscnv(atmp,info,type='csr') - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.true.,rowscale=.true.) - nrt = am4%get_nrows() - call am4%csclip(atmp2,info,lone,nrt,lone,ncol) - call atmp2%cscnv(info,type='CSR') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) - call am4%free() - call atmp2%free() - - ! This is to compute the transpose. It ONLY works if the - ! original A has a symmetric pattern. - call atmp%transc(atmp2) - call atmp2%csclip(dat,info,lone,nrow,lone,ncol) - call dat%cscnv(info,type='csr') - call dat%scal(adinv,info) - - ! Now for the product. - call psb_spspmm(dat,ptilde,datp,info) - - call datp%clone(atmp2,info) - call psb_sphalo(atmp2,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) - if (info == psb_success_) call am4%free() - - - call psb_symbmm(dat,atmp2,datdatp,info) - call psb_numbmm(dat,atmp2,datdatp) - call atmp2%free() - - call datp%mv_to(csc_datp) - call datdatp%mv_to(csc_datdatp) - - call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) - call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - - - ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp - ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden - omp = omp/oden - ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) - ! Compute omega_int - ommx = czero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = czero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - ! Going over the columns of atmp means going over the rows - ! of A^T. Hopefully ;-) - call atmp%cp_to(acsc) - - do i=1, nrow - omf(i) = ommx - do j= acsc%icp(i),acsc%icp(i+1)-1 - if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero - if(psb_minreal(omf(i)) < szero) omf(i) = czero - end do - omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) - call psb_halo(omf,desc_a,info) - call acsc%free() - - - call atmp%mv_to(acsr1) - - do i=1,acsr1%get_nrows() - do j=acsr1%irp(i),acsr1%irp(i+1)-1 - if (acsr1%ja(j) == i) then - acsr1%val(j) = cone - acsr1%val(j)*omf(acsr1%ja(j)) - else - acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) - end if - end do - end do - call atmp%mv_from(acsr1) - - call rtilde%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call rtilde%mv_from(tmpcoo) - call rtilde%cscnv(info,type='csr') - - call psb_spspmm(rtilde,atmp,op_restr,info) - - ! - ! 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,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') - goto 9999 - end if - - ! - ! Now we have to fix this. The only rows of B that are correct - ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) - ! - call op_restr%mv_to(tmpcoo) - - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - call psb_spspmm(la,tmp_prol,am3,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 2' - - call psb_sphalo(am3,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Extend am3') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done sphalo/ rwxtd' - - call psb_spspmm(op_restr,am3,ac,info) - if (info == psb_success_) call am3%free() - if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build ac = op_restr x am3') - goto 9999 - end if +!!$ ! Get the diagonal D +!!$ adiag = a%get_diag(info) +!!$ if (info == psb_success_) & +!!$ & call psb_realloc(ncol,adiag,info) +!!$ if (info == psb_success_) & +!!$ & call psb_halo(adiag,desc_a,info) +!!$ if (info == psb_success_) call a%cp_to_l(la) +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') +!!$ goto 9999 +!!$ end if +!!$ +!!$ do i=1,size(adiag) +!!$ if (adiag(i) /= czero) then +!!$ adinv(i) = cone / adiag(i) +!!$ else +!!$ adinv(i) = cone +!!$ end if +!!$ end do +!!$ +!!$ +!!$ +!!$ ! 1. Allocate Ptilde in sparse matrix form +!!$ call op_prol%mv_to(tmpcoo) +!!$ call ptilde%mv_from(tmpcoo) +!!$ call ptilde%cscnv(info,type='csr') +!!$ +!!$ if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) +!!$ if (info == psb_success_) call la%cscnv(da,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 +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & ' Initial copies done.' +!!$ +!!$ call da%scal(adinv,info) +!!$ +!!$ call psb_spspmm(da,ptilde,dap,info) +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ call dap%clone(atmp,info) +!!$ +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ call psb_spspmm(da,atmp,dadap,info) +!!$ call atmp%free() +!!$ +!!$ ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) +!!$ ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) +!!$ call dap%mv_to(csc_dap) +!!$ call dadap%mv_to(csc_dadap) +!!$ +!!$ call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) +!!$ call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ ! !$ write(0,*) trim(name),' OMP :',omp +!!$ ! !$ write(0,*) trim(name),' ODEN:',oden +!!$ +!!$ omp = omp/oden +!!$ +!!$ ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ call am3%mv_to(acsr3) +!!$ ! Compute omega_int +!!$ ommx = czero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = czero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero +!!$ if(psb_minreal(omf(i)) < szero) omf(i) = czero +!!$ end do +!!$ +!!$ omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) +!!$ +!!$ if (filter_mat) then +!!$ ! +!!$ ! Build the filtered matrix Af from A +!!$ ! +!!$ call la%cscnv(acsrf,info,dupl=psb_dupl_add_) +!!$ +!!$ do i=1,nrow +!!$ tmp = czero +!!$ jd = -1 +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) jd = j +!!$ if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then +!!$ tmp=tmp+acsrf%val(j) +!!$ acsrf%val(j)=czero +!!$ endif +!!$ enddo +!!$ if (jd == -1) then +!!$ write(0,*) 'Wrong input: we need the diagonal!!!!', i +!!$ else +!!$ acsrf%val(jd)=acsrf%val(jd)-tmp +!!$ end if +!!$ enddo +!!$ ! Take out zeroed terms +!!$ call acsrf%clean_zeros(info) +!!$ +!!$ ! +!!$ ! Build the smoothed prolongator using the filtered matrix +!!$ ! +!!$ do i=1,acsrf%get_nrows() +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) then +!!$ acsrf%val(j) = cone - omf(i)*acsrf%val(j) +!!$ else +!!$ acsrf%val(j) = - omf(i)*acsrf%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ +!!$ call af%mv_from(acsrf) +!!$ ! +!!$ ! op_prol = (I-w*D*Af)Ptilde +!!$ ! Doing it this way means to consider diag(Af_i) +!!$ ! +!!$ ! +!!$ call psb_spspmm(af,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 1' +!!$ else +!!$ ! +!!$ ! Build the smoothed prolongator using the original matrix +!!$ ! +!!$ do i=1,acsr3%get_nrows() +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if (acsr3%ja(j) == i) then +!!$ acsr3%val(j) = cone - omf(i)*acsr3%val(j) +!!$ else +!!$ acsr3%val(j) = - omf(i)*acsr3%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ call am3%mv_from(acsr3) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ ! +!!$ ! +!!$ ! op_prol = (I-w*D*A)Ptilde +!!$ ! +!!$ ! +!!$ call psb_spspmm(am3,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ end if +!!$ +!!$ +!!$ ! +!!$ ! Ok, let's start over with the restrictor +!!$ ! +!!$ call ptilde%transc(rtilde) +!!$ call la%cscnv(atmp,info,type='csr') +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.true.,rowscale=.true.) +!!$ nrt = am4%get_nrows() +!!$ call am4%csclip(atmp2,info,lone,nrt,lone,ncol) +!!$ call atmp2%cscnv(info,type='CSR') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) +!!$ call am4%free() +!!$ call atmp2%free() +!!$ +!!$ ! This is to compute the transpose. It ONLY works if the +!!$ ! original A has a symmetric pattern. +!!$ call atmp%transc(atmp2) +!!$ call atmp2%csclip(dat,info,lone,nrow,lone,ncol) +!!$ call dat%cscnv(info,type='csr') +!!$ call dat%scal(adinv,info) +!!$ +!!$ ! Now for the product. +!!$ call psb_spspmm(dat,ptilde,datp,info) +!!$ +!!$ call datp%clone(atmp2,info) +!!$ call psb_sphalo(atmp2,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ +!!$ call psb_symbmm(dat,atmp2,datdatp,info) +!!$ call psb_numbmm(dat,atmp2,datdatp) +!!$ call atmp2%free() +!!$ +!!$ call datp%mv_to(csc_datp) +!!$ call datdatp%mv_to(csc_datdatp) +!!$ +!!$ call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) +!!$ call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ +!!$ +!!$ ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp +!!$ ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden +!!$ omp = omp/oden +!!$ ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) +!!$ ! Compute omega_int +!!$ ommx = czero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = czero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ ! Going over the columns of atmp means going over the rows +!!$ ! of A^T. Hopefully ;-) +!!$ call atmp%cp_to(acsc) +!!$ +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j= acsc%icp(i),acsc%icp(i+1)-1 +!!$ if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero +!!$ if(psb_minreal(omf(i)) < szero) omf(i) = czero +!!$ end do +!!$ omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) +!!$ call psb_halo(omf,desc_a,info) +!!$ call acsc%free() +!!$ +!!$ +!!$ call atmp%mv_to(acsr1) +!!$ +!!$ do i=1,acsr1%get_nrows() +!!$ do j=acsr1%irp(i),acsr1%irp(i+1)-1 +!!$ if (acsr1%ja(j) == i) then +!!$ acsr1%val(j) = cone - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ else +!!$ acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ end if +!!$ end do +!!$ end do +!!$ call atmp%mv_from(acsr1) +!!$ +!!$ call rtilde%mv_to(tmpcoo) +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call rtilde%mv_from(tmpcoo) +!!$ call rtilde%cscnv(info,type='csr') +!!$ +!!$ call psb_spspmm(rtilde,atmp,op_restr,info) +!!$ +!!$ ! +!!$ ! 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,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') +!!$ goto 9999 +!!$ end if +!!$ +!!$ ! +!!$ ! Now we have to fix this. The only rows of B that are correct +!!$ ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) +!!$ ! +!!$ call op_restr%mv_to(tmpcoo) +!!$ +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call op_restr%mv_from(tmpcoo) +!!$ call op_restr%cscnv(info,type='csr') +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(la,tmp_prol,am3,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 2' +!!$ +!!$ call psb_sphalo(am3,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(op_restr,am3,ac,info) +!!$ if (info == psb_success_) call am3%free() +!!$ if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ &a_err='Build ac = op_restr x am3') +!!$ goto 9999 +!!$ end if diff --git a/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 index 03029f40..53e740fe 100644 --- a/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90 @@ -116,7 +116,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms - type(psb_cspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_cspmat_type), intent(inout) :: op_prol,ac,op_restr type(psb_lcspmat_type), intent(inout) :: t_prol type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_d_dec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_d_dec_aggregator_tprol.f90 index e3d1e73c..2edcca6c 100644 --- a/amgprec/impl/aggregator/amg_d_dec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_d_dec_aggregator_tprol.f90 @@ -83,8 +83,8 @@ subroutine amg_d_dec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_d_dec_aggregator_type), target, intent(inout) :: ag type(amg_dml_parms), intent(inout) :: parms type(amg_daggr_data), intent(in) :: ag_data - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_ldspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_d_parmatch_aggregator_tprol.F90 b/amgprec/impl/aggregator/amg_d_parmatch_aggregator_tprol.F90 index 8559fae0..16cc7b2e 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_aggregator_tprol.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_aggregator_tprol.F90 @@ -58,7 +58,7 @@ subroutine amg_d_parmatch_aggregator_build_tprol(ag,parms,ag_data,& type(amg_dml_parms), intent(inout) :: parms type(amg_daggr_data), intent(in) :: ag_data type(psb_dspmat_type), intent(inout) :: a - type(psb_desc_type), intent(inout) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_ldspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 b/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 index 80022b73..2ad6388d 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90 @@ -122,7 +122,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol - type(psb_dspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol,ac,op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld.F90 b/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld.F90 index 0311f274..7faf408a 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld.F90 @@ -108,7 +108,7 @@ subroutine amg_d_parmatch_spmm_bld(a,desc_a,ilaggr,nlaggr,parms,& ! Arguments type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol diff --git a/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_inner.F90 b/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_inner.F90 index a2917e83..04d89b2f 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_inner.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_inner.F90 @@ -108,11 +108,11 @@ subroutine amg_d_parmatch_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,& ! Arguments type(psb_d_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol - type(psb_dspmat_type), intent(out) :: ac, op_prol, op_restr + type(psb_dspmat_type), intent(inout) :: ac, op_prol, op_restr type(psb_desc_type), intent(out) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_ov.F90 b/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_ov.F90 index 2b35c797..e3d3b262 100644 --- a/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_ov.F90 +++ b/amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_ov.F90 @@ -108,7 +108,7 @@ subroutine amg_d_parmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,& ! Arguments type(psb_dspmat_type), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms type(psb_ldspmat_type), intent(inout) :: t_prol diff --git a/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 index 35efe63d..e4b2fd6e 100644 --- a/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90 @@ -86,8 +86,8 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_d_symdec_aggregator_type), target, intent(inout) :: ag type(amg_dml_parms), intent(inout) :: parms type(amg_daggr_data), intent(in) :: ag_data - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_ldspmat_type), intent(out) :: op_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 index fc5728a6..f7b21847 100644 --- a/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90 @@ -105,7 +105,7 @@ ! ! subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,info) + & ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_d_inner_mod, amg_protect_name => amg_daggrmat_minnrg_bld @@ -117,8 +117,8 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms - type(psb_ldspmat_type), intent(inout) :: op_prol - type(psb_ldspmat_type), intent(out) :: ac,op_restr + type(psb_ldspmat_type), intent(inout) :: t_prol + type(psb_dspmat_type), intent(inout) :: op_prol, ac,op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info @@ -171,6 +171,8 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + !NEEDS TO BE REWORKED !! + ! naggr: number of local aggregates ! nrow: local rows. ! @@ -183,361 +185,361 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& goto 9999 end if - ! Get the diagonal D - adiag = a%get_diag(info) - if (info == psb_success_) & - & call psb_realloc(ncol,adiag,info) - if (info == psb_success_) & - & call psb_halo(adiag,desc_a,info) - if (info == psb_success_) call a%cp_to_l(la) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') - goto 9999 - end if - - do i=1,size(adiag) - if (adiag(i) /= dzero) then - adinv(i) = done / adiag(i) - else - adinv(i) = done - end if - end do - - - - ! 1. Allocate Ptilde in sparse matrix form - call op_prol%mv_to(tmpcoo) - call ptilde%mv_from(tmpcoo) - call ptilde%cscnv(info,type='csr') - - if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call la%cscnv(da,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 - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies done.' - - call da%scal(adinv,info) - - call psb_spspmm(da,ptilde,dap,info) - - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') - goto 9999 - end if - - call dap%clone(atmp,info) - - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) - if (info == psb_success_) call am4%free() - - call psb_spspmm(da,atmp,dadap,info) - call atmp%free() - - ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) - ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) - call dap%mv_to(csc_dap) - call dadap%mv_to(csc_dadap) - - call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) - call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - ! !$ write(0,*) trim(name),' OMP :',omp - ! !$ write(0,*) trim(name),' ODEN:',oden - - omp = omp/oden - - ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - call am3%mv_to(acsr3) - ! Compute omega_int - ommx = dzero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = dzero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - do i=1, nrow - omf(i) = ommx - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero - if(psb_minreal(omf(i)) < dzero) omf(i) = dzero - end do - - omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) - - if (filter_mat) then - ! - ! Build the filtered matrix Af from A - ! - call la%cscnv(acsrf,info,dupl=psb_dupl_add_) - - do i=1,nrow - tmp = dzero - jd = -1 - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) jd = j - if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then - tmp=tmp+acsrf%val(j) - acsrf%val(j)=dzero - endif - enddo - if (jd == -1) then - write(0,*) 'Wrong input: we need the diagonal!!!!', i - else - acsrf%val(jd)=acsrf%val(jd)-tmp - end if - enddo - ! Take out zeroed terms - call acsrf%clean_zeros(info) - - ! - ! Build the smoothed prolongator using the filtered matrix - ! - do i=1,acsrf%get_nrows() - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) then - acsrf%val(j) = done - omf(i)*acsrf%val(j) - else - acsrf%val(j) = - omf(i)*acsrf%val(j) - end if - end do - end do - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - - call af%mv_from(acsrf) - ! - ! op_prol = (I-w*D*Af)Ptilde - ! Doing it this way means to consider diag(Af_i) - ! - ! - call psb_spspmm(af,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 1' - else - ! - ! Build the smoothed prolongator using the original matrix - ! - do i=1,acsr3%get_nrows() - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if (acsr3%ja(j) == i) then - acsr3%val(j) = done - omf(i)*acsr3%val(j) - else - acsr3%val(j) = - omf(i)*acsr3%val(j) - end if - end do - end do - - call am3%mv_from(acsr3) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - ! - ! - ! op_prol = (I-w*D*A)Ptilde - ! - ! - call psb_spspmm(am3,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - end if - - - ! - ! Ok, let's start over with the restrictor - ! - call ptilde%transc(rtilde) - call la%cscnv(atmp,info,type='csr') - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.true.,rowscale=.true.) - nrt = am4%get_nrows() - call am4%csclip(atmp2,info,lone,nrt,lone,ncol) - call atmp2%cscnv(info,type='CSR') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) - call am4%free() - call atmp2%free() - - ! This is to compute the transpose. It ONLY works if the - ! original A has a symmetric pattern. - call atmp%transc(atmp2) - call atmp2%csclip(dat,info,lone,nrow,lone,ncol) - call dat%cscnv(info,type='csr') - call dat%scal(adinv,info) - - ! Now for the product. - call psb_spspmm(dat,ptilde,datp,info) - - call datp%clone(atmp2,info) - call psb_sphalo(atmp2,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) - if (info == psb_success_) call am4%free() - - - call psb_symbmm(dat,atmp2,datdatp,info) - call psb_numbmm(dat,atmp2,datdatp) - call atmp2%free() - - call datp%mv_to(csc_datp) - call datdatp%mv_to(csc_datdatp) - - call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) - call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - - - ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp - ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden - omp = omp/oden - ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) - ! Compute omega_int - ommx = dzero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = dzero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - ! Going over the columns of atmp means going over the rows - ! of A^T. Hopefully ;-) - call atmp%cp_to(acsc) - - do i=1, nrow - omf(i) = ommx - do j= acsc%icp(i),acsc%icp(i+1)-1 - if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero - if(psb_minreal(omf(i)) < dzero) omf(i) = dzero - end do - omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) - call psb_halo(omf,desc_a,info) - call acsc%free() - - - call atmp%mv_to(acsr1) - - do i=1,acsr1%get_nrows() - do j=acsr1%irp(i),acsr1%irp(i+1)-1 - if (acsr1%ja(j) == i) then - acsr1%val(j) = done - acsr1%val(j)*omf(acsr1%ja(j)) - else - acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) - end if - end do - end do - call atmp%mv_from(acsr1) - - call rtilde%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call rtilde%mv_from(tmpcoo) - call rtilde%cscnv(info,type='csr') - - call psb_spspmm(rtilde,atmp,op_restr,info) - - ! - ! 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,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') - goto 9999 - end if - - ! - ! Now we have to fix this. The only rows of B that are correct - ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) - ! - call op_restr%mv_to(tmpcoo) - - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - call psb_spspmm(la,tmp_prol,am3,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 2' - - call psb_sphalo(am3,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Extend am3') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done sphalo/ rwxtd' - - call psb_spspmm(op_restr,am3,ac,info) - if (info == psb_success_) call am3%free() - if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build ac = op_restr x am3') - goto 9999 - end if +!!$ ! Get the diagonal D +!!$ adiag = a%get_diag(info) +!!$ if (info == psb_success_) & +!!$ & call psb_realloc(ncol,adiag,info) +!!$ if (info == psb_success_) & +!!$ & call psb_halo(adiag,desc_a,info) +!!$ if (info == psb_success_) call a%cp_to_l(la) +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') +!!$ goto 9999 +!!$ end if +!!$ +!!$ do i=1,size(adiag) +!!$ if (adiag(i) /= dzero) then +!!$ adinv(i) = done / adiag(i) +!!$ else +!!$ adinv(i) = done +!!$ end if +!!$ end do +!!$ +!!$ +!!$ +!!$ ! 1. Allocate Ptilde in sparse matrix form +!!$ call op_prol%mv_to(tmpcoo) +!!$ call ptilde%mv_from(tmpcoo) +!!$ call ptilde%cscnv(info,type='csr') +!!$ +!!$ if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) +!!$ if (info == psb_success_) call la%cscnv(da,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 +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & ' Initial copies done.' +!!$ +!!$ call da%scal(adinv,info) +!!$ +!!$ call psb_spspmm(da,ptilde,dap,info) +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ call dap%clone(atmp,info) +!!$ +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ call psb_spspmm(da,atmp,dadap,info) +!!$ call atmp%free() +!!$ +!!$ ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) +!!$ ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) +!!$ call dap%mv_to(csc_dap) +!!$ call dadap%mv_to(csc_dadap) +!!$ +!!$ call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) +!!$ call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ ! !$ write(0,*) trim(name),' OMP :',omp +!!$ ! !$ write(0,*) trim(name),' ODEN:',oden +!!$ +!!$ omp = omp/oden +!!$ +!!$ ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ call am3%mv_to(acsr3) +!!$ ! Compute omega_int +!!$ ommx = dzero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = dzero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero +!!$ if(psb_minreal(omf(i)) < dzero) omf(i) = dzero +!!$ end do +!!$ +!!$ omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) +!!$ +!!$ if (filter_mat) then +!!$ ! +!!$ ! Build the filtered matrix Af from A +!!$ ! +!!$ call la%cscnv(acsrf,info,dupl=psb_dupl_add_) +!!$ +!!$ do i=1,nrow +!!$ tmp = dzero +!!$ jd = -1 +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) jd = j +!!$ if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then +!!$ tmp=tmp+acsrf%val(j) +!!$ acsrf%val(j)=dzero +!!$ endif +!!$ enddo +!!$ if (jd == -1) then +!!$ write(0,*) 'Wrong input: we need the diagonal!!!!', i +!!$ else +!!$ acsrf%val(jd)=acsrf%val(jd)-tmp +!!$ end if +!!$ enddo +!!$ ! Take out zeroed terms +!!$ call acsrf%clean_zeros(info) +!!$ +!!$ ! +!!$ ! Build the smoothed prolongator using the filtered matrix +!!$ ! +!!$ do i=1,acsrf%get_nrows() +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) then +!!$ acsrf%val(j) = done - omf(i)*acsrf%val(j) +!!$ else +!!$ acsrf%val(j) = - omf(i)*acsrf%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ +!!$ call af%mv_from(acsrf) +!!$ ! +!!$ ! op_prol = (I-w*D*Af)Ptilde +!!$ ! Doing it this way means to consider diag(Af_i) +!!$ ! +!!$ ! +!!$ call psb_spspmm(af,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 1' +!!$ else +!!$ ! +!!$ ! Build the smoothed prolongator using the original matrix +!!$ ! +!!$ do i=1,acsr3%get_nrows() +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if (acsr3%ja(j) == i) then +!!$ acsr3%val(j) = done - omf(i)*acsr3%val(j) +!!$ else +!!$ acsr3%val(j) = - omf(i)*acsr3%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ call am3%mv_from(acsr3) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ ! +!!$ ! +!!$ ! op_prol = (I-w*D*A)Ptilde +!!$ ! +!!$ ! +!!$ call psb_spspmm(am3,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ end if +!!$ +!!$ +!!$ ! +!!$ ! Ok, let's start over with the restrictor +!!$ ! +!!$ call ptilde%transc(rtilde) +!!$ call la%cscnv(atmp,info,type='csr') +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.true.,rowscale=.true.) +!!$ nrt = am4%get_nrows() +!!$ call am4%csclip(atmp2,info,lone,nrt,lone,ncol) +!!$ call atmp2%cscnv(info,type='CSR') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) +!!$ call am4%free() +!!$ call atmp2%free() +!!$ +!!$ ! This is to compute the transpose. It ONLY works if the +!!$ ! original A has a symmetric pattern. +!!$ call atmp%transc(atmp2) +!!$ call atmp2%csclip(dat,info,lone,nrow,lone,ncol) +!!$ call dat%cscnv(info,type='csr') +!!$ call dat%scal(adinv,info) +!!$ +!!$ ! Now for the product. +!!$ call psb_spspmm(dat,ptilde,datp,info) +!!$ +!!$ call datp%clone(atmp2,info) +!!$ call psb_sphalo(atmp2,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ +!!$ call psb_symbmm(dat,atmp2,datdatp,info) +!!$ call psb_numbmm(dat,atmp2,datdatp) +!!$ call atmp2%free() +!!$ +!!$ call datp%mv_to(csc_datp) +!!$ call datdatp%mv_to(csc_datdatp) +!!$ +!!$ call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) +!!$ call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ +!!$ +!!$ ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp +!!$ ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden +!!$ omp = omp/oden +!!$ ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) +!!$ ! Compute omega_int +!!$ ommx = dzero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = dzero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ ! Going over the columns of atmp means going over the rows +!!$ ! of A^T. Hopefully ;-) +!!$ call atmp%cp_to(acsc) +!!$ +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j= acsc%icp(i),acsc%icp(i+1)-1 +!!$ if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero +!!$ if(psb_minreal(omf(i)) < dzero) omf(i) = dzero +!!$ end do +!!$ omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) +!!$ call psb_halo(omf,desc_a,info) +!!$ call acsc%free() +!!$ +!!$ +!!$ call atmp%mv_to(acsr1) +!!$ +!!$ do i=1,acsr1%get_nrows() +!!$ do j=acsr1%irp(i),acsr1%irp(i+1)-1 +!!$ if (acsr1%ja(j) == i) then +!!$ acsr1%val(j) = done - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ else +!!$ acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ end if +!!$ end do +!!$ end do +!!$ call atmp%mv_from(acsr1) +!!$ +!!$ call rtilde%mv_to(tmpcoo) +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call rtilde%mv_from(tmpcoo) +!!$ call rtilde%cscnv(info,type='csr') +!!$ +!!$ call psb_spspmm(rtilde,atmp,op_restr,info) +!!$ +!!$ ! +!!$ ! 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,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') +!!$ goto 9999 +!!$ end if +!!$ +!!$ ! +!!$ ! Now we have to fix this. The only rows of B that are correct +!!$ ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) +!!$ ! +!!$ call op_restr%mv_to(tmpcoo) +!!$ +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call op_restr%mv_from(tmpcoo) +!!$ call op_restr%cscnv(info,type='csr') +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(la,tmp_prol,am3,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 2' +!!$ +!!$ call psb_sphalo(am3,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(op_restr,am3,ac,info) +!!$ if (info == psb_success_) call am3%free() +!!$ if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ &a_err='Build ac = op_restr x am3') +!!$ goto 9999 +!!$ end if diff --git a/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 index 20b24d60..82da3fc7 100644 --- a/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90 @@ -116,7 +116,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms - type(psb_dspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_dspmat_type), intent(inout) :: op_prol,ac,op_restr type(psb_ldspmat_type), intent(inout) :: t_prol type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_s_dec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_s_dec_aggregator_tprol.f90 index 0ab5274e..c52c04f7 100644 --- a/amgprec/impl/aggregator/amg_s_dec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_s_dec_aggregator_tprol.f90 @@ -83,8 +83,8 @@ subroutine amg_s_dec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_s_dec_aggregator_type), target, intent(inout) :: ag type(amg_sml_parms), intent(inout) :: parms type(amg_saggr_data), intent(in) :: ag_data - type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_sspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lsspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_s_parmatch_aggregator_tprol.F90 b/amgprec/impl/aggregator/amg_s_parmatch_aggregator_tprol.F90 index 2cb3d842..cd28893c 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_aggregator_tprol.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_aggregator_tprol.F90 @@ -58,7 +58,7 @@ subroutine amg_s_parmatch_aggregator_build_tprol(ag,parms,ag_data,& type(amg_sml_parms), intent(inout) :: parms type(amg_saggr_data), intent(in) :: ag_data type(psb_sspmat_type), intent(inout) :: a - type(psb_desc_type), intent(inout) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lsspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 b/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 index a89ac411..a3f79fc5 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90 @@ -122,7 +122,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,& integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol - type(psb_sspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_sspmat_type), intent(inout) :: op_prol,ac,op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld.F90 b/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld.F90 index b88930eb..3f9e385e 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld.F90 @@ -108,7 +108,7 @@ subroutine amg_s_parmatch_spmm_bld(a,desc_a,ilaggr,nlaggr,parms,& ! Arguments type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol diff --git a/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_inner.F90 b/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_inner.F90 index cb64c0e1..9f1c81f1 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_inner.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_inner.F90 @@ -108,11 +108,11 @@ subroutine amg_s_parmatch_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,& ! Arguments type(psb_s_csr_sparse_mat), intent(inout) :: a_csr - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol - type(psb_sspmat_type), intent(out) :: ac, op_prol, op_restr + type(psb_sspmat_type), intent(inout) :: ac, op_prol, op_restr type(psb_desc_type), intent(out) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_ov.F90 b/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_ov.F90 index a4fa85e7..442f4225 100644 --- a/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_ov.F90 +++ b/amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_ov.F90 @@ -108,7 +108,7 @@ subroutine amg_s_parmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,& ! Arguments type(psb_sspmat_type), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms type(psb_lsspmat_type), intent(inout) :: t_prol diff --git a/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 index 5a9548eb..c148b89a 100644 --- a/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90 @@ -86,8 +86,8 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_s_symdec_aggregator_type), target, intent(inout) :: ag type(amg_sml_parms), intent(inout) :: parms type(amg_saggr_data), intent(in) :: ag_data - type(psb_sspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_sspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lsspmat_type), intent(out) :: op_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 index 670617fe..66514a66 100644 --- a/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90 @@ -105,7 +105,7 @@ ! ! subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,info) + & ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_s_inner_mod, amg_protect_name => amg_saggrmat_minnrg_bld @@ -117,8 +117,8 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms - type(psb_lsspmat_type), intent(inout) :: op_prol - type(psb_lsspmat_type), intent(out) :: ac,op_restr + type(psb_lsspmat_type), intent(inout) :: t_prol + type(psb_sspmat_type), intent(inout) :: op_prol, ac,op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info @@ -171,6 +171,8 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + !NEEDS TO BE REWORKED !! + ! naggr: number of local aggregates ! nrow: local rows. ! @@ -183,361 +185,361 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& goto 9999 end if - ! Get the diagonal D - adiag = a%get_diag(info) - if (info == psb_success_) & - & call psb_realloc(ncol,adiag,info) - if (info == psb_success_) & - & call psb_halo(adiag,desc_a,info) - if (info == psb_success_) call a%cp_to_l(la) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') - goto 9999 - end if - - do i=1,size(adiag) - if (adiag(i) /= szero) then - adinv(i) = sone / adiag(i) - else - adinv(i) = sone - end if - end do - - - - ! 1. Allocate Ptilde in sparse matrix form - call op_prol%mv_to(tmpcoo) - call ptilde%mv_from(tmpcoo) - call ptilde%cscnv(info,type='csr') - - if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call la%cscnv(da,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 - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies done.' - - call da%scal(adinv,info) - - call psb_spspmm(da,ptilde,dap,info) - - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') - goto 9999 - end if - - call dap%clone(atmp,info) - - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) - if (info == psb_success_) call am4%free() - - call psb_spspmm(da,atmp,dadap,info) - call atmp%free() - - ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) - ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) - call dap%mv_to(csc_dap) - call dadap%mv_to(csc_dadap) - - call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) - call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - ! !$ write(0,*) trim(name),' OMP :',omp - ! !$ write(0,*) trim(name),' ODEN:',oden - - omp = omp/oden - - ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - call am3%mv_to(acsr3) - ! Compute omega_int - ommx = szero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = szero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - do i=1, nrow - omf(i) = ommx - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero - if(psb_minreal(omf(i)) < szero) omf(i) = szero - end do - - omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) - - if (filter_mat) then - ! - ! Build the filtered matrix Af from A - ! - call la%cscnv(acsrf,info,dupl=psb_dupl_add_) - - do i=1,nrow - tmp = szero - jd = -1 - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) jd = j - if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then - tmp=tmp+acsrf%val(j) - acsrf%val(j)=szero - endif - enddo - if (jd == -1) then - write(0,*) 'Wrong input: we need the diagonal!!!!', i - else - acsrf%val(jd)=acsrf%val(jd)-tmp - end if - enddo - ! Take out zeroed terms - call acsrf%clean_zeros(info) - - ! - ! Build the smoothed prolongator using the filtered matrix - ! - do i=1,acsrf%get_nrows() - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) then - acsrf%val(j) = sone - omf(i)*acsrf%val(j) - else - acsrf%val(j) = - omf(i)*acsrf%val(j) - end if - end do - end do - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - - call af%mv_from(acsrf) - ! - ! op_prol = (I-w*D*Af)Ptilde - ! Doing it this way means to consider diag(Af_i) - ! - ! - call psb_spspmm(af,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 1' - else - ! - ! Build the smoothed prolongator using the original matrix - ! - do i=1,acsr3%get_nrows() - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if (acsr3%ja(j) == i) then - acsr3%val(j) = sone - omf(i)*acsr3%val(j) - else - acsr3%val(j) = - omf(i)*acsr3%val(j) - end if - end do - end do - - call am3%mv_from(acsr3) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - ! - ! - ! op_prol = (I-w*D*A)Ptilde - ! - ! - call psb_spspmm(am3,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - end if - - - ! - ! Ok, let's start over with the restrictor - ! - call ptilde%transc(rtilde) - call la%cscnv(atmp,info,type='csr') - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.true.,rowscale=.true.) - nrt = am4%get_nrows() - call am4%csclip(atmp2,info,lone,nrt,lone,ncol) - call atmp2%cscnv(info,type='CSR') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) - call am4%free() - call atmp2%free() - - ! This is to compute the transpose. It ONLY works if the - ! original A has a symmetric pattern. - call atmp%transc(atmp2) - call atmp2%csclip(dat,info,lone,nrow,lone,ncol) - call dat%cscnv(info,type='csr') - call dat%scal(adinv,info) - - ! Now for the product. - call psb_spspmm(dat,ptilde,datp,info) - - call datp%clone(atmp2,info) - call psb_sphalo(atmp2,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) - if (info == psb_success_) call am4%free() - - - call psb_symbmm(dat,atmp2,datdatp,info) - call psb_numbmm(dat,atmp2,datdatp) - call atmp2%free() - - call datp%mv_to(csc_datp) - call datdatp%mv_to(csc_datdatp) - - call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) - call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - - - ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp - ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden - omp = omp/oden - ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) - ! Compute omega_int - ommx = szero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = szero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - ! Going over the columns of atmp means going over the rows - ! of A^T. Hopefully ;-) - call atmp%cp_to(acsc) - - do i=1, nrow - omf(i) = ommx - do j= acsc%icp(i),acsc%icp(i+1)-1 - if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero - if(psb_minreal(omf(i)) < szero) omf(i) = szero - end do - omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) - call psb_halo(omf,desc_a,info) - call acsc%free() - - - call atmp%mv_to(acsr1) - - do i=1,acsr1%get_nrows() - do j=acsr1%irp(i),acsr1%irp(i+1)-1 - if (acsr1%ja(j) == i) then - acsr1%val(j) = sone - acsr1%val(j)*omf(acsr1%ja(j)) - else - acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) - end if - end do - end do - call atmp%mv_from(acsr1) - - call rtilde%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call rtilde%mv_from(tmpcoo) - call rtilde%cscnv(info,type='csr') - - call psb_spspmm(rtilde,atmp,op_restr,info) - - ! - ! 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,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') - goto 9999 - end if - - ! - ! Now we have to fix this. The only rows of B that are correct - ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) - ! - call op_restr%mv_to(tmpcoo) - - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - call psb_spspmm(la,tmp_prol,am3,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 2' - - call psb_sphalo(am3,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Extend am3') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done sphalo/ rwxtd' - - call psb_spspmm(op_restr,am3,ac,info) - if (info == psb_success_) call am3%free() - if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build ac = op_restr x am3') - goto 9999 - end if +!!$ ! Get the diagonal D +!!$ adiag = a%get_diag(info) +!!$ if (info == psb_success_) & +!!$ & call psb_realloc(ncol,adiag,info) +!!$ if (info == psb_success_) & +!!$ & call psb_halo(adiag,desc_a,info) +!!$ if (info == psb_success_) call a%cp_to_l(la) +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') +!!$ goto 9999 +!!$ end if +!!$ +!!$ do i=1,size(adiag) +!!$ if (adiag(i) /= szero) then +!!$ adinv(i) = sone / adiag(i) +!!$ else +!!$ adinv(i) = sone +!!$ end if +!!$ end do +!!$ +!!$ +!!$ +!!$ ! 1. Allocate Ptilde in sparse matrix form +!!$ call op_prol%mv_to(tmpcoo) +!!$ call ptilde%mv_from(tmpcoo) +!!$ call ptilde%cscnv(info,type='csr') +!!$ +!!$ if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) +!!$ if (info == psb_success_) call la%cscnv(da,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 +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & ' Initial copies done.' +!!$ +!!$ call da%scal(adinv,info) +!!$ +!!$ call psb_spspmm(da,ptilde,dap,info) +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ call dap%clone(atmp,info) +!!$ +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ call psb_spspmm(da,atmp,dadap,info) +!!$ call atmp%free() +!!$ +!!$ ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) +!!$ ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) +!!$ call dap%mv_to(csc_dap) +!!$ call dadap%mv_to(csc_dadap) +!!$ +!!$ call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) +!!$ call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ ! !$ write(0,*) trim(name),' OMP :',omp +!!$ ! !$ write(0,*) trim(name),' ODEN:',oden +!!$ +!!$ omp = omp/oden +!!$ +!!$ ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ call am3%mv_to(acsr3) +!!$ ! Compute omega_int +!!$ ommx = szero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = szero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero +!!$ if(psb_minreal(omf(i)) < szero) omf(i) = szero +!!$ end do +!!$ +!!$ omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) +!!$ +!!$ if (filter_mat) then +!!$ ! +!!$ ! Build the filtered matrix Af from A +!!$ ! +!!$ call la%cscnv(acsrf,info,dupl=psb_dupl_add_) +!!$ +!!$ do i=1,nrow +!!$ tmp = szero +!!$ jd = -1 +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) jd = j +!!$ if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then +!!$ tmp=tmp+acsrf%val(j) +!!$ acsrf%val(j)=szero +!!$ endif +!!$ enddo +!!$ if (jd == -1) then +!!$ write(0,*) 'Wrong input: we need the diagonal!!!!', i +!!$ else +!!$ acsrf%val(jd)=acsrf%val(jd)-tmp +!!$ end if +!!$ enddo +!!$ ! Take out zeroed terms +!!$ call acsrf%clean_zeros(info) +!!$ +!!$ ! +!!$ ! Build the smoothed prolongator using the filtered matrix +!!$ ! +!!$ do i=1,acsrf%get_nrows() +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) then +!!$ acsrf%val(j) = sone - omf(i)*acsrf%val(j) +!!$ else +!!$ acsrf%val(j) = - omf(i)*acsrf%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ +!!$ call af%mv_from(acsrf) +!!$ ! +!!$ ! op_prol = (I-w*D*Af)Ptilde +!!$ ! Doing it this way means to consider diag(Af_i) +!!$ ! +!!$ ! +!!$ call psb_spspmm(af,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 1' +!!$ else +!!$ ! +!!$ ! Build the smoothed prolongator using the original matrix +!!$ ! +!!$ do i=1,acsr3%get_nrows() +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if (acsr3%ja(j) == i) then +!!$ acsr3%val(j) = sone - omf(i)*acsr3%val(j) +!!$ else +!!$ acsr3%val(j) = - omf(i)*acsr3%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ call am3%mv_from(acsr3) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ ! +!!$ ! +!!$ ! op_prol = (I-w*D*A)Ptilde +!!$ ! +!!$ ! +!!$ call psb_spspmm(am3,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ end if +!!$ +!!$ +!!$ ! +!!$ ! Ok, let's start over with the restrictor +!!$ ! +!!$ call ptilde%transc(rtilde) +!!$ call la%cscnv(atmp,info,type='csr') +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.true.,rowscale=.true.) +!!$ nrt = am4%get_nrows() +!!$ call am4%csclip(atmp2,info,lone,nrt,lone,ncol) +!!$ call atmp2%cscnv(info,type='CSR') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) +!!$ call am4%free() +!!$ call atmp2%free() +!!$ +!!$ ! This is to compute the transpose. It ONLY works if the +!!$ ! original A has a symmetric pattern. +!!$ call atmp%transc(atmp2) +!!$ call atmp2%csclip(dat,info,lone,nrow,lone,ncol) +!!$ call dat%cscnv(info,type='csr') +!!$ call dat%scal(adinv,info) +!!$ +!!$ ! Now for the product. +!!$ call psb_spspmm(dat,ptilde,datp,info) +!!$ +!!$ call datp%clone(atmp2,info) +!!$ call psb_sphalo(atmp2,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ +!!$ call psb_symbmm(dat,atmp2,datdatp,info) +!!$ call psb_numbmm(dat,atmp2,datdatp) +!!$ call atmp2%free() +!!$ +!!$ call datp%mv_to(csc_datp) +!!$ call datdatp%mv_to(csc_datdatp) +!!$ +!!$ call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) +!!$ call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ +!!$ +!!$ ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp +!!$ ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden +!!$ omp = omp/oden +!!$ ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) +!!$ ! Compute omega_int +!!$ ommx = szero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = szero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ ! Going over the columns of atmp means going over the rows +!!$ ! of A^T. Hopefully ;-) +!!$ call atmp%cp_to(acsc) +!!$ +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j= acsc%icp(i),acsc%icp(i+1)-1 +!!$ if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero +!!$ if(psb_minreal(omf(i)) < szero) omf(i) = szero +!!$ end do +!!$ omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) +!!$ call psb_halo(omf,desc_a,info) +!!$ call acsc%free() +!!$ +!!$ +!!$ call atmp%mv_to(acsr1) +!!$ +!!$ do i=1,acsr1%get_nrows() +!!$ do j=acsr1%irp(i),acsr1%irp(i+1)-1 +!!$ if (acsr1%ja(j) == i) then +!!$ acsr1%val(j) = sone - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ else +!!$ acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ end if +!!$ end do +!!$ end do +!!$ call atmp%mv_from(acsr1) +!!$ +!!$ call rtilde%mv_to(tmpcoo) +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call rtilde%mv_from(tmpcoo) +!!$ call rtilde%cscnv(info,type='csr') +!!$ +!!$ call psb_spspmm(rtilde,atmp,op_restr,info) +!!$ +!!$ ! +!!$ ! 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,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') +!!$ goto 9999 +!!$ end if +!!$ +!!$ ! +!!$ ! Now we have to fix this. The only rows of B that are correct +!!$ ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) +!!$ ! +!!$ call op_restr%mv_to(tmpcoo) +!!$ +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call op_restr%mv_from(tmpcoo) +!!$ call op_restr%cscnv(info,type='csr') +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(la,tmp_prol,am3,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 2' +!!$ +!!$ call psb_sphalo(am3,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(op_restr,am3,ac,info) +!!$ if (info == psb_success_) call am3%free() +!!$ if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ &a_err='Build ac = op_restr x am3') +!!$ goto 9999 +!!$ end if diff --git a/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 index 30532e9c..d96176b2 100644 --- a/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90 @@ -116,7 +116,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_sml_parms), intent(inout) :: parms - type(psb_sspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_sspmat_type), intent(inout) :: op_prol,ac,op_restr type(psb_lsspmat_type), intent(inout) :: t_prol type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_z_dec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_z_dec_aggregator_tprol.f90 index 5a90fcd5..a64e3ebb 100644 --- a/amgprec/impl/aggregator/amg_z_dec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_z_dec_aggregator_tprol.f90 @@ -83,8 +83,8 @@ subroutine amg_z_dec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_z_dec_aggregator_type), target, intent(inout) :: ag type(amg_dml_parms), intent(inout) :: parms type(amg_daggr_data), intent(in) :: ag_data - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_zspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lzspmat_type), intent(out) :: t_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 b/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 index 84de6849..dd4ac4be 100644 --- a/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 +++ b/amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90 @@ -86,8 +86,8 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,& class(amg_z_symdec_aggregator_type), target, intent(inout) :: ag type(amg_dml_parms), intent(inout) :: parms type(amg_daggr_data), intent(in) :: ag_data - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_zspmat_type), intent(inout) :: a + type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:) type(psb_lzspmat_type), intent(out) :: op_prol integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 b/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 index c8a6e227..89d6d89d 100644 --- a/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 +++ b/amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90 @@ -105,7 +105,7 @@ ! ! subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& - & ac,desc_ac,op_prol,op_restr,info) + & ac,desc_ac,op_prol,op_restr,t_prol,info) use psb_base_mod use amg_base_prec_type use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_minnrg_bld @@ -117,8 +117,8 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms - type(psb_lzspmat_type), intent(inout) :: op_prol - type(psb_lzspmat_type), intent(out) :: ac,op_restr + type(psb_lzspmat_type), intent(inout) :: t_prol + type(psb_zspmat_type), intent(inout) :: op_prol, ac,op_restr type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info @@ -171,6 +171,8 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& filter_mat = (parms%aggr_filter == amg_filter_mat_) + !NEEDS TO BE REWORKED !! + ! naggr: number of local aggregates ! nrow: local rows. ! @@ -183,361 +185,361 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,& goto 9999 end if - ! Get the diagonal D - adiag = a%get_diag(info) - if (info == psb_success_) & - & call psb_realloc(ncol,adiag,info) - if (info == psb_success_) & - & call psb_halo(adiag,desc_a,info) - if (info == psb_success_) call a%cp_to_l(la) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') - goto 9999 - end if - - do i=1,size(adiag) - if (adiag(i) /= zzero) then - adinv(i) = zone / adiag(i) - else - adinv(i) = zone - end if - end do - - - - ! 1. Allocate Ptilde in sparse matrix form - call op_prol%mv_to(tmpcoo) - call ptilde%mv_from(tmpcoo) - call ptilde%cscnv(info,type='csr') - - if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) - if (info == psb_success_) call la%cscnv(da,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 - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies done.' - - call da%scal(adinv,info) - - call psb_spspmm(da,ptilde,dap,info) - - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') - goto 9999 - end if - - call dap%clone(atmp,info) - - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) - if (info == psb_success_) call am4%free() - - call psb_spspmm(da,atmp,dadap,info) - call atmp%free() - - ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) - ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) - call dap%mv_to(csc_dap) - call dadap%mv_to(csc_dadap) - - call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) - call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - ! !$ write(0,*) trim(name),' OMP :',omp - ! !$ write(0,*) trim(name),' ODEN:',oden - - omp = omp/oden - - ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - call am3%mv_to(acsr3) - ! Compute omega_int - ommx = zzero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = zzero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - do i=1, nrow - omf(i) = ommx - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero - if(psb_minreal(omf(i)) < dzero) omf(i) = zzero - end do - - omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) - - if (filter_mat) then - ! - ! Build the filtered matrix Af from A - ! - call la%cscnv(acsrf,info,dupl=psb_dupl_add_) - - do i=1,nrow - tmp = zzero - jd = -1 - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) jd = j - if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then - tmp=tmp+acsrf%val(j) - acsrf%val(j)=zzero - endif - enddo - if (jd == -1) then - write(0,*) 'Wrong input: we need the diagonal!!!!', i - else - acsrf%val(jd)=acsrf%val(jd)-tmp - end if - enddo - ! Take out zeroed terms - call acsrf%clean_zeros(info) - - ! - ! Build the smoothed prolongator using the filtered matrix - ! - do i=1,acsrf%get_nrows() - do j=acsrf%irp(i),acsrf%irp(i+1)-1 - if (acsrf%ja(j) == i) then - acsrf%val(j) = zone - omf(i)*acsrf%val(j) - else - acsrf%val(j) = - omf(i)*acsrf%val(j) - end if - end do - end do - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - - call af%mv_from(acsrf) - ! - ! op_prol = (I-w*D*Af)Ptilde - ! Doing it this way means to consider diag(Af_i) - ! - ! - call psb_spspmm(af,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 1' - else - ! - ! Build the smoothed prolongator using the original matrix - ! - do i=1,acsr3%get_nrows() - do j=acsr3%irp(i),acsr3%irp(i+1)-1 - if (acsr3%ja(j) == i) then - acsr3%val(j) = zone - omf(i)*acsr3%val(j) - else - acsr3%val(j) = - omf(i)*acsr3%val(j) - end if - end do - end do - - call am3%mv_from(acsr3) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done gather, going for SYMBMM 1' - ! - ! - ! op_prol = (I-w*D*A)Ptilde - ! - ! - call psb_spspmm(am3,ptilde,op_prol,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done NUMBMM 1' - - end if - - - ! - ! Ok, let's start over with the restrictor - ! - call ptilde%transc(rtilde) - call la%cscnv(atmp,info,type='csr') - call psb_sphalo(atmp,desc_a,am4,info,& - & colcnv=.true.,rowscale=.true.) - nrt = am4%get_nrows() - call am4%csclip(atmp2,info,lone,nrt,lone,ncol) - call atmp2%cscnv(info,type='CSR') - if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) - call am4%free() - call atmp2%free() - - ! This is to compute the transpose. It ONLY works if the - ! original A has a symmetric pattern. - call atmp%transc(atmp2) - call atmp2%csclip(dat,info,lone,nrow,lone,ncol) - call dat%cscnv(info,type='csr') - call dat%scal(adinv,info) - - ! Now for the product. - call psb_spspmm(dat,ptilde,datp,info) - - call datp%clone(atmp2,info) - call psb_sphalo(atmp2,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.,outfmt='CSR ') - if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) - if (info == psb_success_) call am4%free() - - - call psb_symbmm(dat,atmp2,datdatp,info) - call psb_numbmm(dat,atmp2,datdatp) - call atmp2%free() - - call datp%mv_to(csc_datp) - call datdatp%mv_to(csc_datdatp) - - call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) - call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) - call psb_sum(ctxt,omp) - call psb_sum(ctxt,oden) - - - ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp - ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden - omp = omp/oden - ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) - ! Compute omega_int - ommx = zzero - do i=1, ncol - if (ilaggr(i) >0) then - omi(i) = omp(ilaggr(i)) - else - omi(i) = zzero - end if - if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) - end do - ! Compute omega_fine - ! Going over the columns of atmp means going over the rows - ! of A^T. Hopefully ;-) - call atmp%cp_to(acsc) - - do i=1, nrow - omf(i) = ommx - do j= acsc%icp(i),acsc%icp(i+1)-1 - if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) - end do -!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero - if(psb_minreal(omf(i)) < dzero) omf(i) = zzero - end do - omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) - call psb_halo(omf,desc_a,info) - call acsc%free() - - - call atmp%mv_to(acsr1) - - do i=1,acsr1%get_nrows() - do j=acsr1%irp(i),acsr1%irp(i+1)-1 - if (acsr1%ja(j) == i) then - acsr1%val(j) = zone - acsr1%val(j)*omf(acsr1%ja(j)) - else - acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) - end if - end do - end do - call atmp%mv_from(acsr1) - - call rtilde%mv_to(tmpcoo) - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call rtilde%mv_from(tmpcoo) - call rtilde%cscnv(info,type='csr') - - call psb_spspmm(rtilde,atmp,op_restr,info) - - ! - ! 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,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') - goto 9999 - end if - - ! - ! Now we have to fix this. The only rows of B that are correct - ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) - ! - call op_restr%mv_to(tmpcoo) - - nzl = tmpcoo%get_nzeros() - i=0 - do k=1, nzl - if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then - i = i+1 - tmpcoo%val(i) = tmpcoo%val(k) - tmpcoo%ia(i) = tmpcoo%ia(k) - tmpcoo%ja(i) = tmpcoo%ja(k) - end if - end do - call tmpcoo%set_nzeros(i) - call op_restr%mv_from(tmpcoo) - call op_restr%cscnv(info,type='csr') - - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'starting sphalo/ rwxtd' - - call psb_spspmm(la,tmp_prol,am3,info) - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done SPSPMM 2' - - call psb_sphalo(am3,desc_a,am4,info,& - & colcnv=.false.,rowscale=.true.) - if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) - if (info == psb_success_) call am4%free() - - if(info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Extend am3') - goto 9999 - end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Done sphalo/ rwxtd' - - call psb_spspmm(op_restr,am3,ac,info) - if (info == psb_success_) call am3%free() - if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - &a_err='Build ac = op_restr x am3') - goto 9999 - end if +!!$ ! Get the diagonal D +!!$ adiag = a%get_diag(info) +!!$ if (info == psb_success_) & +!!$ & call psb_realloc(ncol,adiag,info) +!!$ if (info == psb_success_) & +!!$ & call psb_halo(adiag,desc_a,info) +!!$ if (info == psb_success_) call a%cp_to_l(la) +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag') +!!$ goto 9999 +!!$ end if +!!$ +!!$ do i=1,size(adiag) +!!$ if (adiag(i) /= zzero) then +!!$ adinv(i) = zone / adiag(i) +!!$ else +!!$ adinv(i) = zone +!!$ end if +!!$ end do +!!$ +!!$ +!!$ +!!$ ! 1. Allocate Ptilde in sparse matrix form +!!$ call op_prol%mv_to(tmpcoo) +!!$ call ptilde%mv_from(tmpcoo) +!!$ call ptilde%cscnv(info,type='csr') +!!$ +!!$ if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_) +!!$ if (info == psb_success_) call la%cscnv(da,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 +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & ' Initial copies done.' +!!$ +!!$ call da%scal(adinv,info) +!!$ +!!$ call psb_spspmm(da,ptilde,dap,info) +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ call dap%clone(atmp,info) +!!$ +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ call psb_spspmm(da,atmp,dadap,info) +!!$ call atmp%free() +!!$ +!!$ ! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap) +!!$ ! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap) +!!$ call dap%mv_to(csc_dap) +!!$ call dadap%mv_to(csc_dadap) +!!$ +!!$ call csc_mat_col_prod(csc_dap,csc_dadap,omp,info) +!!$ call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ ! !$ write(0,*) trim(name),' OMP :',omp +!!$ ! !$ write(0,*) trim(name),' ODEN:',oden +!!$ +!!$ omp = omp/oden +!!$ +!!$ ! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10)) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ call am3%mv_to(acsr3) +!!$ ! Compute omega_int +!!$ ommx = zzero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = zzero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero +!!$ if(psb_minreal(omf(i)) < dzero) omf(i) = zzero +!!$ end do +!!$ +!!$ omf(1:nrow) = omf(1:nrow) * adinv(1:nrow) +!!$ +!!$ if (filter_mat) then +!!$ ! +!!$ ! Build the filtered matrix Af from A +!!$ ! +!!$ call la%cscnv(acsrf,info,dupl=psb_dupl_add_) +!!$ +!!$ do i=1,nrow +!!$ tmp = zzero +!!$ jd = -1 +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) jd = j +!!$ if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then +!!$ tmp=tmp+acsrf%val(j) +!!$ acsrf%val(j)=zzero +!!$ endif +!!$ enddo +!!$ if (jd == -1) then +!!$ write(0,*) 'Wrong input: we need the diagonal!!!!', i +!!$ else +!!$ acsrf%val(jd)=acsrf%val(jd)-tmp +!!$ end if +!!$ enddo +!!$ ! Take out zeroed terms +!!$ call acsrf%clean_zeros(info) +!!$ +!!$ ! +!!$ ! Build the smoothed prolongator using the filtered matrix +!!$ ! +!!$ do i=1,acsrf%get_nrows() +!!$ do j=acsrf%irp(i),acsrf%irp(i+1)-1 +!!$ if (acsrf%ja(j) == i) then +!!$ acsrf%val(j) = zone - omf(i)*acsrf%val(j) +!!$ else +!!$ acsrf%val(j) = - omf(i)*acsrf%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ +!!$ call af%mv_from(acsrf) +!!$ ! +!!$ ! op_prol = (I-w*D*Af)Ptilde +!!$ ! Doing it this way means to consider diag(Af_i) +!!$ ! +!!$ ! +!!$ call psb_spspmm(af,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 1' +!!$ else +!!$ ! +!!$ ! Build the smoothed prolongator using the original matrix +!!$ ! +!!$ do i=1,acsr3%get_nrows() +!!$ do j=acsr3%irp(i),acsr3%irp(i+1)-1 +!!$ if (acsr3%ja(j) == i) then +!!$ acsr3%val(j) = zone - omf(i)*acsr3%val(j) +!!$ else +!!$ acsr3%val(j) = - omf(i)*acsr3%val(j) +!!$ end if +!!$ end do +!!$ end do +!!$ +!!$ call am3%mv_from(acsr3) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done gather, going for SYMBMM 1' +!!$ ! +!!$ ! +!!$ ! op_prol = (I-w*D*A)Ptilde +!!$ ! +!!$ ! +!!$ call psb_spspmm(am3,ptilde,op_prol,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done NUMBMM 1' +!!$ +!!$ end if +!!$ +!!$ +!!$ ! +!!$ ! Ok, let's start over with the restrictor +!!$ ! +!!$ call ptilde%transc(rtilde) +!!$ call la%cscnv(atmp,info,type='csr') +!!$ call psb_sphalo(atmp,desc_a,am4,info,& +!!$ & colcnv=.true.,rowscale=.true.) +!!$ nrt = am4%get_nrows() +!!$ call am4%csclip(atmp2,info,lone,nrt,lone,ncol) +!!$ call atmp2%cscnv(info,type='CSR') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2) +!!$ call am4%free() +!!$ call atmp2%free() +!!$ +!!$ ! This is to compute the transpose. It ONLY works if the +!!$ ! original A has a symmetric pattern. +!!$ call atmp%transc(atmp2) +!!$ call atmp2%csclip(dat,info,lone,nrow,lone,ncol) +!!$ call dat%cscnv(info,type='csr') +!!$ call dat%scal(adinv,info) +!!$ +!!$ ! Now for the product. +!!$ call psb_spspmm(dat,ptilde,datp,info) +!!$ +!!$ call datp%clone(atmp2,info) +!!$ call psb_sphalo(atmp2,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.,outfmt='CSR ') +!!$ if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ +!!$ call psb_symbmm(dat,atmp2,datdatp,info) +!!$ call psb_numbmm(dat,atmp2,datdatp) +!!$ call atmp2%free() +!!$ +!!$ call datp%mv_to(csc_datp) +!!$ call datdatp%mv_to(csc_datdatp) +!!$ +!!$ call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info) +!!$ call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info) +!!$ call psb_sum(ctxt,omp) +!!$ call psb_sum(ctxt,oden) +!!$ +!!$ +!!$ ! !$ write(debug_unit,*) trim(name),' OMP_R :',omp +!!$ ! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden +!!$ omp = omp/oden +!!$ ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) +!!$ ! Compute omega_int +!!$ ommx = zzero +!!$ do i=1, ncol +!!$ if (ilaggr(i) >0) then +!!$ omi(i) = omp(ilaggr(i)) +!!$ else +!!$ omi(i) = zzero +!!$ end if +!!$ if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) +!!$ end do +!!$ ! Compute omega_fine +!!$ ! Going over the columns of atmp means going over the rows +!!$ ! of A^T. Hopefully ;-) +!!$ call atmp%cp_to(acsc) +!!$ +!!$ do i=1, nrow +!!$ omf(i) = ommx +!!$ do j= acsc%icp(i),acsc%icp(i+1)-1 +!!$ if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) +!!$ end do +!!$ ! ! if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero +!!$ if(psb_minreal(omf(i)) < dzero) omf(i) = zzero +!!$ end do +!!$ omf(1:nrow) = omf(1:nrow)*adinv(1:nrow) +!!$ call psb_halo(omf,desc_a,info) +!!$ call acsc%free() +!!$ +!!$ +!!$ call atmp%mv_to(acsr1) +!!$ +!!$ do i=1,acsr1%get_nrows() +!!$ do j=acsr1%irp(i),acsr1%irp(i+1)-1 +!!$ if (acsr1%ja(j) == i) then +!!$ acsr1%val(j) = zone - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ else +!!$ acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j)) +!!$ end if +!!$ end do +!!$ end do +!!$ call atmp%mv_from(acsr1) +!!$ +!!$ call rtilde%mv_to(tmpcoo) +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call rtilde%mv_from(tmpcoo) +!!$ call rtilde%cscnv(info,type='csr') +!!$ +!!$ call psb_spspmm(rtilde,atmp,op_restr,info) +!!$ +!!$ ! +!!$ ! 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,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') +!!$ goto 9999 +!!$ end if +!!$ +!!$ ! +!!$ ! Now we have to fix this. The only rows of B that are correct +!!$ ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:) +!!$ ! +!!$ call op_restr%mv_to(tmpcoo) +!!$ +!!$ nzl = tmpcoo%get_nzeros() +!!$ i=0 +!!$ do k=1, nzl +!!$ if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then +!!$ i = i+1 +!!$ tmpcoo%val(i) = tmpcoo%val(k) +!!$ tmpcoo%ia(i) = tmpcoo%ia(k) +!!$ tmpcoo%ja(i) = tmpcoo%ja(k) +!!$ end if +!!$ end do +!!$ call tmpcoo%set_nzeros(i) +!!$ call op_restr%mv_from(tmpcoo) +!!$ call op_restr%cscnv(info,type='csr') +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(la,tmp_prol,am3,info) +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done SPSPMM 2' +!!$ +!!$ call psb_sphalo(am3,desc_a,am4,info,& +!!$ & colcnv=.false.,rowscale=.true.) +!!$ if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) +!!$ if (info == psb_success_) call am4%free() +!!$ +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'Done sphalo/ rwxtd' +!!$ +!!$ call psb_spspmm(op_restr,am3,ac,info) +!!$ if (info == psb_success_) call am3%free() +!!$ if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ &a_err='Build ac = op_restr x am3') +!!$ goto 9999 +!!$ end if diff --git a/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 b/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 index dc289067..2f944699 100644 --- a/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 +++ b/amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90 @@ -116,7 +116,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,& type(psb_desc_type), intent(inout) :: desc_a integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:) type(amg_dml_parms), intent(inout) :: parms - type(psb_zspmat_type), intent(out) :: op_prol,ac,op_restr + type(psb_zspmat_type), intent(inout) :: op_prol,ac,op_restr type(psb_lzspmat_type), intent(inout) :: t_prol type(psb_desc_type), intent(inout) :: desc_ac integer(psb_ipk_), intent(out) :: info diff --git a/amgprec/impl/amg_cprecset.F90 b/amgprec/impl/amg_cprecset.F90 index 818d51a8..ff96f6cc 100644 --- a/amgprec/impl/amg_cprecset.F90 +++ b/amgprec/impl/amg_cprecset.F90 @@ -45,7 +45,7 @@ subroutine amg_cprecsetsm(p,val,info,ilev,ilmax,pos) implicit none ! Arguments - class(amg_cprec_type), intent(inout) :: p + class(amg_cprec_type), target, intent(inout):: p class(amg_c_base_smoother_type), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: ilev,ilmax diff --git a/amgprec/impl/amg_dprecset.F90 b/amgprec/impl/amg_dprecset.F90 index 607e038f..16928156 100644 --- a/amgprec/impl/amg_dprecset.F90 +++ b/amgprec/impl/amg_dprecset.F90 @@ -45,7 +45,7 @@ subroutine amg_dprecsetsm(p,val,info,ilev,ilmax,pos) implicit none ! Arguments - class(amg_dprec_type), intent(inout) :: p + class(amg_dprec_type), target, intent(inout):: p class(amg_d_base_smoother_type), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: ilev,ilmax diff --git a/amgprec/impl/amg_sprecset.F90 b/amgprec/impl/amg_sprecset.F90 index 1e5daed4..2b2c1c98 100644 --- a/amgprec/impl/amg_sprecset.F90 +++ b/amgprec/impl/amg_sprecset.F90 @@ -45,7 +45,7 @@ subroutine amg_sprecsetsm(p,val,info,ilev,ilmax,pos) implicit none ! Arguments - class(amg_sprec_type), intent(inout) :: p + class(amg_sprec_type), target, intent(inout):: p class(amg_s_base_smoother_type), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: ilev,ilmax diff --git a/amgprec/impl/amg_zprecset.F90 b/amgprec/impl/amg_zprecset.F90 index f5ef7fd7..0337b3c3 100644 --- a/amgprec/impl/amg_zprecset.F90 +++ b/amgprec/impl/amg_zprecset.F90 @@ -45,7 +45,7 @@ subroutine amg_zprecsetsm(p,val,info,ilev,ilmax,pos) implicit none ! Arguments - class(amg_zprec_type), intent(inout) :: p + class(amg_zprec_type), target, intent(inout):: p class(amg_z_base_smoother_type), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), optional, intent(in) :: ilev,ilmax diff --git a/amgprec/impl/solver/amg_c_base_solver_csetr.f90 b/amgprec/impl/solver/amg_c_base_solver_csetr.f90 index f5751876..b78209a4 100644 --- a/amgprec/impl/solver/amg_c_base_solver_csetr.f90 +++ b/amgprec/impl/solver/amg_c_base_solver_csetr.f90 @@ -42,7 +42,7 @@ subroutine amg_c_base_solver_csetr(sv,what,val,info,idx) Implicit None ! Arguments class(amg_c_base_solver_type), intent(inout) :: sv - integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: what real(psb_spk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx diff --git a/amgprec/impl/solver/amg_c_bwgs_solver_bld.f90 b/amgprec/impl/solver/amg_c_bwgs_solver_bld.f90 index 68b87e3a..f760c80f 100644 --- a/amgprec/impl/solver/amg_c_bwgs_solver_bld.f90 +++ b/amgprec/impl/solver/amg_c_bwgs_solver_bld.f90 @@ -44,7 +44,7 @@ subroutine amg_c_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! Arguments type(psb_cspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a + Type(psb_desc_type), Intent(inout) :: desc_a class(amg_c_bwgs_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_cspmat_type), intent(in), target, optional :: b diff --git a/amgprec/impl/solver/amg_c_krm_solver_impl.f90 b/amgprec/impl/solver/amg_c_krm_solver_impl.f90 index 12e6f039..3c93488f 100644 --- a/amgprec/impl/solver/amg_c_krm_solver_impl.f90 +++ b/amgprec/impl/solver/amg_c_krm_solver_impl.f90 @@ -77,7 +77,7 @@ ! This is the implementation file corresponding to amg_c_krm_solver_mod. ! ! -subroutine amg_c_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) +subroutine amg_c_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) use psb_base_mod use amg_c_krm_solver, amg_protect_name => amg_c_krm_solver_bld @@ -85,13 +85,14 @@ subroutine amg_c_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) Implicit None ! Arguments - type(psb_cspmat_type), intent(inout), target :: a + type(psb_cspmat_type), intent(in), target :: a Type(psb_desc_type), Intent(inout) :: desc_a class(amg_c_krm_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_cspmat_type), intent(in), target, optional :: b class(psb_c_base_sparse_mat), intent(in), optional :: amold class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota integer(psb_lpk_) :: lnr diff --git a/amgprec/impl/solver/amg_d_base_solver_csetr.f90 b/amgprec/impl/solver/amg_d_base_solver_csetr.f90 index 5c4affd9..8cae84ee 100644 --- a/amgprec/impl/solver/amg_d_base_solver_csetr.f90 +++ b/amgprec/impl/solver/amg_d_base_solver_csetr.f90 @@ -42,7 +42,7 @@ subroutine amg_d_base_solver_csetr(sv,what,val,info,idx) Implicit None ! Arguments class(amg_d_base_solver_type), intent(inout) :: sv - integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: what real(psb_dpk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx diff --git a/amgprec/impl/solver/amg_d_bwgs_solver_bld.f90 b/amgprec/impl/solver/amg_d_bwgs_solver_bld.f90 index 88b56731..859c8ebe 100644 --- a/amgprec/impl/solver/amg_d_bwgs_solver_bld.f90 +++ b/amgprec/impl/solver/amg_d_bwgs_solver_bld.f90 @@ -44,7 +44,7 @@ subroutine amg_d_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! Arguments type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a + Type(psb_desc_type), Intent(inout) :: desc_a class(amg_d_bwgs_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_dspmat_type), intent(in), target, optional :: b diff --git a/amgprec/impl/solver/amg_d_krm_solver_impl.f90 b/amgprec/impl/solver/amg_d_krm_solver_impl.f90 index 63331533..dd308157 100644 --- a/amgprec/impl/solver/amg_d_krm_solver_impl.f90 +++ b/amgprec/impl/solver/amg_d_krm_solver_impl.f90 @@ -77,7 +77,7 @@ ! This is the implementation file corresponding to amg_d_krm_solver_mod. ! ! -subroutine amg_d_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) +subroutine amg_d_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) use psb_base_mod use amg_d_krm_solver, amg_protect_name => amg_d_krm_solver_bld @@ -85,13 +85,14 @@ subroutine amg_d_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) Implicit None ! Arguments - type(psb_dspmat_type), intent(inout), target :: a + type(psb_dspmat_type), intent(in), target :: a Type(psb_desc_type), Intent(inout) :: desc_a class(amg_d_krm_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_dspmat_type), intent(in), target, optional :: b class(psb_d_base_sparse_mat), intent(in), optional :: amold class(psb_d_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota integer(psb_lpk_) :: lnr diff --git a/amgprec/impl/solver/amg_s_base_solver_csetr.f90 b/amgprec/impl/solver/amg_s_base_solver_csetr.f90 index a7dcf1a4..74f59e0e 100644 --- a/amgprec/impl/solver/amg_s_base_solver_csetr.f90 +++ b/amgprec/impl/solver/amg_s_base_solver_csetr.f90 @@ -42,7 +42,7 @@ subroutine amg_s_base_solver_csetr(sv,what,val,info,idx) Implicit None ! Arguments class(amg_s_base_solver_type), intent(inout) :: sv - integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: what real(psb_spk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx diff --git a/amgprec/impl/solver/amg_s_bwgs_solver_bld.f90 b/amgprec/impl/solver/amg_s_bwgs_solver_bld.f90 index 3df9df7b..e96e1229 100644 --- a/amgprec/impl/solver/amg_s_bwgs_solver_bld.f90 +++ b/amgprec/impl/solver/amg_s_bwgs_solver_bld.f90 @@ -44,7 +44,7 @@ subroutine amg_s_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! Arguments type(psb_sspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a + Type(psb_desc_type), Intent(inout) :: desc_a class(amg_s_bwgs_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_sspmat_type), intent(in), target, optional :: b diff --git a/amgprec/impl/solver/amg_s_krm_solver_impl.f90 b/amgprec/impl/solver/amg_s_krm_solver_impl.f90 index c51aaccb..1b0efd1b 100644 --- a/amgprec/impl/solver/amg_s_krm_solver_impl.f90 +++ b/amgprec/impl/solver/amg_s_krm_solver_impl.f90 @@ -77,7 +77,7 @@ ! This is the implementation file corresponding to amg_s_krm_solver_mod. ! ! -subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) +subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) use psb_base_mod use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_bld @@ -85,13 +85,14 @@ subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) Implicit None ! Arguments - type(psb_sspmat_type), intent(inout), target :: a + type(psb_sspmat_type), intent(in), target :: a Type(psb_desc_type), Intent(inout) :: desc_a class(amg_s_krm_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_sspmat_type), intent(in), target, optional :: b class(psb_s_base_sparse_mat), intent(in), optional :: amold class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota integer(psb_lpk_) :: lnr diff --git a/amgprec/impl/solver/amg_z_base_solver_csetr.f90 b/amgprec/impl/solver/amg_z_base_solver_csetr.f90 index 80875820..5c14d396 100644 --- a/amgprec/impl/solver/amg_z_base_solver_csetr.f90 +++ b/amgprec/impl/solver/amg_z_base_solver_csetr.f90 @@ -42,7 +42,7 @@ subroutine amg_z_base_solver_csetr(sv,what,val,info,idx) Implicit None ! Arguments class(amg_z_base_solver_type), intent(inout) :: sv - integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: what real(psb_dpk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: idx diff --git a/amgprec/impl/solver/amg_z_bwgs_solver_bld.f90 b/amgprec/impl/solver/amg_z_bwgs_solver_bld.f90 index 3db433fe..dec629f5 100644 --- a/amgprec/impl/solver/amg_z_bwgs_solver_bld.f90 +++ b/amgprec/impl/solver/amg_z_bwgs_solver_bld.f90 @@ -44,7 +44,7 @@ subroutine amg_z_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) ! Arguments type(psb_zspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a + Type(psb_desc_type), Intent(inout) :: desc_a class(amg_z_bwgs_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_zspmat_type), intent(in), target, optional :: b diff --git a/amgprec/impl/solver/amg_z_krm_solver_impl.f90 b/amgprec/impl/solver/amg_z_krm_solver_impl.f90 index e4c4e308..33972c4b 100644 --- a/amgprec/impl/solver/amg_z_krm_solver_impl.f90 +++ b/amgprec/impl/solver/amg_z_krm_solver_impl.f90 @@ -77,7 +77,7 @@ ! This is the implementation file corresponding to amg_z_krm_solver_mod. ! ! -subroutine amg_z_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) +subroutine amg_z_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold) use psb_base_mod use amg_z_krm_solver, amg_protect_name => amg_z_krm_solver_bld @@ -85,13 +85,14 @@ subroutine amg_z_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold) Implicit None ! Arguments - type(psb_zspmat_type), intent(inout), target :: a + type(psb_zspmat_type), intent(in), target :: a Type(psb_desc_type), Intent(inout) :: desc_a class(amg_z_krm_solver_type), intent(inout) :: sv integer(psb_ipk_), intent(out) :: info type(psb_zspmat_type), intent(in), target, optional :: b class(psb_z_base_sparse_mat), intent(in), optional :: amold class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota integer(psb_lpk_) :: lnr