diff --git a/mlprec/impl/mld_caggrmat_asb.f90 b/mlprec/impl/mld_caggrmat_asb.f90 index 4b9637f7..28c67e8d 100644 --- a/mlprec/impl/mld_caggrmat_asb.f90 +++ b/mlprec/impl/mld_caggrmat_asb.f90 @@ -113,7 +113,7 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer :: ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act character(len=20) :: name name='mld_aggrmat_asb' @@ -122,7 +122,6 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_erractionsave(err_act) ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() call psb_info(ictxt, me, np) diff --git a/mlprec/impl/mld_caggrmat_minnrg_asb.F90 b/mlprec/impl/mld_caggrmat_minnrg_asb.F90 index d4e84297..2cfeb59c 100644 --- a/mlprec/impl/mld_caggrmat_minnrg_asb.F90 +++ b/mlprec/impl/mld_caggrmat_minnrg_asb.F90 @@ -59,7 +59,7 @@ ! of A, and omega is a suitable smoothing parameter. An estimate of the spectral ! radius of D^(-1)A, to be used in the computation of omega, is provided, ! according to the value of p%parms%aggr_omega_alg, specified by the user -! through mld_cprecinit and mld_zprecset. +! through mld_cprecinit and mld_cprecset. ! ! This routine can also build A_C according to a "bizarre" aggregation algorithm, ! using a "naive" prolongator proposed by the authors of MLD2P4. However, this @@ -119,7 +119,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! Local variables type(psb_cspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt integer :: ictxt,np,me, err_act, icomm @@ -133,6 +133,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) complex(psb_spk_), allocatable :: adiag(:), adinv(:) complex(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:) logical :: filter_mat + integer(psb_ipk_) :: ierr(5) integer :: debug_level, debug_unit integer, parameter :: ncmax=16 real(psb_spk_) :: anorm, theta @@ -162,9 +163,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -188,9 +188,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),& - & a_err='real(psb_spk_)') + info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr; + call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_spk_)') goto 9999 end if @@ -236,7 +235,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies cone.' + & ' Initial copies done.' call da%scal(adinv,info) @@ -292,7 +291,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) 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))) .lt. dzero) omf(i) = zzero +!!$ 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) @@ -411,7 +411,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! ! Ok, let's start over with the restrictor ! - call ptilde%transp(rtilde) + call ptilde%transc(rtilde) call a%cscnv(atmp,info,type='csr') call psb_sphalo(atmp,desc_a,am4,info,& & colcnv=.true.,rowscale=.true.) @@ -424,7 +424,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! This is to compute the transpose. It ONLY works if the ! original A has a symmetric pattern. - call atmp%transp(atmp2) + call atmp%transc(atmp2) call atmp2%csclip(dat,info,1,nrow,1,ncol) call dat%cscnv(info,type='csr') call dat%scal(adinv,info) @@ -473,7 +473,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) 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))) .lt. dzero) omf(i) = zzero +!!$ 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) @@ -671,14 +672,14 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) enddo ndx = nzbr(me+1) - call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,tmpcoo%val,nzbr,idisp,& - & mpi_double_complex,icomm,info) + call mpi_allgatherv(bcoo%val,ndx,mpi_complex,tmpcoo%val,nzbr,idisp,& + & mpi_complex,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/impl/mld_caggrmat_nosmth_asb.F90 b/mlprec/impl/mld_caggrmat_nosmth_asb.F90 index b3278dbe..e96b1d55 100644 --- a/mlprec/impl/mld_caggrmat_nosmth_asb.F90 +++ b/mlprec/impl/mld_caggrmat_nosmth_asb.F90 @@ -101,13 +101,15 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer ::ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act + integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name type(psb_cspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) type(psb_cspmat_type) :: am1,am2 type(psb_c_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo - integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & & naggr, nzt, naggrm1, i name='mld_aggrmat_nosmth_asb' @@ -128,9 +130,8 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) ntaggr = sum(nlaggr) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -201,12 +202,12 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) enddo ndx = nzbr(me+1) - call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,ac_coo%val,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) - call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + call mpi_allgatherv(bcoo%val,ndx,mpi_complex,ac_coo%val,nzbr,idisp,& + & mpi_complex,icomm,minfo) + call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) + call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) if(info /= psb_success_) then info=-1 call psb_errpush(info,name) diff --git a/mlprec/impl/mld_caggrmat_smth_asb.F90 b/mlprec/impl/mld_caggrmat_smth_asb.F90 index d51a0886..7ec9e44a 100644 --- a/mlprec/impl/mld_caggrmat_smth_asb.F90 +++ b/mlprec/impl/mld_caggrmat_smth_asb.F90 @@ -122,12 +122,13 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw - integer ::ictxt,np,me, err_act, icomm + integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_cspmat_type) :: am1,am2, am3, am4 type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_spk_), allocatable :: adiag(:) + integer(psb_ipk_) :: ierr(5) logical :: ml_global_nmb, filter_mat integer :: debug_level, debug_unit integer, parameter :: ncmax=16 @@ -141,7 +142,6 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) debug_level = psb_get_debug_level() ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) @@ -157,9 +157,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -187,9 +186,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(adiag(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& - & a_err='complex(psb_spk_)') + info=psb_err_alloc_request_; ierr(1)=nrow; + call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_spk_)') goto 9999 end if @@ -519,7 +517,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') goto 9999 end if if (debug_level >= psb_debug_outer_) & diff --git a/mlprec/impl/mld_daggrmat_asb.f90 b/mlprec/impl/mld_daggrmat_asb.f90 index 03208ac1..53df4f98 100644 --- a/mlprec/impl/mld_daggrmat_asb.f90 +++ b/mlprec/impl/mld_daggrmat_asb.f90 @@ -52,7 +52,7 @@ ! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine. ! The prolongator P_C is built here from this mapping, according to the ! value of p%iprcparm(mld_aggr_kind_), specified by the user through -! mld_dprecinit and mld_dprecset. +! mld_dprecinit and mld_zprecset. ! ! Currently three different prolongators are implemented, corresponding to ! three aggregation algorithms: @@ -78,7 +78,7 @@ ! ! ! Arguments: -! a - type(psb_dspmat_type), input. +! a - type(psb_cspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. ! desc_a - type(psb_desc_type), input. @@ -113,7 +113,7 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer :: ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act character(len=20) :: name name='mld_aggrmat_asb' @@ -122,7 +122,6 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_erractionsave(err_act) ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() call psb_info(ictxt, me, np) diff --git a/mlprec/impl/mld_daggrmat_minnrg_asb.F90 b/mlprec/impl/mld_daggrmat_minnrg_asb.F90 index 6fed8b3a..4b04e3bb 100644 --- a/mlprec/impl/mld_daggrmat_minnrg_asb.F90 +++ b/mlprec/impl/mld_daggrmat_minnrg_asb.F90 @@ -68,7 +68,7 @@ ! ! The coarse-level matrix A_C is distributed among the parallel processes or ! replicated on each of them, according to the value of p%parms%coarse_mat, -! specified by the user through mld_dprecinit and mld_dprecset. +! specified by the user through mld_dprecinit and mld_zprecset. ! ! For more details see ! M. Brezina and P. Vanek, A black-box iterative solver based on a @@ -119,23 +119,25 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! Local variables type(psb_dspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt integer :: ictxt,np,me, err_act, icomm character(len=20) :: name - type(psb_dspmat_type) :: am1, am2, af, ptilde, rtilde, atran, atp, atdatp - type(psb_dspmat_type) :: am3, am4, ap, adap, atmp, rada, ra, atmp2, dap, dadap, da + type(psb_dspmat_type) :: am1,am2, af, ptilde, rtilde, atran, atp, atdatp + type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_dspmat_type) :: dat, datp, datdatp, atmp3 type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo - type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf + type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc - real(psb_dpk_), allocatable :: adiag(:), omf(:),omp(:),omi(:),& - & oden(:), adinv(:) - logical :: filter_mat - integer :: debug_level, debug_unit - integer, parameter :: ncmax=16 - real(psb_dpk_) :: omega, anorm, tmp, dg, theta, alpha,beta, ommx + real(psb_dpk_), allocatable :: adiag(:), adinv(:) + real(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:) + logical :: filter_mat + integer(psb_ipk_) :: ierr(5) + integer :: debug_level, debug_unit + integer, parameter :: ncmax=16 + real(psb_dpk_) :: anorm, theta + real(psb_dpk_) :: tmp, alpha, beta, ommx name='mld_aggrmat_minnrg' if(psb_get_errstatus().ne.0) return @@ -161,9 +163,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -187,9 +188,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),& - & a_err='real(psb_dpk_)') + info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr; + call psb_errpush(info,name,i_err=ierr,a_err='real(psb_dpk_)') goto 9999 end if @@ -249,14 +249,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call dap%clone(atmp,info) - if (info == psb_success_) call psb_sphalo(atmp,desc_a,am4,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() - if(info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='sphalo 2') - goto 9999 - end if call psb_symbmm(da,atmp,dadap,info) call psb_numbmm(da,atmp,dadap) @@ -284,18 +280,19 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call am3%mv_to(acsr3) ! Compute omega_int - ommx = -huge(done) + ommx = cmplx(szero,szero) do i=1, ncol omi(i) = omp(ilaggr(i)) - ommx = max(ommx,omi(i)) + 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 - omf(i) = min(omf(i),omi(acsr3%ja(j))) + if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) end do - omf(i) = max(dzero,omf(i)) +!!$ 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) @@ -414,7 +411,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! ! Ok, let's start over with the restrictor ! - call ptilde%transp(rtilde) + call ptilde%transc(rtilde) call a%cscnv(atmp,info,type='csr') call psb_sphalo(atmp,desc_a,am4,info,& & colcnv=.true.,rowscale=.true.) @@ -427,7 +424,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! This is to compute the transpose. It ONLY works if the ! original A has a symmetric pattern. - call atmp%transp(atmp2) + call atmp%transc(atmp2) call atmp2%csclip(dat,info,1,nrow,1,ncol) call dat%cscnv(info,type='csr') call dat%scal(adinv,info) @@ -461,10 +458,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = -huge(done) + ommx = cmplx(szero,szero) do i=1, ncol omi(i) = omp(ilaggr(i)) - ommx = max(ommx,omi(i)) + 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 @@ -474,9 +471,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) do i=1, nrow omf(i) = ommx do j= acsc%icp(i),acsc%icp(i+1)-1 - omf(i) = min(omf(i),omi(acsc%ia(j))) + if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) end do - omf(i) = max(dzero,omf(i)) +!!$ 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) @@ -629,14 +627,14 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call p%ac%cscnv(info,type='csr') if (np>1) then - call am1%mv_to(acsr1) - nzl = acsr1%get_nzeros() - call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') + call am1%mv_to(acsr) + nzl = acsr%get_nzeros() + call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I') if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') goto 9999 end if - call am1%mv_from(acsr1) + call am1%mv_from(acsr) endif call am1%set_ncols(p%desc_ac%get_local_cols()) @@ -677,11 +675,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,tmpcoo%val,nzbr,idisp,& & mpi_double_precision,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -817,7 +815,7 @@ contains if (ip1 > nv1) exit if (ip2 > nv2) exit if (iv1(ip1) == iv2(ip2)) then - dot = dot + v1(ip1)*v2(ip2) + dot = dot + (v1(ip1))*v2(ip2) ip1 = ip1 + 1 ip2 = ip2 + 1 else if (iv1(ip1) < iv2(ip2)) then diff --git a/mlprec/impl/mld_daggrmat_nosmth_asb.F90 b/mlprec/impl/mld_daggrmat_nosmth_asb.F90 index 6919f164..715f3601 100644 --- a/mlprec/impl/mld_daggrmat_nosmth_asb.F90 +++ b/mlprec/impl/mld_daggrmat_nosmth_asb.F90 @@ -51,7 +51,7 @@ ! ! The coarse-level matrix A_C is distributed among the parallel processes or ! replicated on each of them, according to the value of p%parms%coarse_mat -! specified by the user through mld_dprecinit and mld_dprecset. +! specified by the user through mld_dprecinit and mld_zprecset. ! ! For details see ! P. D'Ambra, D. di Serafino and S. Filippone, On the development of @@ -101,13 +101,15 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer ::ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act + integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name type(psb_dspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) type(psb_dspmat_type) :: am1,am2 type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo - integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & & naggr, nzt, naggrm1, i name='mld_aggrmat_nosmth_asb' @@ -128,9 +130,8 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) ntaggr = sum(nlaggr) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -202,11 +203,11 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) ndx = nzbr(me+1) call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,ac_coo%val,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) - call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & mpi_double_precision,icomm,minfo) + call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) + call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) if(info /= psb_success_) then info=-1 call psb_errpush(info,name) diff --git a/mlprec/impl/mld_daggrmat_smth_asb.F90 b/mlprec/impl/mld_daggrmat_smth_asb.F90 index cf905100..fb917201 100644 --- a/mlprec/impl/mld_daggrmat_smth_asb.F90 +++ b/mlprec/impl/mld_daggrmat_smth_asb.F90 @@ -59,7 +59,7 @@ ! of A, and omega is a suitable smoothing parameter. An estimate of the spectral ! radius of D^(-1)A, to be used in the computation of omega, is provided, ! according to the value of p%parms%aggr_omega_alg, specified by the user -! through mld_dprecinit and mld_dprecset. +! through mld_dprecinit and mld_zprecset. ! ! This routine can also build A_C according to a "bizarre" aggregation algorithm, ! using a "naive" prolongator proposed by the authors of MLD2P4. However, this @@ -68,7 +68,7 @@ ! ! The coarse-level matrix A_C is distributed among the parallel processes or ! replicated on each of them, according to the value of p%parms%coarse_mat, -! specified by the user through mld_dprecinit and mld_dprecset. +! specified by the user through mld_dprecinit and mld_zprecset. ! ! For more details see ! M. Brezina and P. Vanek, A black-box iterative solver based on a @@ -122,16 +122,17 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw - integer ::ictxt,np,me, err_act, icomm + integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_dspmat_type) :: am1,am2, am3, am4 type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_dpk_), allocatable :: adiag(:) + integer(psb_ipk_) :: ierr(5) logical :: ml_global_nmb, filter_mat integer :: debug_level, debug_unit integer, parameter :: ncmax=16 - real(psb_dpk_) :: omega, anorm, tmp, dg, theta + real(psb_dpk_) :: anorm, omega, tmp, dg, theta name='mld_aggrmat_smth_asb' if(psb_get_errstatus().ne.0) return @@ -141,7 +142,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) debug_level = psb_get_debug_level() ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) @@ -157,9 +157,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -187,9 +186,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(adiag(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& - & a_err='real(psb_dpk_)') + info=psb_err_alloc_request_; ierr(1)=nrow; + call psb_errpush(info,name,i_err=ierr,a_err='real(psb_dpk_)') goto 9999 end if @@ -228,7 +226,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies done.' + & ' Initial copies sone.' if (filter_mat) then ! @@ -296,7 +294,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) dg = done nrw = acsr3%get_nrows() do i=1, nrw - tmp = dzero + tmp = szero do j=acsr3%irp(i),acsr3%irp(i+1)-1 if (acsr3%ja(j) <= nrw) then tmp = tmp + abs(acsr3%val(j)) @@ -519,7 +517,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') goto 9999 end if if (debug_level >= psb_debug_outer_) & diff --git a/mlprec/impl/mld_saggrmat_asb.f90 b/mlprec/impl/mld_saggrmat_asb.f90 index 9045e32f..526d5cbd 100644 --- a/mlprec/impl/mld_saggrmat_asb.f90 +++ b/mlprec/impl/mld_saggrmat_asb.f90 @@ -52,7 +52,7 @@ ! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine. ! The prolongator P_C is built here from this mapping, according to the ! value of p%iprcparm(mld_aggr_kind_), specified by the user through -! mld_sprecinit and mld_dprecset. +! mld_sprecinit and mld_zprecset. ! ! Currently three different prolongators are implemented, corresponding to ! three aggregation algorithms: @@ -78,7 +78,7 @@ ! ! ! Arguments: -! a - type(psb_sspmat_type), input. +! a - type(psb_cspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. ! desc_a - type(psb_desc_type), input. @@ -113,7 +113,7 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer :: ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act character(len=20) :: name name='mld_aggrmat_asb' @@ -122,7 +122,6 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_erractionsave(err_act) ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() call psb_info(ictxt, me, np) diff --git a/mlprec/impl/mld_saggrmat_minnrg_asb.F90 b/mlprec/impl/mld_saggrmat_minnrg_asb.F90 index af31b57e..d3829933 100644 --- a/mlprec/impl/mld_saggrmat_minnrg_asb.F90 +++ b/mlprec/impl/mld_saggrmat_minnrg_asb.F90 @@ -59,7 +59,7 @@ ! of A, and omega is a suitable smoothing parameter. An estimate of the spectral ! radius of D^(-1)A, to be used in the computation of omega, is provided, ! according to the value of p%parms%aggr_omega_alg, specified by the user -! through mld_sprecinit and mld_dprecset. +! through mld_sprecinit and mld_sprecset. ! ! This routine can also build A_C according to a "bizarre" aggregation algorithm, ! using a "naive" prolongator proposed by the authors of MLD2P4. However, this @@ -68,7 +68,7 @@ ! ! The coarse-level matrix A_C is distributed among the parallel processes or ! replicated on each of them, according to the value of p%parms%coarse_mat, -! specified by the user through mld_sprecinit and mld_dprecset. +! specified by the user through mld_sprecinit and mld_zprecset. ! ! For more details see ! M. Brezina and P. Vanek, A black-box iterative solver based on a @@ -119,7 +119,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! Local variables type(psb_sspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt integer :: ictxt,np,me, err_act, icomm @@ -130,12 +130,14 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc - real(psb_spk_), allocatable :: adiag(:), omf(:),omp(:),omi(:),& - & oden(:), adinv(:) - logical :: filter_mat - integer :: debug_level, debug_unit - integer, parameter :: ncmax=16 - real(psb_spk_) :: omega, anorm, tmp, dg, theta, alpha,beta, ommx + real(psb_spk_), allocatable :: adiag(:), adinv(:) + real(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:) + logical :: filter_mat + integer(psb_ipk_) :: ierr(5) + integer :: debug_level, debug_unit + integer, parameter :: ncmax=16 + real(psb_spk_) :: anorm, theta + real(psb_spk_) :: tmp, alpha, beta, ommx name='mld_aggrmat_minnrg' if(psb_get_errstatus().ne.0) return @@ -161,9 +163,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -187,9 +188,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),& - & a_err='real(psb_spk_)') + info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr; + call psb_errpush(info,name,i_err=ierr,a_err='real(psb_spk_)') goto 9999 end if @@ -235,7 +235,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies sone.' + & ' Initial copies done.' call da%scal(adinv,info) @@ -280,18 +280,19 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call am3%mv_to(acsr3) ! Compute omega_int - ommx = -huge(sone) + ommx = cmplx(szero,szero) do i=1, ncol omi(i) = omp(ilaggr(i)) - ommx = max(ommx,omi(i)) + 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 - omf(i) = min(omf(i),omi(acsr3%ja(j))) + if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j)) end do - omf(i) = max(szero,omf(i)) +!!$ 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) @@ -410,7 +411,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! ! Ok, let's start over with the restrictor ! - call ptilde%transp(rtilde) + call ptilde%transc(rtilde) call a%cscnv(atmp,info,type='csr') call psb_sphalo(atmp,desc_a,am4,info,& & colcnv=.true.,rowscale=.true.) @@ -423,7 +424,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! This is to compute the transpose. It ONLY works if the ! original A has a symmetric pattern. - call atmp%transp(atmp2) + call atmp%transc(atmp2) call atmp2%csclip(dat,info,1,nrow,1,ncol) call dat%cscnv(info,type='csr') call dat%scal(adinv,info) @@ -457,10 +458,10 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = -huge(sone) + ommx = cmplx(szero,szero) do i=1, ncol omi(i) = omp(ilaggr(i)) - ommx = max(ommx,omi(i)) + 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 @@ -470,9 +471,10 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) do i=1, nrow omf(i) = ommx do j= acsc%icp(i),acsc%icp(i+1)-1 - omf(i) = min(omf(i),omi(acsc%ia(j))) + if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j)) end do - omf(i) = max(szero,omf(i)) +!!$ 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) @@ -670,14 +672,14 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) enddo ndx = nzbr(me+1) - call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,tmpcoo%val,nzbr,idisp,& - & mpi_double_precision,icomm,info) + call mpi_allgatherv(bcoo%val,ndx,mpi_real,tmpcoo%val,nzbr,idisp,& + & mpi_real,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& @@ -813,7 +815,7 @@ contains if (ip1 > nv1) exit if (ip2 > nv2) exit if (iv1(ip1) == iv2(ip2)) then - dot = dot + v1(ip1)*v2(ip2) + dot = dot + (v1(ip1))*v2(ip2) ip1 = ip1 + 1 ip2 = ip2 + 1 else if (iv1(ip1) < iv2(ip2)) then diff --git a/mlprec/impl/mld_saggrmat_nosmth_asb.F90 b/mlprec/impl/mld_saggrmat_nosmth_asb.F90 index 03d1ae85..86995791 100644 --- a/mlprec/impl/mld_saggrmat_nosmth_asb.F90 +++ b/mlprec/impl/mld_saggrmat_nosmth_asb.F90 @@ -51,7 +51,7 @@ ! ! The coarse-level matrix A_C is distributed among the parallel processes or ! replicated on each of them, according to the value of p%parms%coarse_mat -! specified by the user through mld_sprecinit and mld_dprecset. +! specified by the user through mld_sprecinit and mld_zprecset. ! ! For details see ! P. D'Ambra, D. di Serafino and S. Filippone, On the development of @@ -101,13 +101,15 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer ::ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act + integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name type(psb_sspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) type(psb_sspmat_type) :: am1,am2 type(psb_s_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo - integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & & naggr, nzt, naggrm1, i name='mld_aggrmat_nosmth_asb' @@ -128,9 +130,8 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) ntaggr = sum(nlaggr) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -201,12 +202,12 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) enddo ndx = nzbr(me+1) - call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,ac_coo%val,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) - call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + call mpi_allgatherv(bcoo%val,ndx,mpi_real,ac_coo%val,nzbr,idisp,& + & mpi_real,icomm,minfo) + call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) + call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) if(info /= psb_success_) then info=-1 call psb_errpush(info,name) diff --git a/mlprec/impl/mld_saggrmat_smth_asb.F90 b/mlprec/impl/mld_saggrmat_smth_asb.F90 index 04a87cc1..47519967 100644 --- a/mlprec/impl/mld_saggrmat_smth_asb.F90 +++ b/mlprec/impl/mld_saggrmat_smth_asb.F90 @@ -59,7 +59,7 @@ ! of A, and omega is a suitable smoothing parameter. An estimate of the spectral ! radius of D^(-1)A, to be used in the computation of omega, is provided, ! according to the value of p%parms%aggr_omega_alg, specified by the user -! through mld_sprecinit and mld_dprecset. +! through mld_sprecinit and mld_zprecset. ! ! This routine can also build A_C according to a "bizarre" aggregation algorithm, ! using a "naive" prolongator proposed by the authors of MLD2P4. However, this @@ -68,7 +68,7 @@ ! ! The coarse-level matrix A_C is distributed among the parallel processes or ! replicated on each of them, according to the value of p%parms%coarse_mat, -! specified by the user through mld_sprecinit and mld_dprecset. +! specified by the user through mld_sprecinit and mld_zprecset. ! ! For more details see ! M. Brezina and P. Vanek, A black-box iterative solver based on a @@ -122,16 +122,17 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw - integer ::ictxt,np,me, err_act, icomm + integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_sspmat_type) :: am1,am2, am3, am4 type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde real(psb_spk_), allocatable :: adiag(:) + integer(psb_ipk_) :: ierr(5) logical :: ml_global_nmb, filter_mat integer :: debug_level, debug_unit integer, parameter :: ncmax=16 - real(psb_spk_) :: omega, anorm, tmp, dg, theta + real(psb_spk_) :: anorm, omega, tmp, dg, theta name='mld_aggrmat_smth_asb' if(psb_get_errstatus().ne.0) return @@ -141,7 +142,6 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) debug_level = psb_get_debug_level() ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) @@ -157,9 +157,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -187,9 +186,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(adiag(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& - & a_err='real(psb_spk_)') + info=psb_err_alloc_request_; ierr(1)=nrow; + call psb_errpush(info,name,i_err=ierr,a_err='real(psb_spk_)') goto 9999 end if @@ -519,7 +517,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') goto 9999 end if if (debug_level >= psb_debug_outer_) & diff --git a/mlprec/impl/mld_zaggrmat_asb.f90 b/mlprec/impl/mld_zaggrmat_asb.f90 index 11d302a6..91f14443 100644 --- a/mlprec/impl/mld_zaggrmat_asb.f90 +++ b/mlprec/impl/mld_zaggrmat_asb.f90 @@ -78,7 +78,7 @@ ! ! ! Arguments: -! a - type(psb_zspmat_type), input. +! a - type(psb_cspmat_type), input. ! The sparse matrix structure containing the local part of ! the fine-level matrix. ! desc_a - type(psb_desc_type), input. @@ -113,7 +113,7 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer :: ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act character(len=20) :: name name='mld_aggrmat_asb' @@ -122,7 +122,6 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) call psb_erractionsave(err_act) ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() call psb_info(ictxt, me, np) diff --git a/mlprec/impl/mld_zaggrmat_minnrg_asb.F90 b/mlprec/impl/mld_zaggrmat_minnrg_asb.F90 index 053ddbb9..fd72fe7e 100644 --- a/mlprec/impl/mld_zaggrmat_minnrg_asb.F90 +++ b/mlprec/impl/mld_zaggrmat_minnrg_asb.F90 @@ -119,7 +119,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! Local variables type(psb_zspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt integer :: ictxt,np,me, err_act, icomm @@ -133,6 +133,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) complex(psb_dpk_), allocatable :: adiag(:), adinv(:) complex(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:) logical :: filter_mat + integer(psb_ipk_) :: ierr(5) integer :: debug_level, debug_unit integer, parameter :: ncmax=16 real(psb_dpk_) :: anorm, theta @@ -162,9 +163,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -188,9 +188,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) & omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),& - & a_err='real(psb_dpk_)') + info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr; + call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_dpk_)') goto 9999 end if @@ -236,7 +235,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies zone.' + & ' Initial copies done.' call da%scal(adinv,info) @@ -281,7 +280,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call am3%mv_to(acsr3) ! Compute omega_int - ommx = cmplx(dzero,dzero) + ommx = cmplx(szero,szero) do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) @@ -292,7 +291,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) 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))) .lt. dzero) omf(i) = zzero +!!$ 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) @@ -411,7 +411,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! ! Ok, let's start over with the restrictor ! - call ptilde%transp(rtilde) + call ptilde%transc(rtilde) call a%cscnv(atmp,info,type='csr') call psb_sphalo(atmp,desc_a,am4,info,& & colcnv=.true.,rowscale=.true.) @@ -424,7 +424,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) ! This is to compute the transpose. It ONLY works if the ! original A has a symmetric pattern. - call atmp%transp(atmp2) + call atmp%transc(atmp2) call atmp2%csclip(dat,info,1,nrow,1,ncol) call dat%cscnv(info,type='csr') call dat%scal(adinv,info) @@ -458,7 +458,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) omp = omp/oden ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! Compute omega_int - ommx = cmplx(dzero,dzero) + ommx = cmplx(szero,szero) do i=1, ncol omi(i) = omp(ilaggr(i)) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) @@ -473,7 +473,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) 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))) .lt. dzero) omf(i) = zzero +!!$ 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) @@ -674,11 +675,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,tmpcoo%val,nzbr,idisp,& & mpi_double_complex,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info == psb_success_)& - & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,info) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/impl/mld_zaggrmat_nosmth_asb.F90 b/mlprec/impl/mld_zaggrmat_nosmth_asb.F90 index aebd75ad..40be6002 100644 --- a/mlprec/impl/mld_zaggrmat_nosmth_asb.F90 +++ b/mlprec/impl/mld_zaggrmat_nosmth_asb.F90 @@ -101,13 +101,15 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, intent(out) :: info ! Local variables - integer ::ictxt,np,me, err_act, icomm + integer :: ictxt,np,me, err_act + integer(psb_mpik_) :: icomm, ndx, minfo character(len=20) :: name type(psb_zspmat_type) :: b - integer, allocatable :: nzbr(:), idisp(:) + integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) + integer(psb_ipk_) :: ierr(5) type(psb_zspmat_type) :: am1,am2 type(psb_z_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo - integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& + integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & & naggr, nzt, naggrm1, i name='mld_aggrmat_nosmth_asb' @@ -128,9 +130,8 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) ntaggr = sum(nlaggr) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -202,11 +203,11 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) ndx = nzbr(me+1) call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,ac_coo%val,nzbr,idisp,& - & mpi_double_precision,icomm,info) - call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,& - & psb_mpi_integer,icomm,info) - call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,& - & psb_mpi_integer,icomm,info) + & mpi_double_complex,icomm,minfo) + call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) + call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,& + & psb_mpi_ipk_integer,icomm,minfo) if(info /= psb_success_) then info=-1 call psb_errpush(info,name) diff --git a/mlprec/impl/mld_zaggrmat_smth_asb.F90 b/mlprec/impl/mld_zaggrmat_smth_asb.F90 index 805970dc..9eb44d94 100644 --- a/mlprec/impl/mld_zaggrmat_smth_asb.F90 +++ b/mlprec/impl/mld_zaggrmat_smth_asb.F90 @@ -122,12 +122,13 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) integer, allocatable :: nzbr(:), idisp(:) integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw - integer ::ictxt,np,me, err_act, icomm + integer ::ictxt, np, me, err_act character(len=20) :: name type(psb_zspmat_type) :: am1,am2, am3, am4 type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde complex(psb_dpk_), allocatable :: adiag(:) + integer(psb_ipk_) :: ierr(5) logical :: ml_global_nmb, filter_mat integer :: debug_level, debug_unit integer, parameter :: ncmax=16 @@ -141,7 +142,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) debug_level = psb_get_debug_level() ictxt = desc_a%get_context() - icomm = desc_a%get_mpic() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) @@ -157,9 +157,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(nzbr(np), idisp(np),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),& - & a_err='integer') + info=psb_err_alloc_request_; ierr(1)=2*np; + call psb_errpush(info,name,i_err=ierr,a_err='integer') goto 9999 end if @@ -187,9 +186,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) allocate(adiag(ncol),stat=info) if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),& - & a_err='complex(psb_dpk_)') + info=psb_err_alloc_request_; ierr(1)=nrow; + call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_dpk_)') goto 9999 end if @@ -228,7 +226,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' Initial copies done.' + & ' Initial copies sone.' if (filter_mat) then ! @@ -296,7 +294,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) dg = done nrw = acsr3%get_nrows() do i=1, nrw - tmp = dzero + tmp = szero do j=acsr3%irp(i),acsr3%irp(i+1)-1 if (acsr3%ja(j) <= nrw) then tmp = tmp + abs(acsr3%val(j)) @@ -519,7 +517,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I') if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I') if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac') + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Creating p%desc_ac and converting ac') goto 9999 end if if (debug_level >= psb_debug_outer_) & diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index 46fe8d3d..40165150 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -5,16 +5,16 @@ CSR ! Storage format CSR COO JAD 0100 ! ITMAX -1 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) -1.d-9 ! EPS -3L-M-RAS-I-D4 ! Longer descriptive name for preconditioner (up to 20 chars) +1.d-6 ! EPS +3L-MUL-RAS-BJAC4-ILU ! Descriptive name for preconditioner (up to 40 chars) ML ! Preconditioner NONE JACOBI BJAC AS ML -0 ! Number of overlap layers for AS preconditioner at finest level +1 ! Number of overlap layers for AS preconditioner at finest level HALO ! Restriction operator NONE HALO NONE ! Prolongation operator NONE SUM AVG ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU -0 ! Level-set N for ILU(N) -1.d-5 ! Threshold T for ILU(T,P) -1 ! Smoother/Jacobi sweeps +0 ! Level-set N for ILU(N), and P for ILUT +1.d-4 ! Threshold T for ILU(T,P) +4 ! Smoother/Jacobi sweeps BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML 3 ! Number of levels in a multilevel preconditioner SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED