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

Rework unsmoothed and "bizarre" bld to align with smoothed.

Browse Source
stopcriterion
Salvatore Filippone 5 years ago
parent 85a03d87ba
commit b339d867da
8 changed files with 640 additions and 828 deletions
Show Stats Download Patch File Download Diff File

251
mlprec/impl/aggregator/mld_caggrmat_biz_bld.f90
Unescape Escape View File

@ -84,34 +84,40 @@ subroutine mld_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
use psb_base_mod
use mld_base_prec_type
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_biz_bld
use mld_c_base_aggregator_mod
implicit none
! Arguments
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lcspmat_type), intent(inout) :: op_prol
type(psb_lcspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
integer(psb_ipk_) :: inaggr, nzlp, inrow
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_lcspmat_type) :: am3, am4,tmp_prol, la
type(psb_lc_coo_sparse_mat) :: tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_desc_type) :: tmp_desc
type(psb_c_coo_sparse_mat) :: icoo
type(psb_lc_coo_sparse_mat) :: coo_prol, coo_restr, tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr
complex(psb_spk_), allocatable :: adiag(:)
real(psb_spk_), allocatable :: arwsum(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta
logical, parameter :: debug_new=.false.
character(len=80) :: filename
name='mld_aggrmat_biz_bld'
name='mld_aggrmat_smth_bld'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
@ -120,7 +126,6 @@ subroutine mld_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -130,40 +135,45 @@ subroutine mld_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (parms%aggr_filter == mld_filter_mat_)
!
! naggr: number of local aggregates
! nrow: local rows.
!
! Get the diagonal D
adiag = a%get_diag(info)
adiag = a%get_diag(info)
if (info == psb_success_) &
& call psb_realloc(ncol,adiag,info)
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
if (info == psb_success_) call a%csclip(icoo,info,jmax=a%get_nrows())
call icoo%mv_to_lfmt(acsr,info)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies sone.'
& ' Initial copies done.'
call acsr%cp_to_fmt(acsrf,info)
if (filter_mat) then
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1, nrow
tmp = czero
@ -195,39 +205,13 @@ subroutine mld_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
end if
end do
if (filter_mat) call acsrf%scal(adiag,info)
if (info == psb_success_) call acsr3%scal(adiag,info)
if (info /= psb_success_) goto 9999
if (parms%aggr_omega_alg == mld_eig_est_) then
if (parms%aggr_eig == mld_max_norm_) then
!
! This only works with CSR
!
anorm = szero
dg = sone
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
call psb_amx(ictxt,anorm)
omega = 4.d0/(3.d0*anorm)
parms%aggr_omega_val = omega
@ -247,118 +231,51 @@ subroutine mld_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (filter_mat) then
!
! Build the smoothed prolongator using the filtered matrix
!
do i=1,acsrf%get_nrows()
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) then
acsrf%val(j) = cone - omega*acsrf%val(j)
else
acsrf%val(j) = - omega*acsrf%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
!
do i=1,acsr3%get_nrows()
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) == i) then
acsr3%val(j) = cone - omega*acsr3%val(j)
else
acsr3%val(j) = - omega*acsr3%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*A)Ptilde
! Doing it this way means to consider diag(A_i)
!
!
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
end if
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call acsrf%scal(adiag,info)
if (info /= psb_success_) goto 9999
call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call coo_prol%mv_to_fmt(csr_prol,info)
!
! Build the smoothed prolongator using either A or Af
! acsr1 = (I-w*D*A) Prol acsr1 = (I-w*D*Af) Prol
! This is always done through the variable acsrf which
! is a bit less readable, butsaves space and one extra matrix copy
!
call omega_smooth(omega,acsrf)
call psb_par_spspmm(acsrf,desc_a,csr_prol,acsr1,tmp_desc,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call op_prol%mv_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -368,8 +285,26 @@ subroutine mld_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
9999 continue
call psb_errpush(info,name)
call psb_error_handler(err_act)
return
contains
subroutine omega_smooth(omega,acsr)
implicit none
real(psb_spk_),intent(in) :: omega
type(psb_lc_csr_sparse_mat), intent(inout) :: acsr
!
integer(psb_lpk_) :: i,j
do i=1,acsr%get_nrows()
do j=acsr%irp(i),acsr%irp(i+1)-1
if (acsr%ja(j) == i) then
acsr%val(j) = cone - omega*acsr%val(j)
else
acsr%val(j) = - omega*acsr%val(j)
end if
end do
end do
end subroutine omega_smooth
end subroutine mld_caggrmat_biz_bld

116
mlprec/impl/aggregator/mld_caggrmat_nosmth_bld.f90
Unescape Escape View File

@ -100,7 +100,7 @@ subroutine mld_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
use psb_base_mod
use mld_base_prec_type
use mld_c_inner_mod, mld_protect_name => mld_caggrmat_nosmth_bld
use mld_c_base_aggregator_mod
implicit none
! Arguments
@ -118,11 +118,13 @@ subroutine mld_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_lcspmat_type) :: la
type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2
type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo, coo_prol, coo_restr
type(psb_lc_csr_sparse_mat) :: acsr1, acsr2, acsr
type(psb_desc_type) :: tmp_desc
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
integer(psb_ipk_) :: inaggr, nzlp
name = 'mld_aggrmat_nosmth_bld'
info = psb_success_
@ -143,53 +145,69 @@ subroutine mld_caggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
if (.false.) then
call a%cp_to(acsr)
call op_prol%mv_to(coo_prol)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
end if
call psb_erractionrestore(err_act)
return

251
mlprec/impl/aggregator/mld_daggrmat_biz_bld.f90
Unescape Escape View File

@ -84,34 +84,40 @@ subroutine mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
use psb_base_mod
use mld_base_prec_type
use mld_d_inner_mod, mld_protect_name => mld_daggrmat_biz_bld
use mld_d_base_aggregator_mod
implicit none
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ldspmat_type), intent(inout) :: op_prol
type(psb_ldspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
integer(psb_ipk_) :: inaggr, nzlp, inrow
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_ldspmat_type) :: am3, am4,tmp_prol, la
type(psb_ld_coo_sparse_mat) :: tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_desc_type) :: tmp_desc
type(psb_d_coo_sparse_mat) :: icoo
type(psb_ld_coo_sparse_mat) :: coo_prol, coo_restr, tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr
real(psb_dpk_), allocatable :: adiag(:)
real(psb_dpk_), allocatable :: arwsum(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
logical, parameter :: debug_new=.false.
character(len=80) :: filename
name='mld_aggrmat_biz_bld'
name='mld_aggrmat_smth_bld'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
@ -120,7 +126,6 @@ subroutine mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -130,40 +135,45 @@ subroutine mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (parms%aggr_filter == mld_filter_mat_)
!
! naggr: number of local aggregates
! nrow: local rows.
!
! Get the diagonal D
adiag = a%get_diag(info)
adiag = a%get_diag(info)
if (info == psb_success_) &
& call psb_realloc(ncol,adiag,info)
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
if (info == psb_success_) call a%csclip(icoo,info,jmax=a%get_nrows())
call icoo%mv_to_lfmt(acsr,info)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies sone.'
& ' Initial copies done.'
call acsr%cp_to_fmt(acsrf,info)
if (filter_mat) then
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1, nrow
tmp = dzero
@ -195,39 +205,13 @@ subroutine mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
end if
end do
if (filter_mat) call acsrf%scal(adiag,info)
if (info == psb_success_) call acsr3%scal(adiag,info)
if (info /= psb_success_) goto 9999
if (parms%aggr_omega_alg == mld_eig_est_) then
if (parms%aggr_eig == mld_max_norm_) then
!
! This only works with CSR
!
anorm = dzero
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = dzero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
call psb_amx(ictxt,anorm)
omega = 4.d0/(3.d0*anorm)
parms%aggr_omega_val = omega
@ -247,118 +231,51 @@ subroutine mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (filter_mat) then
!
! Build the smoothed prolongator using the filtered matrix
!
do i=1,acsrf%get_nrows()
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) then
acsrf%val(j) = done - omega*acsrf%val(j)
else
acsrf%val(j) = - omega*acsrf%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
!
do i=1,acsr3%get_nrows()
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) == i) then
acsr3%val(j) = done - omega*acsr3%val(j)
else
acsr3%val(j) = - omega*acsr3%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*A)Ptilde
! Doing it this way means to consider diag(A_i)
!
!
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
end if
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call acsrf%scal(adiag,info)
if (info /= psb_success_) goto 9999
call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call coo_prol%mv_to_fmt(csr_prol,info)
!
! Build the smoothed prolongator using either A or Af
! acsr1 = (I-w*D*A) Prol acsr1 = (I-w*D*Af) Prol
! This is always done through the variable acsrf which
! is a bit less readable, butsaves space and one extra matrix copy
!
call omega_smooth(omega,acsrf)
call psb_par_spspmm(acsrf,desc_a,csr_prol,acsr1,tmp_desc,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call op_prol%mv_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -368,8 +285,26 @@ subroutine mld_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
9999 continue
call psb_errpush(info,name)
call psb_error_handler(err_act)
return
contains
subroutine omega_smooth(omega,acsr)
implicit none
real(psb_dpk_),intent(in) :: omega
type(psb_ld_csr_sparse_mat), intent(inout) :: acsr
!
integer(psb_lpk_) :: i,j
do i=1,acsr%get_nrows()
do j=acsr%irp(i),acsr%irp(i+1)-1
if (acsr%ja(j) == i) then
acsr%val(j) = done - omega*acsr%val(j)
else
acsr%val(j) = - omega*acsr%val(j)
end if
end do
end do
end subroutine omega_smooth
end subroutine mld_daggrmat_biz_bld

116
mlprec/impl/aggregator/mld_daggrmat_nosmth_bld.f90
Unescape Escape View File

@ -100,7 +100,7 @@ subroutine mld_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
use psb_base_mod
use mld_base_prec_type
use mld_d_inner_mod, mld_protect_name => mld_daggrmat_nosmth_bld
use mld_d_base_aggregator_mod
implicit none
! Arguments
@ -118,11 +118,13 @@ subroutine mld_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_ldspmat_type) :: la
type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2
type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo, coo_prol, coo_restr
type(psb_ld_csr_sparse_mat) :: acsr1, acsr2, acsr
type(psb_desc_type) :: tmp_desc
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
integer(psb_ipk_) :: inaggr, nzlp
name = 'mld_aggrmat_nosmth_bld'
info = psb_success_
@ -143,53 +145,69 @@ subroutine mld_daggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
if (.false.) then
call a%cp_to(acsr)
call op_prol%mv_to(coo_prol)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
end if
call psb_erractionrestore(err_act)
return

251
mlprec/impl/aggregator/mld_saggrmat_biz_bld.f90
Unescape Escape View File

@ -84,34 +84,40 @@ subroutine mld_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
use psb_base_mod
use mld_base_prec_type
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_biz_bld
use mld_s_base_aggregator_mod
implicit none
! Arguments
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_sspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lsspmat_type), intent(inout) :: op_prol
type(psb_lsspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
integer(psb_ipk_) :: inaggr, nzlp, inrow
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_lsspmat_type) :: am3, am4,tmp_prol, la
type(psb_ls_coo_sparse_mat) :: tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_desc_type) :: tmp_desc
type(psb_s_coo_sparse_mat) :: icoo
type(psb_ls_coo_sparse_mat) :: coo_prol, coo_restr, tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr
real(psb_spk_), allocatable :: adiag(:)
real(psb_spk_), allocatable :: arwsum(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_spk_) :: anorm, omega, tmp, dg, theta
logical, parameter :: debug_new=.false.
character(len=80) :: filename
name='mld_aggrmat_biz_bld'
name='mld_aggrmat_smth_bld'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
@ -120,7 +126,6 @@ subroutine mld_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -130,40 +135,45 @@ subroutine mld_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (parms%aggr_filter == mld_filter_mat_)
!
! naggr: number of local aggregates
! nrow: local rows.
!
! Get the diagonal D
adiag = a%get_diag(info)
adiag = a%get_diag(info)
if (info == psb_success_) &
& call psb_realloc(ncol,adiag,info)
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
if (info == psb_success_) call a%csclip(icoo,info,jmax=a%get_nrows())
call icoo%mv_to_lfmt(acsr,info)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies sone.'
& ' Initial copies done.'
call acsr%cp_to_fmt(acsrf,info)
if (filter_mat) then
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1, nrow
tmp = szero
@ -195,39 +205,13 @@ subroutine mld_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
end if
end do
if (filter_mat) call acsrf%scal(adiag,info)
if (info == psb_success_) call acsr3%scal(adiag,info)
if (info /= psb_success_) goto 9999
if (parms%aggr_omega_alg == mld_eig_est_) then
if (parms%aggr_eig == mld_max_norm_) then
!
! This only works with CSR
!
anorm = szero
dg = sone
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
call psb_amx(ictxt,anorm)
omega = 4.d0/(3.d0*anorm)
parms%aggr_omega_val = omega
@ -247,118 +231,51 @@ subroutine mld_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (filter_mat) then
!
! Build the smoothed prolongator using the filtered matrix
!
do i=1,acsrf%get_nrows()
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) then
acsrf%val(j) = sone - omega*acsrf%val(j)
else
acsrf%val(j) = - omega*acsrf%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
!
do i=1,acsr3%get_nrows()
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) == i) then
acsr3%val(j) = sone - omega*acsr3%val(j)
else
acsr3%val(j) = - omega*acsr3%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*A)Ptilde
! Doing it this way means to consider diag(A_i)
!
!
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
end if
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call acsrf%scal(adiag,info)
if (info /= psb_success_) goto 9999
call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call coo_prol%mv_to_fmt(csr_prol,info)
!
! Build the smoothed prolongator using either A or Af
! acsr1 = (I-w*D*A) Prol acsr1 = (I-w*D*Af) Prol
! This is always done through the variable acsrf which
! is a bit less readable, butsaves space and one extra matrix copy
!
call omega_smooth(omega,acsrf)
call psb_par_spspmm(acsrf,desc_a,csr_prol,acsr1,tmp_desc,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call op_prol%mv_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -368,8 +285,26 @@ subroutine mld_saggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
9999 continue
call psb_errpush(info,name)
call psb_error_handler(err_act)
return
contains
subroutine omega_smooth(omega,acsr)
implicit none
real(psb_spk_),intent(in) :: omega
type(psb_ls_csr_sparse_mat), intent(inout) :: acsr
!
integer(psb_lpk_) :: i,j
do i=1,acsr%get_nrows()
do j=acsr%irp(i),acsr%irp(i+1)-1
if (acsr%ja(j) == i) then
acsr%val(j) = sone - omega*acsr%val(j)
else
acsr%val(j) = - omega*acsr%val(j)
end if
end do
end do
end subroutine omega_smooth
end subroutine mld_saggrmat_biz_bld

116
mlprec/impl/aggregator/mld_saggrmat_nosmth_bld.f90
Unescape Escape View File

@ -100,7 +100,7 @@ subroutine mld_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
use psb_base_mod
use mld_base_prec_type
use mld_s_inner_mod, mld_protect_name => mld_saggrmat_nosmth_bld
use mld_s_base_aggregator_mod
implicit none
! Arguments
@ -118,11 +118,13 @@ subroutine mld_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_lsspmat_type) :: la
type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2
type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo, coo_prol, coo_restr
type(psb_ls_csr_sparse_mat) :: acsr1, acsr2, acsr
type(psb_desc_type) :: tmp_desc
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
integer(psb_ipk_) :: inaggr, nzlp
name = 'mld_aggrmat_nosmth_bld'
info = psb_success_
@ -143,53 +145,69 @@ subroutine mld_saggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
if (.false.) then
call a%cp_to(acsr)
call op_prol%mv_to(coo_prol)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
end if
call psb_erractionrestore(err_act)
return

251
mlprec/impl/aggregator/mld_zaggrmat_biz_bld.f90
Unescape Escape View File

@ -84,34 +84,40 @@ subroutine mld_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
use psb_base_mod
use mld_base_prec_type
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_biz_bld
use mld_z_base_aggregator_mod
implicit none
! Arguments
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
type(psb_zspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lzspmat_type), intent(inout) :: op_prol
type(psb_lzspmat_type), intent(out) :: ac,op_restr
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer(psb_ipk_) ::ictxt, np, me
integer(psb_ipk_) :: inaggr, nzlp, inrow
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_lzspmat_type) :: am3, am4,tmp_prol, la
type(psb_lz_coo_sparse_mat) :: tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
type(psb_desc_type) :: tmp_desc
type(psb_z_coo_sparse_mat) :: icoo
type(psb_lz_coo_sparse_mat) :: coo_prol, coo_restr, tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsrf, csr_prol, acsr
complex(psb_dpk_), allocatable :: adiag(:)
real(psb_dpk_), allocatable :: arwsum(:)
integer(psb_ipk_) :: ierr(5)
logical :: filter_mat
integer(psb_ipk_) :: debug_level, debug_unit, err_act
integer(psb_ipk_), parameter :: ncmax=16
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
logical, parameter :: debug_new=.false.
character(len=80) :: filename
name='mld_aggrmat_biz_bld'
name='mld_aggrmat_smth_bld'
info=psb_success_
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
@ -120,7 +126,6 @@ subroutine mld_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -130,40 +135,45 @@ subroutine mld_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
ncol = desc_a%get_local_cols()
theta = parms%aggr_thresh
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (parms%aggr_filter == mld_filter_mat_)
!
! naggr: number of local aggregates
! nrow: local rows.
!
! Get the diagonal D
adiag = a%get_diag(info)
adiag = a%get_diag(info)
if (info == psb_success_) &
& call psb_realloc(ncol,adiag,info)
& call psb_realloc(ncol,adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to_l(la)
if (info /= psb_success_) then
if (info == psb_success_) call a%csclip(icoo,info,jmax=a%get_nrows())
call icoo%mv_to_lfmt(acsr,info)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
goto 9999
end if
! 1. Allocate Ptilde in sparse matrix form
call op_prol%mv_to(tmpcoo)
call ptilde%mv_from_coo(tmpcoo,info)
if (info == psb_success_) call la%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies sone.'
& ' Initial copies done.'
call acsr%cp_to_fmt(acsrf,info)
if (filter_mat) then
!
! Build the filtered matrix Af from A
!
if (info == psb_success_) call acsr3%cp_to_fmt(acsrf,info)
do i=1, nrow
tmp = zzero
@ -195,39 +205,13 @@ subroutine mld_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
end if
end do
if (filter_mat) call acsrf%scal(adiag,info)
if (info == psb_success_) call acsr3%scal(adiag,info)
if (info /= psb_success_) goto 9999
if (parms%aggr_omega_alg == mld_eig_est_) then
if (parms%aggr_eig == mld_max_norm_) then
!
! This only works with CSR
!
anorm = dzero
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = dzero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
endif
if (acsr3%ja(j) == i ) then
dg = abs(acsr3%val(j))
end if
end do
anorm = max(anorm,tmp/dg)
enddo
call psb_amx(ictxt,anorm)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
goto 9999
end if
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
call psb_amx(ictxt,anorm)
omega = 4.d0/(3.d0*anorm)
parms%aggr_omega_val = omega
@ -247,118 +231,51 @@ subroutine mld_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
goto 9999
end if
if (filter_mat) then
!
! Build the smoothed prolongator using the filtered matrix
!
do i=1,acsrf%get_nrows()
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) then
acsrf%val(j) = zone - omega*acsrf%val(j)
else
acsrf%val(j) = - omega*acsrf%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_spspmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
else
!
! Build the smoothed prolongator using the original matrix
!
do i=1,acsr3%get_nrows()
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) == i) then
acsr3%val(j) = zone - omega*acsr3%val(j)
else
acsr3%val(j) = - omega*acsr3%val(j)
end if
end do
end do
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done gather, going for SYMBMM 1'
!
! Symbmm90 does the allocation for its result.
!
! acsrm1 = (I-w*D*A)Ptilde
! Doing it this way means to consider diag(A_i)
!
!
call psb_spspmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
end if
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call acsrf%scal(adiag,info)
if (info /= psb_success_) goto 9999
call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call coo_prol%mv_to_fmt(csr_prol,info)
!
! Build the smoothed prolongator using either A or Af
! acsr1 = (I-w*D*A) Prol acsr1 = (I-w*D*Af) Prol
! This is always done through the variable acsrf which
! is a bit less readable, butsaves space and one extra matrix copy
!
call omega_smooth(omega,acsrf)
call psb_par_spspmm(acsrf,desc_a,csr_prol,acsr1,tmp_desc,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
goto 9999
end if
call psb_spspmm(la,tmp_prol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 2')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 2',parms%aggr_prol, mld_smooth_prol_
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
call tmp_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
call psb_spspmm(op_restr,am3,ac,info)
if (info == psb_success_) call am3%free()
if (info == psb_success_) call ac%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build b = op_restr x am3')
goto 9999
end if
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call op_prol%mv_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -368,8 +285,26 @@ subroutine mld_zaggrmat_biz_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
9999 continue
call psb_errpush(info,name)
call psb_error_handler(err_act)
return
contains
subroutine omega_smooth(omega,acsr)
implicit none
real(psb_dpk_),intent(in) :: omega
type(psb_lz_csr_sparse_mat), intent(inout) :: acsr
!
integer(psb_lpk_) :: i,j
do i=1,acsr%get_nrows()
do j=acsr%irp(i),acsr%irp(i+1)-1
if (acsr%ja(j) == i) then
acsr%val(j) = zone - omega*acsr%val(j)
else
acsr%val(j) = - omega*acsr%val(j)
end if
end do
end do
end subroutine omega_smooth
end subroutine mld_zaggrmat_biz_bld

116
mlprec/impl/aggregator/mld_zaggrmat_nosmth_bld.f90
Unescape Escape View File

@ -100,7 +100,7 @@ subroutine mld_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
use psb_base_mod
use mld_base_prec_type
use mld_z_inner_mod, mld_protect_name => mld_zaggrmat_nosmth_bld
use mld_z_base_aggregator_mod
implicit none
! Arguments
@ -118,11 +118,13 @@ subroutine mld_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_lzspmat_type) :: la
type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2
type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo, coo_prol, coo_restr
type(psb_lz_csr_sparse_mat) :: acsr1, acsr2, acsr
type(psb_desc_type) :: tmp_desc
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& naggr, nzt, naggrm1, naggrp1, i, k
integer(psb_ipk_) :: inaggr, nzlp
name = 'mld_aggrmat_nosmth_bld'
info = psb_success_
@ -143,53 +145,69 @@ subroutine mld_zaggrmat_nosmth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
if (.false.) then
call a%cp_to(acsr)
call op_prol%mv_to(coo_prol)
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = coo_prol%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(coo_prol%ja(1:nzlp),info)
call coo_prol%set_ncols(tmp_desc%get_local_cols())
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call op_prol%mv_from(coo_prol)
call op_restr%mv_from(coo_restr)
else
call a%cp_to(ac_coo)
nzt = ac_coo%get_nzeros()
k = 0
do i = 1, nzt
k = k + 1
ac_coo%ia(k) = ilaggr(ac_coo%ia(i))
ac_coo%ja(k) = ilaggr(ac_coo%ja(i))
ac_coo%val(k) = ac_coo%val(i)
! At this point, there may be negative entries,
! because that's how ILAGGR marks singletons
! If this is the case, roll back K
if ((ac_coo%ia(k)<=0).or.(ac_coo%ja(k)<=0)) k = k-1
enddo
call ac_coo%set_nrows(naggr)
call ac_coo%set_ncols(naggr)
call ac_coo%set_nzeros(k)
call ac_coo%set_dupl(psb_dupl_add_)
call ac_coo%fix(info)
call ac_coo%trim()
call ac%mv_from(ac_coo)
call op_prol%cp_to(tmpcoo)
call op_prol%cscnv(info,type='csr',dupl=psb_dupl_add_)
call tmpcoo%transp()
!
! Now we have to fix this. The only rows of tmpcoo/op_restr that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
nzl = tmpcoo%get_nzeros()
i = 0
do k = 1, nzl
if ((naggrm1 < tmpcoo%ia(k)) .and.(tmpcoo%ia(k) <= naggrp1)) then
i = i+1
tmpcoo%val(i) = tmpcoo%val(k)
tmpcoo%ia(i) = tmpcoo%ia(k)
tmpcoo%ja(i) = tmpcoo%ja(k)
end if
end do
call tmpcoo%set_nzeros(i)
! call tmpcoo%trim()
call op_restr%mv_from(tmpcoo)
call op_restr%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) goto 9999
end if
call psb_erractionrestore(err_act)
return

Write Preview
Loading…
Cancel
Save