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

mld2p4-NewNL:

Browse Source 
mlprec/impl/mld_caggrmat_asb.f90
 mlprec/impl/mld_caggrmat_biz_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_biz_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_saggrmat_asb.f90
 mlprec/impl/mld_saggrmat_biz_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_biz_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_c_inner_mod.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_s_inner_mod.f90
 mlprec/mld_z_inner_mod.f90

A bit of internal cleanup.
stopcriterion
Salvatore Filippone 13 years ago
parent acb96d42cf
commit 788d394f58
24 changed files with 972 additions and 2348 deletions
Show Stats Download Patch File Download Diff File

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

@ -113,7 +113,11 @@ 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 type(psb_cspmat_type) :: ac, op_prol,op_restr
type(psb_c_coo_sparse_mat) :: acoo, bcoo
type(psb_c_csr_sparse_mat) :: acsr1
integer :: nzl,ntaggr
integer :: debug_level, debug_unit
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -121,6 +125,9 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(psb_get_errstatus().ne.0) return if(psb_get_errstatus().ne.0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context() ictxt = desc_a%get_context()
@ -129,42 +136,138 @@ 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_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,&
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999
end if
case(mld_smooth_prol_) case(mld_smooth_prol_)
call mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_) case(mld_biz_prol_)
call mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
case(mld_min_energy_)
call mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, &
& p%parms,ac,op_prol,op_restr,info)
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Invalid aggr kind')
goto 9999
end select
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb')
goto 9999
end if
ntaggr = sum(nlaggr)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = ac%get_nzeros()
call ac%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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
case(mld_min_energy_) 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 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='psb_sp_free')
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
end if 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
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 '
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_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999 goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

165
mlprec/impl/mld_caggrmat_biz_asb.f90
Unescape Escape View File

@ -78,7 +78,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_biz_asb use mld_c_inner_mod, mld_protect_name => mld_caggrmat_biz_asb
@ -88,9 +88,9 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_cspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -98,7 +98,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: am3, am4 type(psb_cspmat_type) :: am3, am4
type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_c_coo_sparse_mat) :: tmpcoo
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)
@ -123,7 +123,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -131,7 +131,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -163,16 +163,16 @@ subroutine mld_caggrmat_biz_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
call acoo%allocate(ncol,naggr,ncol) call tmpcoo%allocate(ncol,naggr,ncol)
do i=1,nrow do i=1,nrow
acoo%val(i) = cone tmpcoo%val(i) = cone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(nrow) call tmpcoo%set_nzeros(nrow)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -203,19 +203,19 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= czero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= czero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -232,9 +232,9 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
! !
! This only works with CSR ! This only works with CSR
@ -261,7 +261,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -269,11 +269,11 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -372,7 +372,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -389,10 +389,10 @@ subroutine mld_caggrmat_biz_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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
@ -401,107 +401,6 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

184
mlprec/impl/mld_caggrmat_minnrg_asb.f90
Unescape Escape View File

@ -98,7 +98,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_minnrg_asb use mld_c_inner_mod, mld_protect_name => mld_caggrmat_minnrg_asb
@ -108,9 +108,9 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_cspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
@ -121,8 +121,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_cspmat_type) :: 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) :: tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
complex(psb_spk_), allocatable :: adiag(:), adinv(:) complex(psb_spk_), allocatable :: adiag(:), adinv(:)
complex(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:) complex(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
@ -150,7 +150,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -165,7 +165,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -207,16 +207,16 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = cone tmpcoo%val(i) = cone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acoo%set_asb() call tmpcoo%set_asb()
call ptilde%mv_from(acoo) call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr') call ptilde%cscnv(info,type='csr')
!!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1')
@ -570,164 +570,18 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then 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 = op_restr x am3') &a_err='Build ac = op_restr x am3')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done mv_to_coo'
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' B matrix nzl: ',nzl
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Assembld aux descr. distr.'
call bcoo%set_nrows(p%desc_ac%get_local_rows())
call bcoo%set_ncols(p%desc_ac%get_local_cols())
call bcoo%fix(info)
call p%ac%mv_from(bcoo)
call p%ac%set_asb()
call p%ac%cscnv(info,type='csr')
if (np>1) then
call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%mv_to(tmpcoo)
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_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr')
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 '
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
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
!!$ nzbr(:) = 0
!!$ 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
!!$ 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
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => R i.e. restriction operator
! op_prol => P i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

144
mlprec/impl/mld_caggrmat_nosmth_asb.f90
Unescape Escape View File

@ -81,7 +81,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_nosmth_asb use mld_c_inner_mod, mld_protect_name => mld_caggrmat_nosmth_asb
@ -91,16 +91,16 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
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
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_c_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo type(psb_c_coo_sparse_mat) :: ac_coo, acoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2 type(psb_c_csr_sparse_mat) :: acsr1, acsr2
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
@ -134,136 +134,36 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call acoo1%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,nrow do i=1,nrow
acoo1%val(i) = cone acoo%val(i) = cone
acoo1%ia(i) = i acoo%ia(i) = i
acoo1%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo1%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call acoo1%set_nzeros(nrow) call acoo%set_nzeros(nrow)
call acoo1%set_asb() call acoo%set_asb()
call acoo1%fix(info) call acoo%fix(info)
call op_prol%mv_from(acoo1) call op_prol%mv_from(acoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr) if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow) call a%csclip(ac_coo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = ac_coo%get_nzeros()
do i=1, nzt do i=1, nzt
bcoo%ia(i) = ilaggr(bcoo%ia(i)) ac_coo%ia(i) = ilaggr(ac_coo%ia(i))
bcoo%ja(i) = ilaggr(bcoo%ja(i)) ac_coo%ja(i) = ilaggr(ac_coo%ja(i))
enddo enddo
call bcoo%set_nrows(naggr) call ac_coo%set_nrows(naggr)
call bcoo%set_ncols(naggr) call ac_coo%set_ncols(naggr)
call bcoo%set_dupl(psb_dupl_add_) call ac_coo%set_dupl(psb_dupl_add_)
call bcoo%fix(info) call ac_coo%fix(info)
call ac%mv_from(ac_coo)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then
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
else if (p%parms%coarse_mat == mld_distr_mat_) then
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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
info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='cscnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

185
mlprec/impl/mld_caggrmat_smth_asb.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_smth_asb use mld_c_inner_mod, mld_protect_name => mld_caggrmat_smth_asb
@ -103,9 +103,9 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_cspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -113,7 +113,7 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: am3, am4 type(psb_cspmat_type) :: am3, am4
type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_c_coo_sparse_mat) :: tmpcoo
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)
@ -138,14 +138,14 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -177,16 +177,16 @@ 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
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = cone tmpcoo%val(i) = cone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -217,19 +217,19 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= czero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= czero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -246,13 +246,13 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -260,11 +260,11 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -369,27 +369,27 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
call op_restr%mv_to(acoo) call op_restr%mv_to(tmpcoo)
nzl = acoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
! !
! 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(:)
! !
do k=1, nzl do k=1, nzl
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1 i = i+1
acoo%val(i) = acoo%val(k) tmpcoo%val(i) = tmpcoo%val(k)
acoo%ia(i) = acoo%ia(k) tmpcoo%ia(i) = tmpcoo%ia(k)
acoo%ja(i) = acoo%ja(k) tmpcoo%ja(i) = tmpcoo%ja(k)
end if end if
end do end do
call acoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call acoo%trim() call tmpcoo%trim()
call op_restr%mv_from(acoo) call op_restr%mv_from(tmpcoo)
call op_restr%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 op_restr') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
@ -413,113 +413,12 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3')
goto 9999
end if
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999
end if end if

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

@ -113,7 +113,11 @@ 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 type(psb_dspmat_type) :: ac, op_prol,op_restr
type(psb_d_coo_sparse_mat) :: acoo, bcoo
type(psb_d_csr_sparse_mat) :: acsr1
integer :: nzl,ntaggr
integer :: debug_level, debug_unit
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -121,6 +125,9 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(psb_get_errstatus().ne.0) return if(psb_get_errstatus().ne.0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context() ictxt = desc_a%get_context()
@ -129,42 +136,138 @@ 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_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,&
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999
end if
case(mld_smooth_prol_) case(mld_smooth_prol_)
call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_) case(mld_biz_prol_)
call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
case(mld_min_energy_)
call mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, &
& p%parms,ac,op_prol,op_restr,info)
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Invalid aggr kind')
goto 9999
end select
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb')
goto 9999
end if
ntaggr = sum(nlaggr)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = ac%get_nzeros()
call ac%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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
case(mld_min_energy_) 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 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='psb_sp_free')
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
end if 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
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 '
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_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999 goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

165
mlprec/impl/mld_daggrmat_biz_asb.f90
Unescape Escape View File

@ -78,7 +78,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_daggrmat_biz_asb use mld_d_inner_mod, mld_protect_name => mld_daggrmat_biz_asb
@ -88,9 +88,9 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -98,7 +98,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: am3, am4 type(psb_dspmat_type) :: am3, am4
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_d_coo_sparse_mat) :: tmpcoo
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)
@ -123,7 +123,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -131,7 +131,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -163,16 +163,16 @@ subroutine mld_daggrmat_biz_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
call acoo%allocate(ncol,naggr,ncol) call tmpcoo%allocate(ncol,naggr,ncol)
do i=1,nrow do i=1,nrow
acoo%val(i) = done tmpcoo%val(i) = done
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(nrow) call tmpcoo%set_nzeros(nrow)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -203,19 +203,19 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= dzero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= dzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -232,9 +232,9 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
! !
! This only works with CSR ! This only works with CSR
@ -261,7 +261,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -269,11 +269,11 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -372,7 +372,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -389,10 +389,10 @@ subroutine mld_daggrmat_biz_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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
@ -401,107 +401,6 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

184
mlprec/impl/mld_daggrmat_minnrg_asb.f90
Unescape Escape View File

@ -98,7 +98,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_daggrmat_minnrg_asb use mld_d_inner_mod, mld_protect_name => mld_daggrmat_minnrg_asb
@ -108,9 +108,9 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
@ -121,8 +121,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_dspmat_type) :: 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) :: tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
real(psb_dpk_), allocatable :: adiag(:), adinv(:) real(psb_dpk_), allocatable :: adiag(:), adinv(:)
real(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:) real(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
@ -150,7 +150,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -165,7 +165,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -207,16 +207,16 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = done tmpcoo%val(i) = done
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acoo%set_asb() call tmpcoo%set_asb()
call ptilde%mv_from(acoo) call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr') call ptilde%cscnv(info,type='csr')
!!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1')
@ -570,164 +570,18 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then 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 = op_restr x am3') &a_err='Build ac = op_restr x am3')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done mv_to_coo'
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' B matrix nzl: ',nzl
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Assembld aux descr. distr.'
call bcoo%set_nrows(p%desc_ac%get_local_rows())
call bcoo%set_ncols(p%desc_ac%get_local_cols())
call bcoo%fix(info)
call p%ac%mv_from(bcoo)
call p%ac%set_asb()
call p%ac%cscnv(info,type='csr')
if (np>1) then
call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%mv_to(tmpcoo)
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_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr')
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 '
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
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
!!$ nzbr(:) = 0
!!$ 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
!!$ 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
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => R i.e. restriction operator
! op_prol => P i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

144
mlprec/impl/mld_daggrmat_nosmth_asb.f90
Unescape Escape View File

@ -81,7 +81,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_daggrmat_nosmth_asb use mld_d_inner_mod, mld_protect_name => mld_daggrmat_nosmth_asb
@ -91,16 +91,16 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
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
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo type(psb_d_coo_sparse_mat) :: ac_coo, acoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2 type(psb_d_csr_sparse_mat) :: acsr1, acsr2
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
@ -134,136 +134,36 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call acoo1%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,nrow do i=1,nrow
acoo1%val(i) = done acoo%val(i) = done
acoo1%ia(i) = i acoo%ia(i) = i
acoo1%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo1%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call acoo1%set_nzeros(nrow) call acoo%set_nzeros(nrow)
call acoo1%set_asb() call acoo%set_asb()
call acoo1%fix(info) call acoo%fix(info)
call op_prol%mv_from(acoo1) call op_prol%mv_from(acoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr) if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow) call a%csclip(ac_coo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = ac_coo%get_nzeros()
do i=1, nzt do i=1, nzt
bcoo%ia(i) = ilaggr(bcoo%ia(i)) ac_coo%ia(i) = ilaggr(ac_coo%ia(i))
bcoo%ja(i) = ilaggr(bcoo%ja(i)) ac_coo%ja(i) = ilaggr(ac_coo%ja(i))
enddo enddo
call bcoo%set_nrows(naggr) call ac_coo%set_nrows(naggr)
call bcoo%set_ncols(naggr) call ac_coo%set_ncols(naggr)
call bcoo%set_dupl(psb_dupl_add_) call ac_coo%set_dupl(psb_dupl_add_)
call bcoo%fix(info) call ac_coo%fix(info)
call ac%mv_from(ac_coo)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then
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
else if (p%parms%coarse_mat == mld_distr_mat_) then
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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
info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='cscnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

185
mlprec/impl/mld_daggrmat_smth_asb.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_daggrmat_smth_asb use mld_d_inner_mod, mld_protect_name => mld_daggrmat_smth_asb
@ -103,9 +103,9 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_dspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -113,7 +113,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: am3, am4 type(psb_dspmat_type) :: am3, am4
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_d_coo_sparse_mat) :: tmpcoo
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)
@ -138,14 +138,14 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -177,16 +177,16 @@ 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
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = done tmpcoo%val(i) = done
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -217,19 +217,19 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= dzero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= dzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -246,13 +246,13 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -260,11 +260,11 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -369,27 +369,27 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
call op_restr%mv_to(acoo) call op_restr%mv_to(tmpcoo)
nzl = acoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
! !
! 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(:)
! !
do k=1, nzl do k=1, nzl
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1 i = i+1
acoo%val(i) = acoo%val(k) tmpcoo%val(i) = tmpcoo%val(k)
acoo%ia(i) = acoo%ia(k) tmpcoo%ia(i) = tmpcoo%ia(k)
acoo%ja(i) = acoo%ja(k) tmpcoo%ja(i) = tmpcoo%ja(k)
end if end if
end do end do
call acoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call acoo%trim() call tmpcoo%trim()
call op_restr%mv_from(acoo) call op_restr%mv_from(tmpcoo)
call op_restr%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 op_restr') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
@ -413,113 +413,12 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3')
goto 9999
end if
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999
end if end if

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

@ -113,7 +113,11 @@ 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 type(psb_sspmat_type) :: ac, op_prol,op_restr
type(psb_s_coo_sparse_mat) :: acoo, bcoo
type(psb_s_csr_sparse_mat) :: acsr1
integer :: nzl,ntaggr
integer :: debug_level, debug_unit
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -121,6 +125,9 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(psb_get_errstatus().ne.0) return if(psb_get_errstatus().ne.0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context() ictxt = desc_a%get_context()
@ -129,42 +136,138 @@ 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_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,&
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999
end if
case(mld_smooth_prol_) case(mld_smooth_prol_)
call mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_) case(mld_biz_prol_)
call mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
case(mld_min_energy_)
call mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, &
& p%parms,ac,op_prol,op_restr,info)
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Invalid aggr kind')
goto 9999
end select
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb')
goto 9999
end if
ntaggr = sum(nlaggr)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = ac%get_nzeros()
call ac%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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
case(mld_min_energy_) 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 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='psb_sp_free')
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
end if 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
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 '
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_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999 goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

165
mlprec/impl/mld_saggrmat_biz_asb.f90
Unescape Escape View File

@ -78,7 +78,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_biz_asb use mld_s_inner_mod, mld_protect_name => mld_saggrmat_biz_asb
@ -88,9 +88,9 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_sspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -98,7 +98,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: am3, am4 type(psb_sspmat_type) :: am3, am4
type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_s_coo_sparse_mat) :: tmpcoo
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)
@ -123,7 +123,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -131,7 +131,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -163,16 +163,16 @@ subroutine mld_saggrmat_biz_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
call acoo%allocate(ncol,naggr,ncol) call tmpcoo%allocate(ncol,naggr,ncol)
do i=1,nrow do i=1,nrow
acoo%val(i) = sone tmpcoo%val(i) = sone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(nrow) call tmpcoo%set_nzeros(nrow)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -203,19 +203,19 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= szero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= szero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -232,9 +232,9 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
! !
! This only works with CSR ! This only works with CSR
@ -261,7 +261,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -269,11 +269,11 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -372,7 +372,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -389,10 +389,10 @@ subroutine mld_saggrmat_biz_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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
@ -401,107 +401,6 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

184
mlprec/impl/mld_saggrmat_minnrg_asb.f90
Unescape Escape View File

@ -98,7 +98,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_minnrg_asb use mld_s_inner_mod, mld_protect_name => mld_saggrmat_minnrg_asb
@ -108,9 +108,9 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_sspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
@ -121,8 +121,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_sspmat_type) :: 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) :: tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
real(psb_spk_), allocatable :: adiag(:), adinv(:) real(psb_spk_), allocatable :: adiag(:), adinv(:)
real(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:) real(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
@ -150,7 +150,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -165,7 +165,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -207,16 +207,16 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = sone tmpcoo%val(i) = sone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acoo%set_asb() call tmpcoo%set_asb()
call ptilde%mv_from(acoo) call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr') call ptilde%cscnv(info,type='csr')
!!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1')
@ -570,164 +570,18 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then 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 = op_restr x am3') &a_err='Build ac = op_restr x am3')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done mv_to_coo'
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' B matrix nzl: ',nzl
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Assembld aux descr. distr.'
call bcoo%set_nrows(p%desc_ac%get_local_rows())
call bcoo%set_ncols(p%desc_ac%get_local_cols())
call bcoo%fix(info)
call p%ac%mv_from(bcoo)
call p%ac%set_asb()
call p%ac%cscnv(info,type='csr')
if (np>1) then
call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%mv_to(tmpcoo)
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_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr')
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 '
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
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
!!$ nzbr(:) = 0
!!$ 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
!!$ 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
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => R i.e. restriction operator
! op_prol => P i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

144
mlprec/impl/mld_saggrmat_nosmth_asb.f90
Unescape Escape View File

@ -81,7 +81,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_nosmth_asb use mld_s_inner_mod, mld_protect_name => mld_saggrmat_nosmth_asb
@ -91,16 +91,16 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
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
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_s_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo type(psb_s_coo_sparse_mat) :: ac_coo, acoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2 type(psb_s_csr_sparse_mat) :: acsr1, acsr2
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
@ -134,136 +134,36 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call acoo1%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,nrow do i=1,nrow
acoo1%val(i) = sone acoo%val(i) = sone
acoo1%ia(i) = i acoo%ia(i) = i
acoo1%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo1%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call acoo1%set_nzeros(nrow) call acoo%set_nzeros(nrow)
call acoo1%set_asb() call acoo%set_asb()
call acoo1%fix(info) call acoo%fix(info)
call op_prol%mv_from(acoo1) call op_prol%mv_from(acoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr) if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow) call a%csclip(ac_coo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = ac_coo%get_nzeros()
do i=1, nzt do i=1, nzt
bcoo%ia(i) = ilaggr(bcoo%ia(i)) ac_coo%ia(i) = ilaggr(ac_coo%ia(i))
bcoo%ja(i) = ilaggr(bcoo%ja(i)) ac_coo%ja(i) = ilaggr(ac_coo%ja(i))
enddo enddo
call bcoo%set_nrows(naggr) call ac_coo%set_nrows(naggr)
call bcoo%set_ncols(naggr) call ac_coo%set_ncols(naggr)
call bcoo%set_dupl(psb_dupl_add_) call ac_coo%set_dupl(psb_dupl_add_)
call bcoo%fix(info) call ac_coo%fix(info)
call ac%mv_from(ac_coo)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then
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
else if (p%parms%coarse_mat == mld_distr_mat_) then
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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
info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='cscnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

185
mlprec/impl/mld_saggrmat_smth_asb.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_smth_asb use mld_s_inner_mod, mld_protect_name => mld_saggrmat_smth_asb
@ -103,9 +103,9 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_sspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -113,7 +113,7 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: am3, am4 type(psb_sspmat_type) :: am3, am4
type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_s_coo_sparse_mat) :: tmpcoo
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)
@ -138,14 +138,14 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -177,16 +177,16 @@ 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
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = sone tmpcoo%val(i) = sone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -217,19 +217,19 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= szero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= szero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -246,13 +246,13 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -260,11 +260,11 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -369,27 +369,27 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
call op_restr%mv_to(acoo) call op_restr%mv_to(tmpcoo)
nzl = acoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
! !
! 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(:)
! !
do k=1, nzl do k=1, nzl
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1 i = i+1
acoo%val(i) = acoo%val(k) tmpcoo%val(i) = tmpcoo%val(k)
acoo%ia(i) = acoo%ia(k) tmpcoo%ia(i) = tmpcoo%ia(k)
acoo%ja(i) = acoo%ja(k) tmpcoo%ja(i) = tmpcoo%ja(k)
end if end if
end do end do
call acoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call acoo%trim() call tmpcoo%trim()
call op_restr%mv_from(acoo) call op_restr%mv_from(tmpcoo)
call op_restr%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 op_restr') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
@ -413,113 +413,12 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3')
goto 9999
end if
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999
end if end if

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

@ -113,7 +113,11 @@ 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 type(psb_zspmat_type) :: ac, op_prol,op_restr
type(psb_z_coo_sparse_mat) :: acoo, bcoo
type(psb_z_csr_sparse_mat) :: acsr1
integer :: nzl,ntaggr
integer :: debug_level, debug_unit
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
character(len=20) :: name character(len=20) :: name
@ -121,6 +125,9 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
if(psb_get_errstatus().ne.0) return if(psb_get_errstatus().ne.0) return
info=psb_success_ info=psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context() ictxt = desc_a%get_context()
@ -129,42 +136,138 @@ 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_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,&
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_nosmth_asb')
goto 9999
end if
case(mld_smooth_prol_) case(mld_smooth_prol_)
call mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999
end if
case(mld_biz_prol_) case(mld_biz_prol_)
call mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) call mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr, &
if(info /= psb_success_) then & p%parms,ac,op_prol,op_restr,info)
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
case(mld_min_energy_)
call mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr, &
& p%parms,ac,op_prol,op_restr,info)
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Invalid aggr kind')
goto 9999
end select
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat asb')
goto 9999
end if
ntaggr = sum(nlaggr)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = ac%get_nzeros()
call ac%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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
case(mld_min_energy_) 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 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='psb_sp_free')
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_smth_asb')
goto 9999 goto 9999
end if 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
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 '
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_) &
& call psb_gather(p%ac,ac,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid aggr kind') call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999 goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

165
mlprec/impl/mld_zaggrmat_biz_asb.f90
Unescape Escape View File

@ -78,7 +78,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_biz_asb use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_biz_asb
@ -88,9 +88,9 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -98,7 +98,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: am3, am4 type(psb_zspmat_type) :: am3, am4
type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_z_coo_sparse_mat) :: tmpcoo
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)
@ -123,7 +123,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -131,7 +131,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -163,16 +163,16 @@ subroutine mld_zaggrmat_biz_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
call acoo%allocate(ncol,naggr,ncol) call tmpcoo%allocate(ncol,naggr,ncol)
do i=1,nrow do i=1,nrow
acoo%val(i) = zone tmpcoo%val(i) = zone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(nrow) call tmpcoo%set_nzeros(nrow)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -203,19 +203,19 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= zzero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= zzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -232,9 +232,9 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
! !
! This only works with CSR ! This only works with CSR
@ -261,7 +261,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -269,11 +269,11 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -372,7 +372,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -389,10 +389,10 @@ subroutine mld_zaggrmat_biz_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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999 goto 9999
@ -401,107 +401,6 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,p,info)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

184
mlprec/impl/mld_zaggrmat_minnrg_asb.f90
Unescape Escape View File

@ -98,7 +98,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_minnrg_asb use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_minnrg_asb
@ -108,9 +108,9 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: op_prol,op_restr, b
! Local variables ! Local variables
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
@ -121,8 +121,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_zspmat_type) :: 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) :: tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
complex(psb_dpk_), allocatable :: adiag(:), adinv(:) complex(psb_dpk_), allocatable :: adiag(:), adinv(:)
complex(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:) complex(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
@ -150,7 +150,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
@ -165,7 +165,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -207,16 +207,16 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form ! 1. Allocate Ptilde in sparse matrix form
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = zone tmpcoo%val(i) = zone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acoo%set_asb() call tmpcoo%set_asb()
call ptilde%mv_from(acoo) call ptilde%mv_from(tmpcoo)
call ptilde%cscnv(info,type='csr') call ptilde%cscnv(info,type='csr')
!!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1') !!$ call local_dump(me,ptilde,'csr-ptilde','Ptilde-1')
@ -570,164 +570,18 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
if (info == psb_success_) call am3%free() if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
if (info /= psb_success_) then 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 = op_restr x am3') &a_err='Build ac = op_restr x am3')
goto 9999 goto 9999
end if end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done mv_to_coo'
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
call b%mv_to(bcoo)
nzl = bcoo%get_nzeros()
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' B matrix nzl: ',nzl
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Assembld aux descr. distr.'
call bcoo%set_nrows(p%desc_ac%get_local_rows())
call bcoo%set_ncols(p%desc_ac%get_local_cols())
call bcoo%fix(info)
call p%ac%mv_from(bcoo)
call p%ac%set_asb()
call p%ac%cscnv(info,type='csr')
if (np>1) then
call op_prol%mv_to(acsr)
nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call op_prol%mv_from(acsr)
endif
call op_prol%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call op_restr%mv_to(tmpcoo)
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_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting op_restr to local')
goto 9999
end if
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr')
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 '
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
!!$ call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
!!$ nzbr(:) = 0
!!$ 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
!!$ 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
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => R i.e. restriction operator
! op_prol => P i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&

144
mlprec/impl/mld_zaggrmat_nosmth_asb.f90
Unescape Escape View File

@ -81,7 +81,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_nosmth_asb use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_nosmth_asb
@ -91,16 +91,16 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr
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
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
type(psb_z_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo, acoo type(psb_z_coo_sparse_mat) :: ac_coo, acoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2 type(psb_z_csr_sparse_mat) :: acsr1, acsr2
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
integer :: nrow, nglob, ncol, ntaggr, nzl, ip, & integer :: nrow, nglob, ncol, ntaggr, nzl, ip, &
@ -134,136 +134,36 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
call acoo1%allocate(ncol,ntaggr,ncol) call acoo%allocate(ncol,ntaggr,ncol)
do i=1,nrow do i=1,nrow
acoo1%val(i) = zone acoo%val(i) = zone
acoo1%ia(i) = i acoo%ia(i) = i
acoo1%ja(i) = ilaggr(i) acoo%ja(i) = ilaggr(i)
end do end do
call acoo1%set_dupl(psb_dupl_add_) call acoo%set_dupl(psb_dupl_add_)
call acoo1%set_nzeros(nrow) call acoo%set_nzeros(nrow)
call acoo1%set_asb() call acoo%set_asb()
call acoo1%fix(info) call acoo%fix(info)
call op_prol%mv_from(acoo1) call op_prol%mv_from(acoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_) call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) call op_prol%transp(op_restr) if (info == psb_success_) call op_prol%transp(op_restr)
call a%csclip(bcoo,info,jmax=nrow) call a%csclip(ac_coo,info,jmax=nrow)
nzt = bcoo%get_nzeros() nzt = ac_coo%get_nzeros()
do i=1, nzt do i=1, nzt
bcoo%ia(i) = ilaggr(bcoo%ia(i)) ac_coo%ia(i) = ilaggr(ac_coo%ia(i))
bcoo%ja(i) = ilaggr(bcoo%ja(i)) ac_coo%ja(i) = ilaggr(ac_coo%ja(i))
enddo enddo
call bcoo%set_nrows(naggr) call ac_coo%set_nrows(naggr)
call bcoo%set_ncols(naggr) call ac_coo%set_ncols(naggr)
call bcoo%set_dupl(psb_dupl_add_) call ac_coo%set_dupl(psb_dupl_add_)
call bcoo%fix(info) call ac_coo%fix(info)
call ac%mv_from(ac_coo)
call b%mv_from(bcoo)
if (p%parms%coarse_mat == mld_repl_mat_) then
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
else if (p%parms%coarse_mat == mld_distr_mat_) then
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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
info = psb_err_internal_error_
call psb_errpush(psb_err_internal_error_,name,a_err='invalid mld_coarse_mat_')
goto 9999
end if
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='cscnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
if (info == psb_success_) &
& p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='linmap build')
goto 9999
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

185
mlprec/impl/mld_zaggrmat_smth_asb.f90
Unescape Escape View File

@ -93,7 +93,7 @@
! info - integer, output. ! info - integer, output.
! Error code. ! Error code.
! !
subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,info)
use psb_base_mod use psb_base_mod
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_smth_asb use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_smth_asb
@ -103,9 +103,9 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
type(psb_zspmat_type) :: op_prol, op_restr, b
! Local variables ! Local variables
integer :: nrow, nglob, ncol, ntaggr, ip, ndx,& integer :: nrow, nglob, ncol, ntaggr, ip, ndx,&
@ -113,7 +113,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt, np, me, err_act integer ::ictxt, np, me, err_act
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: am3, am4 type(psb_zspmat_type) :: am3, am4
type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo type(psb_z_coo_sparse_mat) :: tmpcoo
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)
@ -138,14 +138,14 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
nrow = desc_a%get_local_rows() nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols() ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh theta = parms%aggr_thresh
naggr = nlaggr(me+1) naggr = nlaggr(me+1)
ntaggr = sum(nlaggr) ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me)) naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1)) naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_) filter_mat = (parms%aggr_filter == mld_filter_mat_)
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)
@ -177,16 +177,16 @@ 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
call acoo%allocate(ncol,ntaggr,ncol) call tmpcoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol do i=1,ncol
acoo%val(i) = zone tmpcoo%val(i) = zone
acoo%ia(i) = i tmpcoo%ia(i) = i
acoo%ja(i) = ilaggr(i) tmpcoo%ja(i) = ilaggr(i)
end do end do
call acoo%set_nzeros(ncol) call tmpcoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info) call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_) if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
@ -217,19 +217,19 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if end if
enddo enddo
! Take out zeroed terms ! Take out zeroed terms
call acsrf%mv_to_coo(acoof,info) call acsrf%mv_to_coo(tmpcoo,info)
k = 0 k = 0
do j=1,acoof%get_nzeros() do j=1,tmpcoo%get_nzeros()
if ((acoof%val(j) /= zzero) .or. (acoof%ia(j) == acoof%ja(j))) then if ((tmpcoo%val(j) /= zzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1 k = k + 1
acoof%val(k) = acoof%val(j) tmpcoo%val(k) = tmpcoo%val(j)
acoof%ia(k) = acoof%ia(j) tmpcoo%ia(k) = tmpcoo%ia(j)
acoof%ja(k) = acoof%ja(j) tmpcoo%ja(k) = tmpcoo%ja(j)
end if end if
end do end do
call acoof%set_nzeros(k) call tmpcoo%set_nzeros(k)
call acoof%set_dupl(psb_dupl_add_) call tmpcoo%set_dupl(psb_dupl_add_)
call acsrf%mv_from_coo(acoof,info) call acsrf%mv_from_coo(tmpcoo,info)
end if end if
@ -246,13 +246,13 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info /= psb_success_) goto 9999 if (info /= psb_success_) goto 9999
if (p%parms%aggr_omega_alg == mld_eig_est_) then if (parms%aggr_omega_alg == mld_eig_est_) then
if (p%parms%aggr_eig == mld_max_norm_) then if (parms%aggr_eig == mld_max_norm_) then
anorm = acsr3%csnmi() anorm = acsr3%csnmi()
omega = 4.d0/(3.d0*anorm) omega = 4.d0/(3.d0*anorm)
p%parms%aggr_omega_val = omega parms%aggr_omega_val = omega
else else
info = psb_err_internal_error_ info = psb_err_internal_error_
@ -260,11 +260,11 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
goto 9999 goto 9999
end if end if
else if (p%parms%aggr_omega_alg == mld_user_choice_) then else if (parms%aggr_omega_alg == mld_user_choice_) then
omega = p%parms%aggr_omega_val omega = parms%aggr_omega_val
else if (p%parms%aggr_omega_alg /= mld_user_choice_) then else if (parms%aggr_omega_alg /= mld_user_choice_) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_') call psb_errpush(info,name,a_err='invalid mld_aggr_omega_alg_')
goto 9999 goto 9999
@ -369,27 +369,27 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_numbmm(a,op_prol,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',parms%aggr_kind, mld_smooth_prol_
call op_prol%transp(op_restr) call op_prol%transp(op_restr)
call op_restr%mv_to(acoo) call op_restr%mv_to(tmpcoo)
nzl = acoo%get_nzeros() nzl = tmpcoo%get_nzeros()
i=0 i=0
! !
! 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(:)
! !
do k=1, nzl do k=1, nzl
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1 i = i+1
acoo%val(i) = acoo%val(k) tmpcoo%val(i) = tmpcoo%val(k)
acoo%ia(i) = acoo%ia(k) tmpcoo%ia(i) = tmpcoo%ia(k)
acoo%ja(i) = acoo%ja(k) tmpcoo%ja(i) = tmpcoo%ja(k)
end if end if
end do end do
call acoo%set_nzeros(i) call tmpcoo%set_nzeros(i)
call acoo%trim() call tmpcoo%trim()
call op_restr%mv_from(acoo) call op_restr%mv_from(tmpcoo)
call op_restr%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 op_restr') call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv op_restr')
@ -413,113 +413,12 @@ 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(op_restr,am3,b,info) call psb_symbmm(op_restr,am3,ac,info)
if (info == psb_success_) call psb_numbmm(op_restr,am3,b) if (info == psb_success_) call psb_numbmm(op_restr,am3,ac)
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 ac%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 = op_restr x am3') call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = op_restr x am3')
goto 9999
end if
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = b%get_nzeros()
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_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call p%ac%mv_from(bcoo)
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()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_sp_free')
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
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 '
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_) &
& call psb_gather(p%ac,b,p%desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! op_restr => PR^T i.e. restriction operator
! op_prol => PR i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,op_restr,op_prol,ilaggr,nlaggr)
if (info == psb_success_) call op_prol%free()
if (info == psb_success_) call op_restr%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999 goto 9999
end if end if

7
mlprec/mld_c_inner_mod.f90
Unescape Escape View File

@ -146,13 +146,14 @@ module mld_c_inner_mod
abstract interface abstract interface
subroutine mld_caggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_caggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,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, mld_sml_parms
type(psb_cspmat_type), intent(in) :: a type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_c_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_cspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_caggrmat_var_asb end subroutine mld_caggrmat_var_asb
end interface end interface

7
mlprec/mld_d_inner_mod.f90
Unescape Escape View File

@ -146,13 +146,14 @@ module mld_d_inner_mod
abstract interface abstract interface
subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_daggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,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, mld_dml_parms
type(psb_dspmat_type), intent(in) :: a type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_d_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_dspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_daggrmat_var_asb end subroutine mld_daggrmat_var_asb
end interface end interface

7
mlprec/mld_s_inner_mod.f90
Unescape Escape View File

@ -146,13 +146,14 @@ module mld_s_inner_mod
abstract interface abstract interface
subroutine mld_saggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_saggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,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, mld_sml_parms
type(psb_sspmat_type), intent(in) :: a type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_s_onelev_type), intent(inout), target :: p type(mld_sml_parms), intent(inout) :: parms
type(psb_sspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_saggrmat_var_asb end subroutine mld_saggrmat_var_asb
end interface end interface

7
mlprec/mld_z_inner_mod.f90
Unescape Escape View File

@ -146,13 +146,14 @@ module mld_z_inner_mod
abstract interface abstract interface
subroutine mld_zaggrmat_var_asb(a,desc_a,ilaggr,nlaggr,p,info) subroutine mld_zaggrmat_var_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr,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, mld_dml_parms
type(psb_zspmat_type), intent(in) :: a type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:) integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_z_onelev_type), intent(inout), target :: p type(mld_dml_parms), intent(inout) :: parms
type(psb_zspmat_type), intent(out) :: ac,op_prol,op_restr
integer, intent(out) :: info integer, intent(out) :: info
end subroutine mld_zaggrmat_var_asb end subroutine mld_zaggrmat_var_asb
end interface end interface

Write Preview
Loading…
Cancel
Save