fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

mld2p4-NewML:

Browse Source 
mlprec/impl/Makefile
 mlprec/impl/mld_caggrmat_asb.f90
 mlprec/impl/mld_caggrmat_minnrg_asb.F90
 mlprec/impl/mld_caggrmat_nosmth_asb.F90
 mlprec/impl/mld_caggrmat_smth_asb.F90
 mlprec/impl/mld_daggrmat_asb.f90
 mlprec/impl/mld_daggrmat_minnrg_asb.F90
 mlprec/impl/mld_daggrmat_nosmth_asb.F90
 mlprec/impl/mld_daggrmat_smth_asb.F90
 mlprec/impl/mld_dprecinit.F90
 mlprec/impl/mld_dprecset.F90
 mlprec/impl/mld_saggrmat_asb.f90
 mlprec/impl/mld_saggrmat_minnrg_asb.F90
 mlprec/impl/mld_saggrmat_nosmth_asb.F90
 mlprec/impl/mld_saggrmat_smth_asb.F90
 mlprec/impl/mld_zaggrmat_asb.f90
 mlprec/impl/mld_zaggrmat_minnrg_asb.F90
 mlprec/impl/mld_zaggrmat_nosmth_asb.F90
 mlprec/impl/mld_zaggrmat_smth_asb.F90
 mlprec/mld_base_prec_type.F90
 mlprec/mld_c_inner_mod.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_d_prec_type.f90
 mlprec/mld_s_inner_mod.f90
 mlprec/mld_z_inner_mod.f90
 tests/pdegen/runs/ppde.inp

1. New _biz_asb routines
2. New interface declarations.
stopcriterion
Salvatore Filippone 13 years ago
parent 3203682a24
commit e5f9b851ce
26 changed files with 1249 additions and 1753 deletions
Show Stats Download Patch File Download Diff File

8
mlprec/impl/Makefile
Unescape Escape View File

@ -7,13 +7,13 @@ HERE=..
FINCLUDES=$(FMFLAG).. $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBINCDIR) $(FMFLAG)$(PSBLIBDIR) FINCLUDES=$(FMFLAG).. $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBINCDIR) $(FMFLAG)$(PSBLIBDIR)
DMPFOBJS=mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o DMPFOBJS=mld_daggrmat_nosmth_asb.o mld_daggrmat_smth_asb.o mld_daggrmat_minnrg_asb.o mld_daggrmat_biz_asb.o
SMPFOBJS=mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o mld_saggrmat_minnrg_asb.o SMPFOBJS=mld_saggrmat_nosmth_asb.o mld_saggrmat_smth_asb.o mld_saggrmat_minnrg_asb.o mld_saggrmat_biz_asb.o
ZMPFOBJS=mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o mld_zaggrmat_minnrg_asb.o ZMPFOBJS=mld_zaggrmat_nosmth_asb.o mld_zaggrmat_smth_asb.o mld_zaggrmat_minnrg_asb.o mld_zaggrmat_biz_asb.o
CMPFOBJS=mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o mld_caggrmat_minnrg_asb.o CMPFOBJS=mld_caggrmat_nosmth_asb.o mld_caggrmat_smth_asb.o mld_caggrmat_minnrg_asb.o mld_caggrmat_biz_asb.o
MPFOBJS=$(SMPFOBJS) $(DMPFOBJS) $(CMPFOBJS) $(ZMPFOBJS) MPFOBJS=$(SMPFOBJS) $(DMPFOBJS) $(CMPFOBJS) $(ZMPFOBJS)

17
mlprec/impl/mld_caggrmat_asb.f90
Unescape Escape View File

@ -113,6 +113,7 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info integer, intent(out) :: info
! Local variables ! Local variables
type(psb_cspmat_type) :: b, op_prol,op_restr
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -128,15 +129,23 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case (p%parms%aggr_kind) select case (p%parms%aggr_kind)
case (mld_no_smooth_) case (mld_no_smooth_)
call mld_aggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999 goto 9999
end if end if
case(mld_smooth_prol_,mld_biz_prol_) case(mld_smooth_prol_)
call mld_aggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_)
call mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
@ -144,7 +153,7 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_min_energy_) case(mld_min_energy_)
call mld_aggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999

152
mlprec/impl/mld_caggrmat_minnrg_asb.F90
Unescape Escape View File

@ -116,15 +116,15 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_c_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_cspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
type(psb_cspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm integer :: ictxt,np,me, err_act, icomm
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: am1,am2, af, ptilde, rtilde, atran, atp, atdatp type(psb_cspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_cspmat_type) :: dat, datp, datdatp, atmp3 type(psb_cspmat_type) :: dat, datp, datdatp, atmp3
type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
@ -354,17 +354,17 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*Af)Ptilde ! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i) ! Doing it this way means to consider diag(Af_i)
! !
! !
call psb_symbmm(af,ptilde,am1,info) call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(af,ptilde,am1) call psb_numbmm(af,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -390,16 +390,16 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*A)Ptilde ! op_prol = (I-w*D*A)Ptilde
! !
! !
call psb_symbmm(am3,ptilde,am1,info) call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(am3,ptilde,am1) call psb_numbmm(am3,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -509,20 +509,20 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call rtilde%mv_from(tmpcoo) call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr') call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,am2,info) call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,am2) call psb_numbmm(rtilde,atmp,op_restr)
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
@ -530,7 +530,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Now we have to fix this. The only rows of B that are correct ! 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(:) ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
! !
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
@ -543,21 +543,21 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
end do end do
call tmpcoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'
@ -576,14 +576,13 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd' & 'Done sphalo/ rwxtd'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
call b%mv_to(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
&a_err='Build b = am2 x am3') &a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
@ -597,6 +596,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_distr_mat_) case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros() nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -627,29 +627,29 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%cscnv(info,type='csr') call p%ac%cscnv(info,type='csr')
if (np>1) then if (np>1) then
call am1%mv_to(acsr) call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros() nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr) call op_prol%mv_from(acsr)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -659,46 +659,54 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! !
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
nzbr(:) = 0 if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
nzbr(me+1) = bcoo%get_nzeros()
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
if (info /= psb_success_) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& 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 if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
& a_err=' from mpi_allgatherv')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call bcoo%free()
call tmpcoo%fix(info)
call p%ac%mv_from(tmpcoo)
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) goto 9999 if(info /= psb_success_) goto 9999
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
deallocate(nzbr,idisp,stat=info) !!$ nzbr(:) = 0
if (info /= psb_success_) then !!$ nzbr(me+1) = bcoo%get_nzeros()
info = psb_err_alloc_dealloc_ !!$
call psb_errpush(info,name) !!$ call psb_sum(ictxt,nzbr(1:np))
goto 9999 !!$ nzac = sum(nzbr)
end if !!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
!!$ if (info /= psb_success_) goto 9999
!!$
!!$ do ip=1,np
!!$ idispip) = sum(nzbr(1:ip-1))
!!$ enddo
!!$ ndx = nzbr(me+1)
!!$
!!$ 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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
!!$ & psb_mpi_ipk_integer,icomm,info)
!!$ if (info == psb_success_)&
!!$ & 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,&
!!$ & a_err=' from mpi_allgatherv')
!!$ goto 9999
!!$ end if
!!$
!!$ call bcoo%free()
!!$ call tmpcoo%fix(info)
!!$ call p%ac%mv_from(tmpcoo)
!!$ call p%ac%cscnv(info,type='csr')
!!$ if(info /= psb_success_) goto 9999
!!$
!!$ deallocate(nzbr,idisp,stat=info)
!!$ if (info /= psb_success_) then
!!$ info = psb_err_alloc_dealloc_
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
@ -715,13 +723,13 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => R i.e. restriction operator ! op_restr => R i.e. restriction operator
! am1 => P i.e. prolongation operator ! op_prol => P i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

125
mlprec/impl/mld_caggrmat_nosmth_asb.F90
Unescape Escape View File

@ -99,17 +99,17 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_c_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_cspmat_type) :: b, op_prol,op_restr
! Local variables ! Local variables
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_cspmat_type) :: am1,am2 type(psb_c_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo
type(psb_c_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo type(psb_c_csr_sparse_mat) :: acsr1, acsr2
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, i & naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb' name='mld_aggrmat_nosmth_asb'
@ -128,32 +128,19 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1=sum(nlaggr(1:me)) naggrm1=sum(nlaggr(1:me))
if (p%parms%coarse_mat == mld_repl_mat_) then
do i=1, nrow do i=1, nrow
ilaggr(i) = ilaggr(i) + naggrm1 ilaggr(i) = ilaggr(i) + naggrm1
end do end do
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
end if
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
if (p%parms%coarse_mat == mld_repl_mat_) then
call acoo1%allocate(ncol,ntaggr,ncol) call acoo1%allocate(ncol,ntaggr,ncol)
else
call acoo1%allocate(ncol,naggr,ncol)
end if
do i=1,nrow do i=1,nrow
acoo1%val(i) = cone acoo1%val(i) = cone
@ -165,10 +152,13 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call acoo1%set_nzeros(nrow) call acoo1%set_nzeros(nrow)
call acoo1%set_asb() call acoo1%set_asb()
call acoo1%fix(info) call acoo1%fix(info)
call acoo1%transp(acoo2)
call a%csclip(bcoo,info,jmax=nrow)
call op_prol%mv_from(acoo1)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = bcoo%get_nzeros()
do i=1, nzt do i=1, nzt
@ -181,6 +171,8 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call bcoo%fix(info) call bcoo%fix(info)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then if (p%parms%coarse_mat == mld_repl_mat_) then
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
@ -189,55 +181,74 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
nzbr(:) = 0 else if (p%parms%coarse_mat == mld_distr_mat_) then
nzbr(me+1) = nzt
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call ac_coo%allocate(ntaggr,ntaggr,nzac) nzl = b%get_nzeros()
call b%mv_to(bcoo)
do ip=1,np if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
idisp(ip) = sum(nzbr(1:ip-1)) if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
enddo if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
ndx = nzbr(me+1) 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')
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 if (info /= psb_success_) then
info=-1 call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name) & a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
call ac_coo%set_nzeros(nzac) if (debug_level >= psb_debug_outer_) &
call ac_coo%set_dupl(psb_dupl_add_) & write(debug_unit,*) me,' ',trim(name),&
call ac_coo%fix(info) & 'Assembld aux descr. distr.'
call p%ac%mv_from(ac_coo) call p%ac%mv_from(bcoo)
else if (p%parms%coarse_mat == mld_distr_mat_) then call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
call p%ac%mv_from(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac, desc_ac') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
end if
call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end if end if
call bcoo%free()
deallocate(nzbr,idisp)
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
@ -245,20 +256,16 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call am1%mv_from(acoo1)
call am1%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
if (info == psb_success_) call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
if (info == psb_success_) & if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,& & p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build') call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999 goto 9999

206
mlprec/impl/mld_caggrmat_smth_asb.F90
Unescape Escape View File

@ -61,11 +61,6 @@
! according to the value of p%parms%aggr_omega_alg, specified by the user ! 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_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
! prolongator still requires a deep analysis and testing and its use is not
! recommended.
!
! The coarse-level matrix A_C is distributed among the parallel processes or ! 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, ! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_cprecinit and mld_zprecset. ! specified by the user through mld_cprecinit and mld_zprecset.
@ -116,20 +111,19 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_c_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_cspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
type(psb_cspmat_type) :: b integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: am1,am2, am3, am4 type(psb_cspmat_type) :: am3, am4
type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_spk_), allocatable :: adiag(:) complex(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat logical :: filter_mat
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer, parameter :: ncmax=16 integer, parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta real(psb_spk_) :: anorm, omega, tmp, dg, theta
@ -155,22 +149,10 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
ml_global_nmb = ( (p%parms%aggr_kind == mld_smooth_prol_).or.&
& ( (p%parms%aggr_kind == mld_biz_prol_).and.&
& (p%parms%coarse_mat == mld_repl_mat_)) )
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (p%parms%aggr_filter == mld_filter_mat_)
if (ml_global_nmb) then
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
@ -178,7 +160,6 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
end if
! naggr: number of local aggregates ! naggr: number of local aggregates
! nrow: local rows. ! nrow: local rows.
@ -202,7 +183,6 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = cone acoo%val(i) = cone
@ -210,15 +190,6 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
acoo%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call acoo%set_nzeros(ncol)
else
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo%val(i) = cone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo%set_nzeros(nrow)
endif
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(acoo,info)
@ -285,35 +256,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_eig == mld_max_norm_) then if (p%parms%aggr_eig == mld_max_norm_) then
if (p%parms%aggr_kind == mld_biz_prol_) then
!
! This only works with CSR
!
anorm = szero
dg = sone
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega p%parms%aggr_omega_val = omega
@ -409,38 +352,33 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call ptilde%free() call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
if (ml_global_nmb) then
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else
call psb_rwextd(ncol,am1,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_
if (p%parms%aggr_kind == mld_smooth_prol_) then call op_prol%transp(op_restr)
call am1%transp(am2) call op_restr%mv_to(acoo)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
i=0 i=0
! !
@ -457,28 +395,21 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end do end do
call acoo%set_nzeros(i) call acoo%set_nzeros(i)
call acoo%trim() call acoo%trim()
call am2%mv_from(acoo) call op_restr%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
goto 9999 goto 9999
end if end if
else
call am1%transp(am2)
endif
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
if (p%parms%aggr_kind == mld_smooth_prol_) then ! op_restr = ((i-wDA)Ptilde)^T
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,& call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else if (p%parms%aggr_kind == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999 goto 9999
@ -488,27 +419,22 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3' & 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = am2 x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
select case(p%parms%aggr_kind)
case(mld_smooth_prol_)
select case(p%parms%coarse_mat) select case(p%parms%coarse_mat)
case(mld_distr_mat_) case(mld_distr_mat_)
nzac = b%get_nzeros() nzl = b%get_nzeros()
nzl = nzac
call b%mv_to(bcoo) call b%mv_to(bcoo)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
@ -530,38 +456,37 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%set_ncols(p%desc_ac%get_local_cols()) call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb() call p%ac%set_asb()
if (info == psb_success_) deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if end if
if (np>1) then if (np>1) then
call am1%mv_to(acsr1) call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros() nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo) call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo) if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr') if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -577,73 +502,12 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
case(mld_biz_prol_)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call psb_move_alloc(b,p%ac,info)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Build desc_ac, ac')
goto 9999
end if
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999
end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end select end select
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_smooth_prol_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
@ -652,14 +516,14 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

17
mlprec/impl/mld_daggrmat_asb.f90
Unescape Escape View File

@ -113,6 +113,7 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info integer, intent(out) :: info
! Local variables ! Local variables
type(psb_dspmat_type) :: b, op_prol,op_restr
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -128,15 +129,23 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case (p%parms%aggr_kind) select case (p%parms%aggr_kind)
case (mld_no_smooth_) case (mld_no_smooth_)
call mld_aggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999 goto 9999
end if end if
case(mld_smooth_prol_,mld_biz_prol_) case(mld_smooth_prol_)
call mld_aggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_)
call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
@ -144,7 +153,7 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_min_energy_) case(mld_min_energy_)
call mld_aggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999

156
mlprec/impl/mld_daggrmat_minnrg_asb.F90
Unescape Escape View File

@ -116,15 +116,15 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_d_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
type(psb_dspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm integer :: ictxt,np,me, err_act, icomm
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: am1,am2, af, ptilde, rtilde, atran, atp, atdatp type(psb_dspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_dspmat_type) :: dat, datp, datdatp, atmp3 type(psb_dspmat_type) :: dat, datp, datdatp, atmp3
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
@ -280,7 +280,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call am3%mv_to(acsr3) call am3%mv_to(acsr3)
! Compute omega_int ! Compute omega_int
ommx = cmplx(szero,szero) ommx = cmplx(dzero,dzero)
do i=1, ncol do i=1, ncol
omi(i) = omp(ilaggr(i)) omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@ -354,17 +354,17 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*Af)Ptilde ! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i) ! Doing it this way means to consider diag(Af_i)
! !
! !
call psb_symbmm(af,ptilde,am1,info) call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(af,ptilde,am1) call psb_numbmm(af,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -390,16 +390,16 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*A)Ptilde ! op_prol = (I-w*D*A)Ptilde
! !
! !
call psb_symbmm(am3,ptilde,am1,info) call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(am3,ptilde,am1) call psb_numbmm(am3,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -458,7 +458,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
omp = omp/oden omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int ! Compute omega_int
ommx = cmplx(szero,szero) ommx = cmplx(dzero,dzero)
do i=1, ncol do i=1, ncol
omi(i) = omp(ilaggr(i)) omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@ -509,20 +509,20 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call rtilde%mv_from(tmpcoo) call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr') call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,am2,info) call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,am2) call psb_numbmm(rtilde,atmp,op_restr)
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
@ -530,7 +530,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Now we have to fix this. The only rows of B that are correct ! 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(:) ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
! !
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
@ -543,21 +543,21 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
end do end do
call tmpcoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'
@ -576,14 +576,13 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd' & 'Done sphalo/ rwxtd'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
call b%mv_to(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
&a_err='Build b = am2 x am3') &a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
@ -597,6 +596,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_distr_mat_) case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros() nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -627,29 +627,29 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%cscnv(info,type='csr') call p%ac%cscnv(info,type='csr')
if (np>1) then if (np>1) then
call am1%mv_to(acsr) call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros() nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr) call op_prol%mv_from(acsr)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -659,46 +659,54 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! !
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
nzbr(:) = 0 if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
nzbr(me+1) = bcoo%get_nzeros()
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
if (info /= psb_success_) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& 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 if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
& a_err=' from mpi_allgatherv')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call bcoo%free()
call tmpcoo%fix(info)
call p%ac%mv_from(tmpcoo)
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) goto 9999 if(info /= psb_success_) goto 9999
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
deallocate(nzbr,idisp,stat=info) !!$ nzbr(:) = 0
if (info /= psb_success_) then !!$ nzbr(me+1) = bcoo%get_nzeros()
info = psb_err_alloc_dealloc_ !!$
call psb_errpush(info,name) !!$ call psb_sum(ictxt,nzbr(1:np))
goto 9999 !!$ nzac = sum(nzbr)
end if !!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
!!$ if (info /= psb_success_) goto 9999
!!$
!!$ do ip=1,np
!!$ idispip) = sum(nzbr(1:ip-1))
!!$ enddo
!!$ ndx = nzbr(me+1)
!!$
!!$ 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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
!!$ & psb_mpi_ipk_integer,icomm,info)
!!$ if (info == psb_success_)&
!!$ & 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,&
!!$ & a_err=' from mpi_allgatherv')
!!$ goto 9999
!!$ end if
!!$
!!$ call bcoo%free()
!!$ call tmpcoo%fix(info)
!!$ call p%ac%mv_from(tmpcoo)
!!$ call p%ac%cscnv(info,type='csr')
!!$ if(info /= psb_success_) goto 9999
!!$
!!$ deallocate(nzbr,idisp,stat=info)
!!$ if (info /= psb_success_) then
!!$ info = psb_err_alloc_dealloc_
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
@ -715,13 +723,13 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => R i.e. restriction operator ! op_restr => R i.e. restriction operator
! am1 => P i.e. prolongation operator ! op_prol => P i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

125
mlprec/impl/mld_daggrmat_nosmth_asb.F90
Unescape Escape View File

@ -99,17 +99,17 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_d_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: b, op_prol,op_restr
! Local variables ! Local variables
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_dspmat_type) :: am1,am2 type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo
type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo type(psb_d_csr_sparse_mat) :: acsr1, acsr2
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, i & naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb' name='mld_aggrmat_nosmth_asb'
@ -128,32 +128,19 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1=sum(nlaggr(1:me)) naggrm1=sum(nlaggr(1:me))
if (p%parms%coarse_mat == mld_repl_mat_) then
do i=1, nrow do i=1, nrow
ilaggr(i) = ilaggr(i) + naggrm1 ilaggr(i) = ilaggr(i) + naggrm1
end do end do
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
end if
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
if (p%parms%coarse_mat == mld_repl_mat_) then
call acoo1%allocate(ncol,ntaggr,ncol) call acoo1%allocate(ncol,ntaggr,ncol)
else
call acoo1%allocate(ncol,naggr,ncol)
end if
do i=1,nrow do i=1,nrow
acoo1%val(i) = done acoo1%val(i) = done
@ -165,10 +152,13 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call acoo1%set_nzeros(nrow) call acoo1%set_nzeros(nrow)
call acoo1%set_asb() call acoo1%set_asb()
call acoo1%fix(info) call acoo1%fix(info)
call acoo1%transp(acoo2)
call a%csclip(bcoo,info,jmax=nrow)
call op_prol%mv_from(acoo1)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = bcoo%get_nzeros()
do i=1, nzt do i=1, nzt
@ -181,6 +171,8 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call bcoo%fix(info) call bcoo%fix(info)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then if (p%parms%coarse_mat == mld_repl_mat_) then
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
@ -189,55 +181,74 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
nzbr(:) = 0 else if (p%parms%coarse_mat == mld_distr_mat_) then
nzbr(me+1) = nzt
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call ac_coo%allocate(ntaggr,ntaggr,nzac) nzl = b%get_nzeros()
call b%mv_to(bcoo)
do ip=1,np if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
idisp(ip) = sum(nzbr(1:ip-1)) if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
enddo if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
ndx = nzbr(me+1) 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')
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,ac_coo%val,nzbr,idisp,&
& 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 if (info /= psb_success_) then
info=-1 call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name) & a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
call ac_coo%set_nzeros(nzac) if (debug_level >= psb_debug_outer_) &
call ac_coo%set_dupl(psb_dupl_add_) & write(debug_unit,*) me,' ',trim(name),&
call ac_coo%fix(info) & 'Assembld aux descr. distr.'
call p%ac%mv_from(ac_coo) call p%ac%mv_from(bcoo)
else if (p%parms%coarse_mat == mld_distr_mat_) then call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
call p%ac%mv_from(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac, desc_ac') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
end if
call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end if end if
call bcoo%free()
deallocate(nzbr,idisp)
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
@ -245,20 +256,16 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call am1%mv_from(acoo1)
call am1%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
if (info == psb_success_) call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
if (info == psb_success_) & if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,& & p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build') call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999 goto 9999

206
mlprec/impl/mld_daggrmat_smth_asb.F90
Unescape Escape View File

@ -61,11 +61,6 @@
! according to the value of p%parms%aggr_omega_alg, specified by the user ! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_dprecinit and mld_zprecset. ! 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
! prolongator still requires a deep analysis and testing and its use is not
! recommended.
!
! The coarse-level matrix A_C is distributed among the parallel processes or ! 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, ! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_dprecinit and mld_zprecset. ! specified by the user through mld_dprecinit and mld_zprecset.
@ -116,20 +111,19 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_d_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
type(psb_dspmat_type) :: b integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: am1,am2, am3, am4 type(psb_dspmat_type) :: am3, am4
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_dpk_), allocatable :: adiag(:) real(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat logical :: filter_mat
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer, parameter :: ncmax=16 integer, parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta real(psb_dpk_) :: anorm, omega, tmp, dg, theta
@ -155,22 +149,10 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
ml_global_nmb = ( (p%parms%aggr_kind == mld_smooth_prol_).or.&
& ( (p%parms%aggr_kind == mld_biz_prol_).and.&
& (p%parms%coarse_mat == mld_repl_mat_)) )
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (p%parms%aggr_filter == mld_filter_mat_)
if (ml_global_nmb) then
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
@ -178,7 +160,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
end if
! naggr: number of local aggregates ! naggr: number of local aggregates
! nrow: local rows. ! nrow: local rows.
@ -202,7 +183,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = done acoo%val(i) = done
@ -210,15 +190,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
acoo%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call acoo%set_nzeros(ncol)
else
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo%val(i) = done
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo%set_nzeros(nrow)
endif
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(acoo,info)
@ -285,35 +256,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_eig == mld_max_norm_) then if (p%parms%aggr_eig == mld_max_norm_) then
if (p%parms%aggr_kind == mld_biz_prol_) then
!
! This only works with CSR
!
anorm = dzero
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega p%parms%aggr_omega_val = omega
@ -409,38 +352,33 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call ptilde%free() call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
if (ml_global_nmb) then
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else
call psb_rwextd(ncol,am1,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_
if (p%parms%aggr_kind == mld_smooth_prol_) then call op_prol%transp(op_restr)
call am1%transp(am2) call op_restr%mv_to(acoo)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
i=0 i=0
! !
@ -457,28 +395,21 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end do end do
call acoo%set_nzeros(i) call acoo%set_nzeros(i)
call acoo%trim() call acoo%trim()
call am2%mv_from(acoo) call op_restr%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
goto 9999 goto 9999
end if end if
else
call am1%transp(am2)
endif
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
if (p%parms%aggr_kind == mld_smooth_prol_) then ! op_restr = ((i-wDA)Ptilde)^T
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,& call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else if (p%parms%aggr_kind == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999 goto 9999
@ -488,27 +419,22 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3' & 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = am2 x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
select case(p%parms%aggr_kind)
case(mld_smooth_prol_)
select case(p%parms%coarse_mat) select case(p%parms%coarse_mat)
case(mld_distr_mat_) case(mld_distr_mat_)
nzac = b%get_nzeros() nzl = b%get_nzeros()
nzl = nzac
call b%mv_to(bcoo) call b%mv_to(bcoo)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
@ -530,38 +456,37 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%set_ncols(p%desc_ac%get_local_cols()) call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb() call p%ac%set_asb()
if (info == psb_success_) deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if end if
if (np>1) then if (np>1) then
call am1%mv_to(acsr1) call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros() nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo) call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo) if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr') if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -577,73 +502,12 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
case(mld_biz_prol_)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call psb_move_alloc(b,p%ac,info)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Build desc_ac, ac')
goto 9999
end if
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999
end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end select end select
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_smooth_prol_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
@ -652,14 +516,14 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

1
mlprec/impl/mld_dprecinit.F90
Unescape Escape View File

@ -125,6 +125,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
! Do we want to do something? ! Do we want to do something?
endif endif
endif endif
p%coarse_aggr_size = -1
select case(psb_toupper(ptype(1:len_trim(ptype)))) select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NOPREC','NONE') case ('NOPREC','NONE')

5
mlprec/impl/mld_dprecset.F90
Unescape Escape View File

@ -129,6 +129,11 @@ subroutine mld_dprecseti(p,what,val,info,ilev)
return return
endif endif
if (what == mld_coarse_aggr_size_) then
p%coarse_aggr_size = max(val,-1)
return
end if
! !
! Set preconditioner parameters at level ilev. ! Set preconditioner parameters at level ilev.
! !

17
mlprec/impl/mld_saggrmat_asb.f90
Unescape Escape View File

@ -113,6 +113,7 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info integer, intent(out) :: info
! Local variables ! Local variables
type(psb_sspmat_type) :: b, op_prol,op_restr
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -128,15 +129,23 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case (p%parms%aggr_kind) select case (p%parms%aggr_kind)
case (mld_no_smooth_) case (mld_no_smooth_)
call mld_aggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999 goto 9999
end if end if
case(mld_smooth_prol_,mld_biz_prol_) case(mld_smooth_prol_)
call mld_aggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_)
call mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
@ -144,7 +153,7 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_min_energy_) case(mld_min_energy_)
call mld_aggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999

152
mlprec/impl/mld_saggrmat_minnrg_asb.F90
Unescape Escape View File

@ -116,15 +116,15 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_s_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_sspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
type(psb_sspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm integer :: ictxt,np,me, err_act, icomm
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: am1,am2, af, ptilde, rtilde, atran, atp, atdatp type(psb_sspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_sspmat_type) :: dat, datp, datdatp, atmp3 type(psb_sspmat_type) :: dat, datp, datdatp, atmp3
type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
@ -354,17 +354,17 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*Af)Ptilde ! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i) ! Doing it this way means to consider diag(Af_i)
! !
! !
call psb_symbmm(af,ptilde,am1,info) call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(af,ptilde,am1) call psb_numbmm(af,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -390,16 +390,16 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*A)Ptilde ! op_prol = (I-w*D*A)Ptilde
! !
! !
call psb_symbmm(am3,ptilde,am1,info) call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(am3,ptilde,am1) call psb_numbmm(am3,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -509,20 +509,20 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call rtilde%mv_from(tmpcoo) call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr') call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,am2,info) call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,am2) call psb_numbmm(rtilde,atmp,op_restr)
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
@ -530,7 +530,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Now we have to fix this. The only rows of B that are correct ! 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(:) ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
! !
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
@ -543,21 +543,21 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
end do end do
call tmpcoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'
@ -576,14 +576,13 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd' & 'Done sphalo/ rwxtd'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
call b%mv_to(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
&a_err='Build b = am2 x am3') &a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
@ -597,6 +596,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_distr_mat_) case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros() nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -627,29 +627,29 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%cscnv(info,type='csr') call p%ac%cscnv(info,type='csr')
if (np>1) then if (np>1) then
call am1%mv_to(acsr) call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros() nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr) call op_prol%mv_from(acsr)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -659,46 +659,54 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! !
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
nzbr(:) = 0 if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
nzbr(me+1) = bcoo%get_nzeros()
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
if (info /= psb_success_) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& 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 if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
& a_err=' from mpi_allgatherv')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call bcoo%free()
call tmpcoo%fix(info)
call p%ac%mv_from(tmpcoo)
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) goto 9999 if(info /= psb_success_) goto 9999
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
deallocate(nzbr,idisp,stat=info) !!$ nzbr(:) = 0
if (info /= psb_success_) then !!$ nzbr(me+1) = bcoo%get_nzeros()
info = psb_err_alloc_dealloc_ !!$
call psb_errpush(info,name) !!$ call psb_sum(ictxt,nzbr(1:np))
goto 9999 !!$ nzac = sum(nzbr)
end if !!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
!!$ if (info /= psb_success_) goto 9999
!!$
!!$ do ip=1,np
!!$ idispip) = sum(nzbr(1:ip-1))
!!$ enddo
!!$ ndx = nzbr(me+1)
!!$
!!$ 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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
!!$ & psb_mpi_ipk_integer,icomm,info)
!!$ if (info == psb_success_)&
!!$ & 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,&
!!$ & a_err=' from mpi_allgatherv')
!!$ goto 9999
!!$ end if
!!$
!!$ call bcoo%free()
!!$ call tmpcoo%fix(info)
!!$ call p%ac%mv_from(tmpcoo)
!!$ call p%ac%cscnv(info,type='csr')
!!$ if(info /= psb_success_) goto 9999
!!$
!!$ deallocate(nzbr,idisp,stat=info)
!!$ if (info /= psb_success_) then
!!$ info = psb_err_alloc_dealloc_
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
@ -715,13 +723,13 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => R i.e. restriction operator ! op_restr => R i.e. restriction operator
! am1 => P i.e. prolongation operator ! op_prol => P i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

125
mlprec/impl/mld_saggrmat_nosmth_asb.F90
Unescape Escape View File

@ -99,17 +99,17 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_s_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_sspmat_type) :: b, op_prol,op_restr
! Local variables ! Local variables
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_sspmat_type) :: am1,am2 type(psb_s_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo
type(psb_s_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo type(psb_s_csr_sparse_mat) :: acsr1, acsr2
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, i & naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb' name='mld_aggrmat_nosmth_asb'
@ -128,32 +128,19 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1=sum(nlaggr(1:me)) naggrm1=sum(nlaggr(1:me))
if (p%parms%coarse_mat == mld_repl_mat_) then
do i=1, nrow do i=1, nrow
ilaggr(i) = ilaggr(i) + naggrm1 ilaggr(i) = ilaggr(i) + naggrm1
end do end do
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
end if
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
if (p%parms%coarse_mat == mld_repl_mat_) then
call acoo1%allocate(ncol,ntaggr,ncol) call acoo1%allocate(ncol,ntaggr,ncol)
else
call acoo1%allocate(ncol,naggr,ncol)
end if
do i=1,nrow do i=1,nrow
acoo1%val(i) = sone acoo1%val(i) = sone
@ -165,10 +152,13 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call acoo1%set_nzeros(nrow) call acoo1%set_nzeros(nrow)
call acoo1%set_asb() call acoo1%set_asb()
call acoo1%fix(info) call acoo1%fix(info)
call acoo1%transp(acoo2)
call a%csclip(bcoo,info,jmax=nrow)
call op_prol%mv_from(acoo1)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = bcoo%get_nzeros()
do i=1, nzt do i=1, nzt
@ -181,6 +171,8 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call bcoo%fix(info) call bcoo%fix(info)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then if (p%parms%coarse_mat == mld_repl_mat_) then
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
@ -189,55 +181,74 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
nzbr(:) = 0 else if (p%parms%coarse_mat == mld_distr_mat_) then
nzbr(me+1) = nzt
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call ac_coo%allocate(ntaggr,ntaggr,nzac) nzl = b%get_nzeros()
call b%mv_to(bcoo)
do ip=1,np if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
idisp(ip) = sum(nzbr(1:ip-1)) if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
enddo if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
ndx = nzbr(me+1) 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')
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 if (info /= psb_success_) then
info=-1 call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name) & a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
call ac_coo%set_nzeros(nzac) if (debug_level >= psb_debug_outer_) &
call ac_coo%set_dupl(psb_dupl_add_) & write(debug_unit,*) me,' ',trim(name),&
call ac_coo%fix(info) & 'Assembld aux descr. distr.'
call p%ac%mv_from(ac_coo) call p%ac%mv_from(bcoo)
else if (p%parms%coarse_mat == mld_distr_mat_) then call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
call p%ac%mv_from(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac, desc_ac') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
end if
call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end if end if
call bcoo%free()
deallocate(nzbr,idisp)
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
@ -245,20 +256,16 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call am1%mv_from(acoo1)
call am1%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
if (info == psb_success_) call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
if (info == psb_success_) & if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,& & p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build') call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999 goto 9999

206
mlprec/impl/mld_saggrmat_smth_asb.F90
Unescape Escape View File

@ -61,11 +61,6 @@
! according to the value of p%parms%aggr_omega_alg, specified by the user ! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_sprecinit and mld_zprecset. ! 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
! prolongator still requires a deep analysis and testing and its use is not
! recommended.
!
! The coarse-level matrix A_C is distributed among the parallel processes or ! 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, ! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_sprecinit and mld_zprecset. ! specified by the user through mld_sprecinit and mld_zprecset.
@ -116,20 +111,19 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_s_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_sspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
type(psb_sspmat_type) :: b integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: am1,am2, am3, am4 type(psb_sspmat_type) :: am3, am4
type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_spk_), allocatable :: adiag(:) real(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat logical :: filter_mat
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer, parameter :: ncmax=16 integer, parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta real(psb_spk_) :: anorm, omega, tmp, dg, theta
@ -155,22 +149,10 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
ml_global_nmb = ( (p%parms%aggr_kind == mld_smooth_prol_).or.&
& ( (p%parms%aggr_kind == mld_biz_prol_).and.&
& (p%parms%coarse_mat == mld_repl_mat_)) )
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (p%parms%aggr_filter == mld_filter_mat_)
if (ml_global_nmb) then
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
@ -178,7 +160,6 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
end if
! naggr: number of local aggregates ! naggr: number of local aggregates
! nrow: local rows. ! nrow: local rows.
@ -202,7 +183,6 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = sone acoo%val(i) = sone
@ -210,15 +190,6 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
acoo%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call acoo%set_nzeros(ncol)
else
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo%val(i) = sone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo%set_nzeros(nrow)
endif
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(acoo,info)
@ -285,35 +256,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_eig == mld_max_norm_) then if (p%parms%aggr_eig == mld_max_norm_) then
if (p%parms%aggr_kind == mld_biz_prol_) then
!
! This only works with CSR
!
anorm = szero
dg = sone
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega p%parms%aggr_omega_val = omega
@ -409,38 +352,33 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call ptilde%free() call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
if (ml_global_nmb) then
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else
call psb_rwextd(ncol,am1,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_
if (p%parms%aggr_kind == mld_smooth_prol_) then call op_prol%transp(op_restr)
call am1%transp(am2) call op_restr%mv_to(acoo)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
i=0 i=0
! !
@ -457,28 +395,21 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end do end do
call acoo%set_nzeros(i) call acoo%set_nzeros(i)
call acoo%trim() call acoo%trim()
call am2%mv_from(acoo) call op_restr%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
goto 9999 goto 9999
end if end if
else
call am1%transp(am2)
endif
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
if (p%parms%aggr_kind == mld_smooth_prol_) then ! op_restr = ((i-wDA)Ptilde)^T
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,& call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else if (p%parms%aggr_kind == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999 goto 9999
@ -488,27 +419,22 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3' & 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = am2 x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
select case(p%parms%aggr_kind)
case(mld_smooth_prol_)
select case(p%parms%coarse_mat) select case(p%parms%coarse_mat)
case(mld_distr_mat_) case(mld_distr_mat_)
nzac = b%get_nzeros() nzl = b%get_nzeros()
nzl = nzac
call b%mv_to(bcoo) call b%mv_to(bcoo)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
@ -530,38 +456,37 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%set_ncols(p%desc_ac%get_local_cols()) call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb() call p%ac%set_asb()
if (info == psb_success_) deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if end if
if (np>1) then if (np>1) then
call am1%mv_to(acsr1) call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros() nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo) call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo) if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr') if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -577,73 +502,12 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end select end select
case(mld_biz_prol_)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call psb_move_alloc(b,p%ac,info)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Build desc_ac, ac')
goto 9999
end if
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999
end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_smooth_prol_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
@ -652,14 +516,14 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

17
mlprec/impl/mld_zaggrmat_asb.f90
Unescape Escape View File

@ -113,6 +113,7 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info integer, intent(out) :: info
! Local variables ! Local variables
type(psb_zspmat_type) :: b, op_prol,op_restr
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -128,15 +129,23 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case (p%parms%aggr_kind) select case (p%parms%aggr_kind)
case (mld_no_smooth_) case (mld_no_smooth_)
call mld_aggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999 goto 9999
end if end if
case(mld_smooth_prol_,mld_biz_prol_) case(mld_smooth_prol_)
call mld_aggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_)
call mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
@ -144,7 +153,7 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_min_energy_) case(mld_min_energy_)
call mld_aggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb') call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999

156
mlprec/impl/mld_zaggrmat_minnrg_asb.F90
Unescape Escape View File

@ -116,15 +116,15 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_z_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
type(psb_zspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm integer :: ictxt,np,me, err_act, icomm
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: am1,am2, af, ptilde, rtilde, atran, atp, atdatp type(psb_zspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_zspmat_type) :: dat, datp, datdatp, atmp3 type(psb_zspmat_type) :: dat, datp, datdatp, atmp3
type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
@ -280,7 +280,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call am3%mv_to(acsr3) call am3%mv_to(acsr3)
! Compute omega_int ! Compute omega_int
ommx = cmplx(szero,szero) ommx = cmplx(dzero,dzero)
do i=1, ncol do i=1, ncol
omi(i) = omp(ilaggr(i)) omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@ -354,17 +354,17 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*Af)Ptilde ! op_prol = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i) ! Doing it this way means to consider diag(Af_i)
! !
! !
call psb_symbmm(af,ptilde,am1,info) call psb_symbmm(af,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(af,ptilde,am1) call psb_numbmm(af,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -390,16 +390,16 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Symbmm90 does the allocation for its result. ! Symbmm90 does the allocation for its result.
! !
! am1 = (I-w*D*A)Ptilde ! op_prol = (I-w*D*A)Ptilde
! !
! !
call psb_symbmm(am3,ptilde,am1,info) call psb_symbmm(am3,ptilde,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999 goto 9999
end if end if
call psb_numbmm(am3,ptilde,am1) call psb_numbmm(am3,ptilde,op_prol)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -458,7 +458,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
omp = omp/oden omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10)) ! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int ! Compute omega_int
ommx = cmplx(szero,szero) ommx = cmplx(dzero,dzero)
do i=1, ncol do i=1, ncol
omi(i) = omp(ilaggr(i)) omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i) if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@ -509,20 +509,20 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call rtilde%mv_from(tmpcoo) call rtilde%mv_from(tmpcoo)
call rtilde%cscnv(info,type='csr') call rtilde%cscnv(info,type='csr')
call psb_symbmm(rtilde,atmp,am2,info) call psb_symbmm(rtilde,atmp,op_restr,info)
call psb_numbmm(rtilde,atmp,am2) call psb_numbmm(rtilde,atmp,op_restr)
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
@ -530,7 +530,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Now we have to fix this. The only rows of B that are correct ! 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(:) ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
! !
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
@ -543,21 +543,21 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
end do end do
call tmpcoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'
@ -576,14 +576,13 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done sphalo/ rwxtd' & 'Done sphalo/ rwxtd'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
call b%mv_to(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
&a_err='Build b = am2 x am3') &a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
@ -597,6 +596,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
case(mld_distr_mat_) case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros() nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -627,29 +627,29 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%cscnv(info,type='csr') call p%ac%cscnv(info,type='csr')
if (np>1) then if (np>1) then
call am1%mv_to(acsr) call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros() nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr) call op_prol%mv_from(acsr)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%mv_to(tmpcoo) call op_restr%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros() nzl = tmpcoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
call am2%mv_from(tmpcoo) call op_restr%mv_from(tmpcoo)
call am2%cscnv(info,type='csr') call op_restr%cscnv(info,type='csr')
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -659,46 +659,54 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! !
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
nzbr(:) = 0 if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
nzbr(me+1) = bcoo%get_nzeros()
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
if (info /= psb_success_) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& 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 if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
& a_err=' from mpi_allgatherv')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
call bcoo%free()
call tmpcoo%fix(info)
call p%ac%mv_from(tmpcoo)
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) goto 9999 if(info /= psb_success_) goto 9999
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
deallocate(nzbr,idisp,stat=info) !!$ nzbr(:) = 0
if (info /= psb_success_) then !!$ nzbr(me+1) = bcoo%get_nzeros()
info = psb_err_alloc_dealloc_ !!$
call psb_errpush(info,name) !!$ call psb_sum(ictxt,nzbr(1:np))
goto 9999 !!$ nzac = sum(nzbr)
end if !!$ if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
!!$ if (info /= psb_success_) goto 9999
!!$
!!$ do ip=1,np
!!$ idispip) = sum(nzbr(1:ip-1))
!!$ enddo
!!$ ndx = nzbr(me+1)
!!$
!!$ 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_ipk_integer,tmpcoo%ia,nzbr,idisp,&
!!$ & psb_mpi_ipk_integer,icomm,info)
!!$ if (info == psb_success_)&
!!$ & 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,&
!!$ & a_err=' from mpi_allgatherv')
!!$ goto 9999
!!$ end if
!!$
!!$ call bcoo%free()
!!$ call tmpcoo%fix(info)
!!$ call p%ac%mv_from(tmpcoo)
!!$ call p%ac%cscnv(info,type='csr')
!!$ if(info /= psb_success_) goto 9999
!!$
!!$ deallocate(nzbr,idisp,stat=info)
!!$ if (info /= psb_success_) then
!!$ info = psb_err_alloc_dealloc_
!!$ call psb_errpush(info,name)
!!$ goto 9999
!!$ end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
@ -715,13 +723,13 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => R i.e. restriction operator ! op_restr => R i.e. restriction operator
! am1 => P i.e. prolongation operator ! op_prol => P i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

125
mlprec/impl/mld_zaggrmat_nosmth_asb.F90
Unescape Escape View File

@ -99,17 +99,17 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_z_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: b, op_prol,op_restr
! Local variables ! Local variables
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: b
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_zspmat_type) :: am1,am2 type(psb_z_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo
type(psb_z_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo type(psb_z_csr_sparse_mat) :: acsr1, acsr2
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, & integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, i & naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb' name='mld_aggrmat_nosmth_asb'
@ -128,32 +128,19 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1=sum(nlaggr(1:me)) naggrm1=sum(nlaggr(1:me))
if (p%parms%coarse_mat == mld_repl_mat_) then
do i=1, nrow do i=1, nrow
ilaggr(i) = ilaggr(i) + naggrm1 ilaggr(i) = ilaggr(i) + naggrm1
end do end do
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
end if
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
if (p%parms%coarse_mat == mld_repl_mat_) then
call acoo1%allocate(ncol,ntaggr,ncol) call acoo1%allocate(ncol,ntaggr,ncol)
else
call acoo1%allocate(ncol,naggr,ncol)
end if
do i=1,nrow do i=1,nrow
acoo1%val(i) = zone acoo1%val(i) = zone
@ -165,10 +152,13 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call acoo1%set_nzeros(nrow) call acoo1%set_nzeros(nrow)
call acoo1%set_asb() call acoo1%set_asb()
call acoo1%fix(info) call acoo1%fix(info)
call acoo1%transp(acoo2)
call a%csclip(bcoo,info,jmax=nrow)
call op_prol%mv_from(acoo1)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = bcoo%get_nzeros()
do i=1, nzt do i=1, nzt
@ -181,6 +171,8 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call bcoo%fix(info) call bcoo%fix(info)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then if (p%parms%coarse_mat == mld_repl_mat_) then
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.) call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
@ -189,55 +181,74 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999 goto 9999
end if end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
nzbr(:) = 0 else if (p%parms%coarse_mat == mld_distr_mat_) then
nzbr(me+1) = nzt
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
call ac_coo%allocate(ntaggr,ntaggr,nzac) nzl = b%get_nzeros()
call b%mv_to(bcoo)
do ip=1,np if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
idisp(ip) = sum(nzbr(1:ip-1)) if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
enddo if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
ndx = nzbr(me+1) 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')
call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,ac_coo%val,nzbr,idisp,&
& 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 if (info /= psb_success_) then
info=-1 call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name) & a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
call ac_coo%set_nzeros(nzac) if (debug_level >= psb_debug_outer_) &
call ac_coo%set_dupl(psb_dupl_add_) & write(debug_unit,*) me,' ',trim(name),&
call ac_coo%fix(info) & 'Assembld aux descr. distr.'
call p%ac%mv_from(ac_coo) call p%ac%mv_from(bcoo)
else if (p%parms%coarse_mat == mld_distr_mat_) then call p%ac%set_nrows(p%desc_ac%get_local_rows())
call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb()
call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
call p%ac%mv_from(bcoo)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac, desc_ac') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if
if (np>1) then
call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr1)
endif endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
end if
call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_') call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end if end if
call bcoo%free()
deallocate(nzbr,idisp)
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
@ -245,20 +256,16 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call am1%mv_from(acoo1)
call am1%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
if (info == psb_success_) call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
if (info == psb_success_) & if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,& & p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build') call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999 goto 9999

206
mlprec/impl/mld_zaggrmat_smth_asb.F90
Unescape Escape View File

@ -61,11 +61,6 @@
! according to the value of p%parms%aggr_omega_alg, specified by the user ! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_zprecinit and mld_zprecset. ! through mld_zprecinit 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
! prolongator still requires a deep analysis and testing and its use is not
! recommended.
!
! The coarse-level matrix A_C is distributed among the parallel processes or ! 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, ! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_zprecinit and mld_zprecset. ! specified by the user through mld_zprecinit and mld_zprecset.
@ -116,20 +111,19 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_z_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
type(psb_zspmat_type) :: b integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: am1,am2, am3, am4 type(psb_zspmat_type) :: am3, am4
type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_dpk_), allocatable :: adiag(:) complex(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat logical :: filter_mat
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer, parameter :: ncmax=16 integer, parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta real(psb_dpk_) :: anorm, omega, tmp, dg, theta
@ -155,22 +149,10 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
ml_global_nmb = ( (p%parms%aggr_kind == mld_smooth_prol_).or.&
& ( (p%parms%aggr_kind == mld_biz_prol_).and.&
& (p%parms%coarse_mat == mld_repl_mat_)) )
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (p%parms%aggr_filter == mld_filter_mat_)
if (ml_global_nmb) then
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1 ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info) call psb_halo(ilaggr,desc_a,info)
@ -178,7 +160,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999 goto 9999
end if end if
end if
! naggr: number of local aggregates ! naggr: number of local aggregates
! nrow: local rows. ! nrow: local rows.
@ -202,7 +183,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = zone acoo%val(i) = zone
@ -210,15 +190,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
acoo%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call acoo%set_nzeros(ncol)
else
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo%val(i) = zone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo%set_nzeros(nrow)
endif
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(acoo,info)
@ -285,35 +256,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_eig == mld_max_norm_) then if (p%parms%aggr_eig == mld_max_norm_) then
if (p%parms%aggr_kind == mld_biz_prol_) then
!
! This only works with CSR
!
anorm = dzero
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
else
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega p%parms%aggr_omega_val = omega
@ -409,38 +352,33 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call ptilde%free() call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
if (ml_global_nmb) then
! !
! Now we have to gather the halo of am1, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call psb_sphalo(am1,desc_a,am4,info,& call psb_sphalo(op_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else
call psb_rwextd(ncol,am1,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,am1,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,am1,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',p%parms%aggr_kind, mld_smooth_prol_
if (p%parms%aggr_kind == mld_smooth_prol_) then call op_prol%transp(op_restr)
call am1%transp(am2) call op_restr%mv_to(acoo)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
i=0 i=0
! !
@ -457,28 +395,21 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end do end do
call acoo%set_nzeros(i) call acoo%set_nzeros(i)
call acoo%trim() call acoo%trim()
call am2%mv_from(acoo) call op_restr%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
goto 9999 goto 9999
end if end if
else
call am1%transp(am2)
endif
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
if (p%parms%aggr_kind == mld_smooth_prol_) then ! op_restr = ((i-wDA)Ptilde)^T
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,& call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
else if (p%parms%aggr_kind == mld_biz_prol_) then
call psb_rwextd(ncol,am3,info)
endif
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999 goto 9999
@ -488,27 +419,22 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting symbmm 3' & 'starting symbmm 3'
call psb_symbmm(am2,am3,b,info) call psb_symbmm(op_restr,am3,b,info)
if (info == psb_success_) call psb_numbmm(am2,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,b)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_) if (info == psb_success_) call b%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = am2 x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
end if end if
select case(p%parms%aggr_kind)
case(mld_smooth_prol_)
select case(p%parms%coarse_mat) select case(p%parms%coarse_mat)
case(mld_distr_mat_) case(mld_distr_mat_)
nzac = b%get_nzeros() nzl = b%get_nzeros()
nzl = nzac
call b%mv_to(bcoo) call b%mv_to(bcoo)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1)) if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
@ -530,38 +456,37 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%set_ncols(p%desc_ac%get_local_cols()) call p%ac%set_ncols(p%desc_ac%get_local_cols())
call p%ac%set_asb() call p%ac%set_asb()
if (info == psb_success_) deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
goto 9999 goto 9999
end if end if
if (np>1) then if (np>1) then
call am1%mv_to(acsr1) call op_prol%mv_to(acsr1)
nzl = acsr1%get_nzeros() nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I') call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc') call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999 goto 9999
end if end if
call am1%mv_from(acsr1) call op_prol%mv_from(acsr1)
endif endif
call am1%set_ncols(p%desc_ac%get_local_cols()) call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_) call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo) call op_restr%mv_to(acoo)
nzl = acoo%get_nzeros() nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I') if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo) if (info == psb_success_) call op_restr%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr') if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local') call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999 goto 9999
end if end if
end if end if
call am2%set_nrows(p%desc_ac%get_local_cols()) call op_restr%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
@ -577,73 +502,12 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
case(mld_biz_prol_)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call psb_move_alloc(b,p%ac,info)
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=naggr)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Build desc_ac, ac')
goto 9999
end if
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_cdall')
goto 9999
end if
call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if(info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default case default
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_') call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999 goto 9999
end select end select
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_smooth_prol_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_) call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
@ -652,14 +516,14 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! !
! Copy the prolongation/restriction matrices into the descriptor map. ! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => PR^T i.e. restriction operator ! op_restr => PR^T i.e. restriction operator
! am1 => PR i.e. prolongation operator ! op_prol => PR i.e. prolongation operator
! !
p%map = psb_linmap(psb_map_aggr_,desc_a,& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr) & p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call am1%free() if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call am2%free() if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free') call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999

1
mlprec/mld_base_prec_type.F90
Unescape Escape View File

@ -157,6 +157,7 @@ module mld_base_prec_type
integer, parameter :: mld_coarse_fillin_ = 32 integer, parameter :: mld_coarse_fillin_ = 32
integer, parameter :: mld_coarse_subsolve_ = 33 integer, parameter :: mld_coarse_subsolve_ = 33
integer, parameter :: mld_smoother_sweeps_ = 34 integer, parameter :: mld_smoother_sweeps_ = 34
integer, parameter :: mld_coarse_aggr_size_ = 35
integer, parameter :: mld_ifpsz_ = 36 integer, parameter :: mld_ifpsz_ = 36
! !

36
mlprec/mld_c_inner_mod.f90
Unescape Escape View File

@ -130,6 +130,7 @@ module mld_c_inner_mod
end subroutine mld_c_dec_map_bld end subroutine mld_c_dec_map_bld
end interface mld_dec_map_bld end interface mld_dec_map_bld
interface mld_aggrmat_asb interface mld_aggrmat_asb
subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
@ -142,32 +143,10 @@ module mld_c_inner_mod
end subroutine mld_caggrmat_asb end subroutine mld_caggrmat_asb
end interface mld_aggrmat_asb end interface mld_aggrmat_asb
interface mld_aggrmat_nosmth_asb
subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_c_prec_type, only : mld_c_onelev_type
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_caggrmat_nosmth_asb
end interface mld_aggrmat_nosmth_asb
interface mld_aggrmat_smth_asb
subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_c_prec_type, only : mld_c_onelev_type
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_caggrmat_smth_asb
end interface mld_aggrmat_smth_asb
interface mld_aggrmat_minnrg_asb abstract interface
subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_
use mld_c_prec_type, only : mld_c_onelev_type use mld_c_prec_type, only : mld_c_onelev_type
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
@ -175,7 +154,12 @@ module mld_c_inner_mod
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_c_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_caggrmat_minnrg_asb end subroutine mld_caggrmat_var_asb
end interface mld_aggrmat_minnrg_asb end interface
procedure(mld_caggrmat_var_asb) :: mld_caggrmat_nosmth_asb, mld_caggrmat_smth_asb,&
& mld_caggrmat_minnrg_asb, mld_caggrmat_biz_asb
end module mld_c_inner_mod end module mld_c_inner_mod

36
mlprec/mld_d_inner_mod.f90
Unescape Escape View File

@ -130,6 +130,7 @@ module mld_d_inner_mod
end subroutine mld_d_dec_map_bld end subroutine mld_d_dec_map_bld
end interface mld_dec_map_bld end interface mld_dec_map_bld
interface mld_aggrmat_asb interface mld_aggrmat_asb
subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
@ -142,32 +143,10 @@ module mld_d_inner_mod
end subroutine mld_daggrmat_asb end subroutine mld_daggrmat_asb
end interface mld_aggrmat_asb end interface mld_aggrmat_asb
interface mld_aggrmat_nosmth_asb
subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_d_prec_type, only : mld_d_onelev_type
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_daggrmat_nosmth_asb
end interface mld_aggrmat_nosmth_asb
interface mld_aggrmat_smth_asb
subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_d_prec_type, only : mld_d_onelev_type
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_daggrmat_smth_asb
end interface mld_aggrmat_smth_asb
interface mld_aggrmat_minnrg_asb abstract interface
subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
use mld_d_prec_type, only : mld_d_onelev_type use mld_d_prec_type, only : mld_d_onelev_type
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
@ -175,7 +154,12 @@ module mld_d_inner_mod
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_d_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_daggrmat_minnrg_asb end subroutine mld_daggrmat_var_asb
end interface mld_aggrmat_minnrg_asb end interface
procedure(mld_daggrmat_var_asb) :: mld_daggrmat_nosmth_asb, mld_daggrmat_smth_asb,&
& mld_daggrmat_minnrg_asb, mld_daggrmat_biz_asb
end module mld_d_inner_mod end module mld_d_inner_mod

1
mlprec/mld_d_prec_type.f90
Unescape Escape View File

@ -81,6 +81,7 @@ module mld_d_prec_type
type, extends(psb_dprec_type) :: mld_dprec_type type, extends(psb_dprec_type) :: mld_dprec_type
integer :: ictxt integer :: ictxt
integer :: coarse_aggr_size
real(psb_dpk_) :: op_complexity=dzero real(psb_dpk_) :: op_complexity=dzero
type(mld_d_onelev_type), allocatable :: precv(:) type(mld_d_onelev_type), allocatable :: precv(:)
contains contains

36
mlprec/mld_s_inner_mod.f90
Unescape Escape View File

@ -130,6 +130,7 @@ module mld_s_inner_mod
end subroutine mld_s_dec_map_bld end subroutine mld_s_dec_map_bld
end interface mld_dec_map_bld end interface mld_dec_map_bld
interface mld_aggrmat_asb interface mld_aggrmat_asb
subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
@ -142,32 +143,10 @@ module mld_s_inner_mod
end subroutine mld_saggrmat_asb end subroutine mld_saggrmat_asb
end interface mld_aggrmat_asb end interface mld_aggrmat_asb
interface mld_aggrmat_nosmth_asb
subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_s_prec_type, only : mld_s_onelev_type
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_saggrmat_nosmth_asb
end interface mld_aggrmat_nosmth_asb
interface mld_aggrmat_smth_asb
subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_s_prec_type, only : mld_s_onelev_type
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_saggrmat_smth_asb
end interface mld_aggrmat_smth_asb
interface mld_aggrmat_minnrg_asb abstract interface
subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_
use mld_s_prec_type, only : mld_s_onelev_type use mld_s_prec_type, only : mld_s_onelev_type
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
@ -175,7 +154,12 @@ module mld_s_inner_mod
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_s_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_saggrmat_minnrg_asb end subroutine mld_saggrmat_var_asb
end interface mld_aggrmat_minnrg_asb end interface
procedure(mld_saggrmat_var_asb) :: mld_saggrmat_nosmth_asb, mld_saggrmat_smth_asb,&
& mld_saggrmat_minnrg_asb, mld_saggrmat_biz_asb
end module mld_s_inner_mod end module mld_s_inner_mod

36
mlprec/mld_z_inner_mod.f90
Unescape Escape View File

@ -130,6 +130,7 @@ module mld_z_inner_mod
end subroutine mld_z_dec_map_bld end subroutine mld_z_dec_map_bld
end interface mld_dec_map_bld end interface mld_dec_map_bld
interface mld_aggrmat_asb interface mld_aggrmat_asb
subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
@ -142,32 +143,10 @@ module mld_z_inner_mod
end subroutine mld_zaggrmat_asb end subroutine mld_zaggrmat_asb
end interface mld_aggrmat_asb end interface mld_aggrmat_asb
interface mld_aggrmat_nosmth_asb
subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_z_prec_type, only : mld_z_onelev_type
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_zaggrmat_nosmth_asb
end interface mld_aggrmat_nosmth_asb
interface mld_aggrmat_smth_asb
subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_z_prec_type, only : mld_z_onelev_type
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p
integer, intent(out) :: info
end subroutine mld_zaggrmat_smth_asb
end interface mld_aggrmat_smth_asb
interface mld_aggrmat_minnrg_asb abstract interface
subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_
use mld_z_prec_type, only : mld_z_onelev_type use mld_z_prec_type, only : mld_z_onelev_type
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
@ -175,7 +154,12 @@ module mld_z_inner_mod
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_z_onelev_type), intent(inout), target :: p
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_zaggrmat_minnrg_asb end subroutine mld_zaggrmat_var_asb
end interface mld_aggrmat_minnrg_asb end interface
procedure(mld_zaggrmat_var_asb) :: mld_zaggrmat_nosmth_asb, mld_zaggrmat_smth_asb,&
& mld_zaggrmat_minnrg_asb, mld_zaggrmat_biz_asb
end module mld_z_inner_mod end module mld_z_inner_mod

2
tests/pdegen/runs/ppde.inp
Unescape Escape View File

@ -17,7 +17,7 @@ ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU
4 ! Smoother/Jacobi sweeps 4 ! Smoother/Jacobi sweeps
BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML
3 ! Number of levels in a multilevel preconditioner 3 ! Number of levels in a multilevel preconditioner
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED NONSMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED
DEC ! Type of aggregation DEC SYMDEC GLB DEC ! Type of aggregation DEC SYMDEC GLB
MULT ! Type of multilevel correction: ADD MULT MULT ! Type of multilevel correction: ADD MULT
TWOSIDE ! Side of correction PRE POST TWOSIDE (ignored for ADD) TWOSIDE ! Side of correction PRE POST TWOSIDE (ignored for ADD)

Write Preview
Loading…
Cancel
Save