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

Merge branch 'unify_aggr_bld' into development

Browse Source
stopcriterion
Salvatore Filippone 5 years ago
parent 2295f16aa7 850317a62b
commit 497984f9e3
92 changed files with 1988 additions and 1392 deletions
Show Stats Download Patch File Download Diff File

34
mlprec/impl/aggregator/mld_c_dec_aggregator_mat_asb.f90
Unescape Escape View File

@ -99,9 +99,9 @@ subroutine mld_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
integer(psb_ipk_), intent(out) :: info
!
integer(psb_ipk_) :: ictxt, np, me
type(psb_lc_coo_sparse_mat) :: acoo, bcoo
type(psb_lc_coo_sparse_mat) :: lacoo, lbcoo
type(psb_c_coo_sparse_mat) :: acoo
type(psb_lc_csr_sparse_mat) :: acsr1
type(psb_lcspmat_type) :: lac, lac1
type(psb_cspmat_type) :: tmp_ac
integer(psb_ipk_) :: i_nr, i_nc, i_nl, nzl
integer(psb_lpk_) :: ntaggr
@ -125,14 +125,14 @@ subroutine mld_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_distr_mat_)
call ac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
call ac%mv_to(lbcoo)
nzl = lbcoo%get_nzeros()
i_nl = nlaggr(me+1)
if (info == psb_success_) call psb_cdall(ictxt,desc_ac,info,nl=i_nl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdins(nzl,lbcoo%ia,lbcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating desc_ac and converting ac')
@ -141,7 +141,7 @@ subroutine mld_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call ac%mv_from(bcoo)
call ac%mv_from(lbcoo)
call ac%set_nrows(desc_ac%get_local_rows())
call ac%set_ncols(desc_ac%get_local_cols())
call ac%set_asb()
@ -164,12 +164,12 @@ subroutine mld_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
call op_prol%set_ncols(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),desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
!call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(lacoo)
nzl = lacoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(lacoo%ia(1:nzl),desc_ac,info,'I')
call lacoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(lacoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -189,9 +189,13 @@ subroutine mld_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_repl_mat_)
!
!
! If we are here, it means we assume that an IPK version of the
! coarse matrix can hold all indices. User beware!
!
call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call tmp_ac%mv_from_l(ac)
if (info == psb_success_) call ac%mv_to(acoo)
if (info == psb_success_) call tmp_ac%mv_from(acoo)
if (info == psb_success_) &
& call psb_gather(ac,tmp_ac,desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999

186
mlprec/impl/aggregator/mld_c_spmm_bld_inner.f90
Unescape Escape View File

@ -35,13 +35,193 @@
! File: mld_daggrmat_nosmth_bld.F90
!
!
subroutine mld_c_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
subroutine mld_c_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
use psb_base_mod
use mld_c_inner_mod
use mld_c_base_aggregator_mod, mld_protect_name => mld_c_spmm_bld_inner
implicit none
! Arguments
type(psb_lc_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lc_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_lcspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
character(len=40) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_lc_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_lc_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr
integer(psb_ipk_) :: debug_level, debug_unit, naggr
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& nzt, naggrm1, naggrp1, i, k
integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave
logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false.
integer(psb_ipk_), save :: idx_spspmm=-1
name='mld_spmm_bld_inner'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nglob = desc_a%get_global_rows()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
if ((do_timings).and.(idx_spspmm==-1)) &
& idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm")
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
!write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr
nrpsave = coo_prol%get_nrows()
ncpsave = coo_prol%get_ncols()
nzpsave = coo_prol%get_nzeros()
!write(0,*)me,' ',name,' input coo_prol ',nrpsave,ncpsave,nzpsave
!
! Here COO_PROL should be with GLOBAL indices on the cols
! and LOCAL indices on the rows.
!
if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',&
& coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),&
& nrow,ntaggr,naggr
call coo_prol%cp_to_fmt(csr_prol,info)
if (debug) write(0,*) me,trim(name),' Product AxPROL ',&
& a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), &
& desc_a%get_local_rows(),desc_a%get_local_cols(),&
& desc_cprol%get_local_rows(),desc_a%get_local_cols()
if (debug) flush(0)
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
if (debug) write(0,*) me,trim(name),' Done AxPROL ',&
& acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),&
& desc_cprol%get_local_rows(),desc_cprol%get_local_cols()
!
! Ok first product done.
!
! Remember that RESTR must be built from PROL after halo extension,
! which is done above in psb_par_spspmm
if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',&
& csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros()
call csr_prol%cp_to_fmt(coo_restr,info)
!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),&
!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros()
call coo_restr%transp()
nzl = coo_restr%get_nzeros()
call desc_cprol%l2gip(coo_restr%ia(1:nzl),info)
i=0
!
! Now we have to fix this. The only rows of the restrictor that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < coo_restr%ia(k)) .and.(coo_restr%ia(k) <= naggrp1)) then
i = i+1
coo_restr%val(i) = coo_restr%val(k)
coo_restr%ia(i) = coo_restr%ia(k)
coo_restr%ja(i) = coo_restr%ja(k)
end if
end do
call coo_restr%set_nzeros(i)
call coo_restr%fix(info)
call coo_restr%cp_to_coo(tmpcoo,info)
!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
nzl = tmpcoo%get_nzeros()
call psb_glob_to_loc(tmpcoo%ia(1:nzl),desc_cprol,info,iact='I',owned=.true.)
call tmpcoo%clean_negidx(info)
nzl = tmpcoo%get_nzeros()
call tmpcoo%set_nrows(desc_cprol%get_local_rows())
call tmpcoo%set_ncols(desc_a%get_local_cols())
!!$ write(0,*)me,' ',name,' after G2L on rows ',tmpcoo%get_nrows(),tmpcoo%get_ncols(),tmpcoo%get_nzeros()
call csr_restr%mv_from_coo(tmpcoo,info)
if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',&
& csr_restr%get_nrows(),csr_restr%get_ncols(), &
& desc_cprol%get_local_rows(),desc_a%get_local_cols(),&
& acsr3%get_nrows(),acsr3%get_ncols()
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
call csr_restr%free()
call ac_csr%mv_to_coo(ac_coo,info)
nza = ac_coo%get_nzeros()
if (debug) write(0,*) me,trim(name),' Fixing ac ',&
& ac_coo%get_nrows(),ac_coo%get_ncols(), nza
call ac_coo%fix(info)
call desc_cprol%indxmap%l2gip(ac_coo%ia(1:nza),info)
call desc_cprol%indxmap%l2gip(ac_coo%ja(1:nza),info)
call ac_coo%set_nrows(ntaggr)
call ac_coo%set_ncols(ntaggr)
if (debug) write(0,*) me,' ',trim(name),' Before mv_from',psb_get_errstatus()
if (info == 0) call ac%mv_from(ac_coo)
if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus()
if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr
! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros()
nza = coo_prol%get_nzeros()
call desc_cprol%indxmap%l2gip(coo_prol%ja(1:nza),info)
if (debug) then
write(0,*) me,' ',trim(name),' Checkpoint at exit'
call psb_barrier(ictxt)
write(0,*) me,' ',trim(name),' Checkpoint through'
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = coo_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_spmm_bld_inner
subroutine mld_c_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
use psb_base_mod
use mld_c_inner_mod
use mld_c_base_aggregator_mod, mld_protect_name => mld_c_old_spmm_bld_inner
implicit none
! Arguments
type(psb_lc_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
@ -218,4 +398,4 @@ subroutine mld_c_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
return
end subroutine mld_c_spmm_bld_inner
end subroutine mld_c_old_spmm_bld_inner

243
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)
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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
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
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

114
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,67 @@ 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))
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 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)
if (info /= psb_success_) goto 9999
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)
end if
call psb_erractionrestore(err_act)
return

34
mlprec/impl/aggregator/mld_caggrmat_smth_bld.f90
Unescape Escape View File

@ -126,7 +126,7 @@ subroutine mld_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_desc_type) :: tmp_desc
type(psb_lc_coo_sparse_mat) :: tmpcoo
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(:)
@ -176,7 +176,7 @@ subroutine mld_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to(acsr)
call op_prol%cp_to(tmpcoo)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -257,10 +257,10 @@ subroutine mld_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = tmpcoo%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(tmpcoo%ja(1:nzlp),info)
call tmpcoo%set_ncols(tmp_desc%get_local_cols())
call tmpcoo%mv_to_fmt(csr_prol,info)
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
@ -278,16 +278,24 @@ subroutine mld_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call acsr1%set_dupl(psb_dupl_add_)
call acsr1%set_ncols(ntaggr)
call op_prol%cp_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
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),&

34
mlprec/impl/aggregator/mld_d_dec_aggregator_mat_asb.f90
Unescape Escape View File

@ -99,9 +99,9 @@ subroutine mld_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
integer(psb_ipk_), intent(out) :: info
!
integer(psb_ipk_) :: ictxt, np, me
type(psb_ld_coo_sparse_mat) :: acoo, bcoo
type(psb_ld_coo_sparse_mat) :: lacoo, lbcoo
type(psb_d_coo_sparse_mat) :: acoo
type(psb_ld_csr_sparse_mat) :: acsr1
type(psb_ldspmat_type) :: lac, lac1
type(psb_dspmat_type) :: tmp_ac
integer(psb_ipk_) :: i_nr, i_nc, i_nl, nzl
integer(psb_lpk_) :: ntaggr
@ -125,14 +125,14 @@ subroutine mld_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_distr_mat_)
call ac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
call ac%mv_to(lbcoo)
nzl = lbcoo%get_nzeros()
i_nl = nlaggr(me+1)
if (info == psb_success_) call psb_cdall(ictxt,desc_ac,info,nl=i_nl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdins(nzl,lbcoo%ia,lbcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating desc_ac and converting ac')
@ -141,7 +141,7 @@ subroutine mld_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call ac%mv_from(bcoo)
call ac%mv_from(lbcoo)
call ac%set_nrows(desc_ac%get_local_rows())
call ac%set_ncols(desc_ac%get_local_cols())
call ac%set_asb()
@ -164,12 +164,12 @@ subroutine mld_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
call op_prol%set_ncols(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),desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
!call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(lacoo)
nzl = lacoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(lacoo%ia(1:nzl),desc_ac,info,'I')
call lacoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(lacoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -189,9 +189,13 @@ subroutine mld_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_repl_mat_)
!
!
! If we are here, it means we assume that an IPK version of the
! coarse matrix can hold all indices. User beware!
!
call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call tmp_ac%mv_from_l(ac)
if (info == psb_success_) call ac%mv_to(acoo)
if (info == psb_success_) call tmp_ac%mv_from(acoo)
if (info == psb_success_) &
& call psb_gather(ac,tmp_ac,desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999

186
mlprec/impl/aggregator/mld_d_spmm_bld_inner.f90
Unescape Escape View File

@ -35,13 +35,193 @@
! File: mld_daggrmat_nosmth_bld.F90
!
!
subroutine mld_d_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
subroutine mld_d_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
use psb_base_mod
use mld_d_inner_mod
use mld_d_base_aggregator_mod, mld_protect_name => mld_d_spmm_bld_inner
implicit none
! Arguments
type(psb_ld_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ld_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_ldspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
character(len=40) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_ld_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_ld_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr
integer(psb_ipk_) :: debug_level, debug_unit, naggr
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& nzt, naggrm1, naggrp1, i, k
integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave
logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false.
integer(psb_ipk_), save :: idx_spspmm=-1
name='mld_spmm_bld_inner'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nglob = desc_a%get_global_rows()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
if ((do_timings).and.(idx_spspmm==-1)) &
& idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm")
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
!write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr
nrpsave = coo_prol%get_nrows()
ncpsave = coo_prol%get_ncols()
nzpsave = coo_prol%get_nzeros()
!write(0,*)me,' ',name,' input coo_prol ',nrpsave,ncpsave,nzpsave
!
! Here COO_PROL should be with GLOBAL indices on the cols
! and LOCAL indices on the rows.
!
if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',&
& coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),&
& nrow,ntaggr,naggr
call coo_prol%cp_to_fmt(csr_prol,info)
if (debug) write(0,*) me,trim(name),' Product AxPROL ',&
& a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), &
& desc_a%get_local_rows(),desc_a%get_local_cols(),&
& desc_cprol%get_local_rows(),desc_a%get_local_cols()
if (debug) flush(0)
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
if (debug) write(0,*) me,trim(name),' Done AxPROL ',&
& acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),&
& desc_cprol%get_local_rows(),desc_cprol%get_local_cols()
!
! Ok first product done.
!
! Remember that RESTR must be built from PROL after halo extension,
! which is done above in psb_par_spspmm
if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',&
& csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros()
call csr_prol%cp_to_fmt(coo_restr,info)
!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),&
!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros()
call coo_restr%transp()
nzl = coo_restr%get_nzeros()
call desc_cprol%l2gip(coo_restr%ia(1:nzl),info)
i=0
!
! Now we have to fix this. The only rows of the restrictor that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < coo_restr%ia(k)) .and.(coo_restr%ia(k) <= naggrp1)) then
i = i+1
coo_restr%val(i) = coo_restr%val(k)
coo_restr%ia(i) = coo_restr%ia(k)
coo_restr%ja(i) = coo_restr%ja(k)
end if
end do
call coo_restr%set_nzeros(i)
call coo_restr%fix(info)
call coo_restr%cp_to_coo(tmpcoo,info)
!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
nzl = tmpcoo%get_nzeros()
call psb_glob_to_loc(tmpcoo%ia(1:nzl),desc_cprol,info,iact='I',owned=.true.)
call tmpcoo%clean_negidx(info)
nzl = tmpcoo%get_nzeros()
call tmpcoo%set_nrows(desc_cprol%get_local_rows())
call tmpcoo%set_ncols(desc_a%get_local_cols())
!!$ write(0,*)me,' ',name,' after G2L on rows ',tmpcoo%get_nrows(),tmpcoo%get_ncols(),tmpcoo%get_nzeros()
call csr_restr%mv_from_coo(tmpcoo,info)
if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',&
& csr_restr%get_nrows(),csr_restr%get_ncols(), &
& desc_cprol%get_local_rows(),desc_a%get_local_cols(),&
& acsr3%get_nrows(),acsr3%get_ncols()
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
call csr_restr%free()
call ac_csr%mv_to_coo(ac_coo,info)
nza = ac_coo%get_nzeros()
if (debug) write(0,*) me,trim(name),' Fixing ac ',&
& ac_coo%get_nrows(),ac_coo%get_ncols(), nza
call ac_coo%fix(info)
call desc_cprol%indxmap%l2gip(ac_coo%ia(1:nza),info)
call desc_cprol%indxmap%l2gip(ac_coo%ja(1:nza),info)
call ac_coo%set_nrows(ntaggr)
call ac_coo%set_ncols(ntaggr)
if (debug) write(0,*) me,' ',trim(name),' Before mv_from',psb_get_errstatus()
if (info == 0) call ac%mv_from(ac_coo)
if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus()
if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr
! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros()
nza = coo_prol%get_nzeros()
call desc_cprol%indxmap%l2gip(coo_prol%ja(1:nza),info)
if (debug) then
write(0,*) me,' ',trim(name),' Checkpoint at exit'
call psb_barrier(ictxt)
write(0,*) me,' ',trim(name),' Checkpoint through'
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = coo_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_spmm_bld_inner
subroutine mld_d_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
use psb_base_mod
use mld_d_inner_mod
use mld_d_base_aggregator_mod, mld_protect_name => mld_d_old_spmm_bld_inner
implicit none
! Arguments
type(psb_ld_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
@ -218,4 +398,4 @@ subroutine mld_d_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
return
end subroutine mld_d_spmm_bld_inner
end subroutine mld_d_old_spmm_bld_inner

243
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)
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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
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
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

114
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,67 @@ 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))
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 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)
if (info /= psb_success_) goto 9999
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)
end if
call psb_erractionrestore(err_act)
return

34
mlprec/impl/aggregator/mld_daggrmat_smth_bld.f90
Unescape Escape View File

@ -126,7 +126,7 @@ subroutine mld_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_desc_type) :: tmp_desc
type(psb_ld_coo_sparse_mat) :: tmpcoo
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(:)
@ -176,7 +176,7 @@ subroutine mld_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to(acsr)
call op_prol%cp_to(tmpcoo)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -257,10 +257,10 @@ subroutine mld_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = tmpcoo%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(tmpcoo%ja(1:nzlp),info)
call tmpcoo%set_ncols(tmp_desc%get_local_cols())
call tmpcoo%mv_to_fmt(csr_prol,info)
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
@ -278,16 +278,24 @@ subroutine mld_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call acsr1%set_dupl(psb_dupl_add_)
call acsr1%set_ncols(ntaggr)
call op_prol%cp_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
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),&

34
mlprec/impl/aggregator/mld_s_dec_aggregator_mat_asb.f90
Unescape Escape View File

@ -99,9 +99,9 @@ subroutine mld_s_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
integer(psb_ipk_), intent(out) :: info
!
integer(psb_ipk_) :: ictxt, np, me
type(psb_ls_coo_sparse_mat) :: acoo, bcoo
type(psb_ls_coo_sparse_mat) :: lacoo, lbcoo
type(psb_s_coo_sparse_mat) :: acoo
type(psb_ls_csr_sparse_mat) :: acsr1
type(psb_lsspmat_type) :: lac, lac1
type(psb_sspmat_type) :: tmp_ac
integer(psb_ipk_) :: i_nr, i_nc, i_nl, nzl
integer(psb_lpk_) :: ntaggr
@ -125,14 +125,14 @@ subroutine mld_s_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_distr_mat_)
call ac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
call ac%mv_to(lbcoo)
nzl = lbcoo%get_nzeros()
i_nl = nlaggr(me+1)
if (info == psb_success_) call psb_cdall(ictxt,desc_ac,info,nl=i_nl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdins(nzl,lbcoo%ia,lbcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating desc_ac and converting ac')
@ -141,7 +141,7 @@ subroutine mld_s_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call ac%mv_from(bcoo)
call ac%mv_from(lbcoo)
call ac%set_nrows(desc_ac%get_local_rows())
call ac%set_ncols(desc_ac%get_local_cols())
call ac%set_asb()
@ -164,12 +164,12 @@ subroutine mld_s_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
call op_prol%set_ncols(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),desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
!call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(lacoo)
nzl = lacoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(lacoo%ia(1:nzl),desc_ac,info,'I')
call lacoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(lacoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -189,9 +189,13 @@ subroutine mld_s_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_repl_mat_)
!
!
! If we are here, it means we assume that an IPK version of the
! coarse matrix can hold all indices. User beware!
!
call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call tmp_ac%mv_from_l(ac)
if (info == psb_success_) call ac%mv_to(acoo)
if (info == psb_success_) call tmp_ac%mv_from(acoo)
if (info == psb_success_) &
& call psb_gather(ac,tmp_ac,desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999

186
mlprec/impl/aggregator/mld_s_spmm_bld_inner.f90
Unescape Escape View File

@ -35,13 +35,193 @@
! File: mld_daggrmat_nosmth_bld.F90
!
!
subroutine mld_s_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
subroutine mld_s_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
use psb_base_mod
use mld_s_inner_mod
use mld_s_base_aggregator_mod, mld_protect_name => mld_s_spmm_bld_inner
implicit none
! Arguments
type(psb_ls_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_ls_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_lsspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
character(len=40) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_ls_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_ls_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr
integer(psb_ipk_) :: debug_level, debug_unit, naggr
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& nzt, naggrm1, naggrp1, i, k
integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave
logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false.
integer(psb_ipk_), save :: idx_spspmm=-1
name='mld_spmm_bld_inner'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nglob = desc_a%get_global_rows()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
if ((do_timings).and.(idx_spspmm==-1)) &
& idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm")
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
!write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr
nrpsave = coo_prol%get_nrows()
ncpsave = coo_prol%get_ncols()
nzpsave = coo_prol%get_nzeros()
!write(0,*)me,' ',name,' input coo_prol ',nrpsave,ncpsave,nzpsave
!
! Here COO_PROL should be with GLOBAL indices on the cols
! and LOCAL indices on the rows.
!
if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',&
& coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),&
& nrow,ntaggr,naggr
call coo_prol%cp_to_fmt(csr_prol,info)
if (debug) write(0,*) me,trim(name),' Product AxPROL ',&
& a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), &
& desc_a%get_local_rows(),desc_a%get_local_cols(),&
& desc_cprol%get_local_rows(),desc_a%get_local_cols()
if (debug) flush(0)
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
if (debug) write(0,*) me,trim(name),' Done AxPROL ',&
& acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),&
& desc_cprol%get_local_rows(),desc_cprol%get_local_cols()
!
! Ok first product done.
!
! Remember that RESTR must be built from PROL after halo extension,
! which is done above in psb_par_spspmm
if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',&
& csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros()
call csr_prol%cp_to_fmt(coo_restr,info)
!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),&
!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros()
call coo_restr%transp()
nzl = coo_restr%get_nzeros()
call desc_cprol%l2gip(coo_restr%ia(1:nzl),info)
i=0
!
! Now we have to fix this. The only rows of the restrictor that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < coo_restr%ia(k)) .and.(coo_restr%ia(k) <= naggrp1)) then
i = i+1
coo_restr%val(i) = coo_restr%val(k)
coo_restr%ia(i) = coo_restr%ia(k)
coo_restr%ja(i) = coo_restr%ja(k)
end if
end do
call coo_restr%set_nzeros(i)
call coo_restr%fix(info)
call coo_restr%cp_to_coo(tmpcoo,info)
!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
nzl = tmpcoo%get_nzeros()
call psb_glob_to_loc(tmpcoo%ia(1:nzl),desc_cprol,info,iact='I',owned=.true.)
call tmpcoo%clean_negidx(info)
nzl = tmpcoo%get_nzeros()
call tmpcoo%set_nrows(desc_cprol%get_local_rows())
call tmpcoo%set_ncols(desc_a%get_local_cols())
!!$ write(0,*)me,' ',name,' after G2L on rows ',tmpcoo%get_nrows(),tmpcoo%get_ncols(),tmpcoo%get_nzeros()
call csr_restr%mv_from_coo(tmpcoo,info)
if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',&
& csr_restr%get_nrows(),csr_restr%get_ncols(), &
& desc_cprol%get_local_rows(),desc_a%get_local_cols(),&
& acsr3%get_nrows(),acsr3%get_ncols()
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
call csr_restr%free()
call ac_csr%mv_to_coo(ac_coo,info)
nza = ac_coo%get_nzeros()
if (debug) write(0,*) me,trim(name),' Fixing ac ',&
& ac_coo%get_nrows(),ac_coo%get_ncols(), nza
call ac_coo%fix(info)
call desc_cprol%indxmap%l2gip(ac_coo%ia(1:nza),info)
call desc_cprol%indxmap%l2gip(ac_coo%ja(1:nza),info)
call ac_coo%set_nrows(ntaggr)
call ac_coo%set_ncols(ntaggr)
if (debug) write(0,*) me,' ',trim(name),' Before mv_from',psb_get_errstatus()
if (info == 0) call ac%mv_from(ac_coo)
if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus()
if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr
! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros()
nza = coo_prol%get_nzeros()
call desc_cprol%indxmap%l2gip(coo_prol%ja(1:nza),info)
if (debug) then
write(0,*) me,' ',trim(name),' Checkpoint at exit'
call psb_barrier(ictxt)
write(0,*) me,' ',trim(name),' Checkpoint through'
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = coo_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_spmm_bld_inner
subroutine mld_s_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
use psb_base_mod
use mld_s_inner_mod
use mld_s_base_aggregator_mod, mld_protect_name => mld_s_old_spmm_bld_inner
implicit none
! Arguments
type(psb_ls_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
@ -218,4 +398,4 @@ subroutine mld_s_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
return
end subroutine mld_s_spmm_bld_inner
end subroutine mld_s_old_spmm_bld_inner

243
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)
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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
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
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

114
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,67 @@ 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))
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 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)
if (info /= psb_success_) goto 9999
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)
end if
call psb_erractionrestore(err_act)
return

34
mlprec/impl/aggregator/mld_saggrmat_smth_bld.f90
Unescape Escape View File

@ -126,7 +126,7 @@ subroutine mld_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_desc_type) :: tmp_desc
type(psb_ls_coo_sparse_mat) :: tmpcoo
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(:)
@ -176,7 +176,7 @@ subroutine mld_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to(acsr)
call op_prol%cp_to(tmpcoo)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -257,10 +257,10 @@ subroutine mld_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = tmpcoo%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(tmpcoo%ja(1:nzlp),info)
call tmpcoo%set_ncols(tmp_desc%get_local_cols())
call tmpcoo%mv_to_fmt(csr_prol,info)
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
@ -278,16 +278,24 @@ subroutine mld_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call acsr1%set_dupl(psb_dupl_add_)
call acsr1%set_ncols(ntaggr)
call op_prol%cp_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
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),&

34
mlprec/impl/aggregator/mld_z_dec_aggregator_mat_asb.f90
Unescape Escape View File

@ -99,9 +99,9 @@ subroutine mld_z_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
integer(psb_ipk_), intent(out) :: info
!
integer(psb_ipk_) :: ictxt, np, me
type(psb_lz_coo_sparse_mat) :: acoo, bcoo
type(psb_lz_coo_sparse_mat) :: lacoo, lbcoo
type(psb_z_coo_sparse_mat) :: acoo
type(psb_lz_csr_sparse_mat) :: acsr1
type(psb_lzspmat_type) :: lac, lac1
type(psb_zspmat_type) :: tmp_ac
integer(psb_ipk_) :: i_nr, i_nc, i_nl, nzl
integer(psb_lpk_) :: ntaggr
@ -125,14 +125,14 @@ subroutine mld_z_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_distr_mat_)
call ac%mv_to(bcoo)
nzl = bcoo%get_nzeros()
call ac%mv_to(lbcoo)
nzl = lbcoo%get_nzeros()
i_nl = nlaggr(me+1)
if (info == psb_success_) call psb_cdall(ictxt,desc_ac,info,nl=i_nl)
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdins(nzl,lbcoo%ia,lbcoo%ja,desc_ac,info)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ia(1:nzl),desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(lbcoo%ja(1:nzl),desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating desc_ac and converting ac')
@ -141,7 +141,7 @@ subroutine mld_z_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call ac%mv_from(bcoo)
call ac%mv_from(lbcoo)
call ac%set_nrows(desc_ac%get_local_rows())
call ac%set_ncols(desc_ac%get_local_cols())
call ac%set_asb()
@ -164,12 +164,12 @@ subroutine mld_z_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
call op_prol%set_ncols(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),desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(acoo)
!call op_restr%cscnv(info,type='coo',dupl=psb_dupl_add_)
call op_restr%mv_to(lacoo)
nzl = lacoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(lacoo%ia(1:nzl),desc_ac,info,'I')
call lacoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call op_restr%mv_from(lacoo)
if (info == psb_success_) call op_restr%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -189,9 +189,13 @@ subroutine mld_z_dec_aggregator_mat_asb(ag,parms,a,desc_a,ilaggr,nlaggr,&
case(mld_repl_mat_)
!
!
! If we are here, it means we assume that an IPK version of the
! coarse matrix can hold all indices. User beware!
!
call psb_cdall(ictxt,desc_ac,info,mg=ntaggr,repl=.true.)
if (info == psb_success_) call psb_cdasb(desc_ac,info)
if (info == psb_success_) call tmp_ac%mv_from_l(ac)
if (info == psb_success_) call ac%mv_to(acoo)
if (info == psb_success_) call tmp_ac%mv_from(acoo)
if (info == psb_success_) &
& call psb_gather(ac,tmp_ac,desc_ac,info,dupl=psb_dupl_add_,keeploc=.false.)
if (info /= psb_success_) goto 9999

186
mlprec/impl/aggregator/mld_z_spmm_bld_inner.f90
Unescape Escape View File

@ -35,13 +35,193 @@
! File: mld_daggrmat_nosmth_bld.F90
!
!
subroutine mld_z_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
subroutine mld_z_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
use psb_base_mod
use mld_z_inner_mod
use mld_z_base_aggregator_mod, mld_protect_name => mld_z_spmm_bld_inner
implicit none
! Arguments
type(psb_lz_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lz_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_lzspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
! Local variables
integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ictxt,np,me, icomm, ndx, minfo
character(len=40) :: name
integer(psb_ipk_) :: ierr(5)
type(psb_lz_coo_sparse_mat) :: ac_coo, tmpcoo
type(psb_lz_csr_sparse_mat) :: acsr3, csr_prol, ac_csr, csr_restr
integer(psb_ipk_) :: debug_level, debug_unit, naggr
integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, nzl, ip, &
& nzt, naggrm1, naggrp1, i, k
integer(psb_lpk_) :: nrsave, ncsave, nzsave, nza, nrpsave, ncpsave, nzpsave
logical, parameter :: do_timings=.true., oldstyle=.false., debug=.false.
integer(psb_ipk_), save :: idx_spspmm=-1
name='mld_spmm_bld_inner'
if(psb_get_errstatus().ne.0) return
info=psb_success_
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nglob = desc_a%get_global_rows()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
if ((do_timings).and.(idx_spspmm==-1)) &
& idx_spspmm = psb_get_timer_idx("SPMM_BLD: par_spspmm")
naggr = nlaggr(me+1)
ntaggr = sum(nlaggr)
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
!write(0,*)me,' ',name,' input sizes',nlaggr(:),':',naggr
nrpsave = coo_prol%get_nrows()
ncpsave = coo_prol%get_ncols()
nzpsave = coo_prol%get_nzeros()
!write(0,*)me,' ',name,' input coo_prol ',nrpsave,ncpsave,nzpsave
!
! Here COO_PROL should be with GLOBAL indices on the cols
! and LOCAL indices on the rows.
!
if (debug) write(0,*) me,' ',trim(name),' Size check on entry New: ',&
& coo_prol%get_fmt(),coo_prol%get_nrows(),coo_prol%get_ncols(),coo_prol%get_nzeros(),&
& nrow,ntaggr,naggr
call coo_prol%cp_to_fmt(csr_prol,info)
if (debug) write(0,*) me,trim(name),' Product AxPROL ',&
& a_csr%get_nrows(),a_csr%get_ncols(), csr_prol%get_nrows(), &
& desc_a%get_local_rows(),desc_a%get_local_cols(),&
& desc_cprol%get_local_rows(),desc_a%get_local_cols()
if (debug) flush(0)
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(a_csr,desc_a,csr_prol,acsr3,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
if (debug) write(0,*) me,trim(name),' Done AxPROL ',&
& acsr3%get_nrows(),acsr3%get_ncols(), acsr3%get_nzeros(),&
& desc_cprol%get_local_rows(),desc_cprol%get_local_cols()
!
! Ok first product done.
!
! Remember that RESTR must be built from PROL after halo extension,
! which is done above in psb_par_spspmm
if (debug) write(0,*)me,' ',name,' No inp_restr, transposing prol ',&
& csr_prol%get_nrows(),csr_prol%get_ncols(),csr_prol%get_nzeros()
call csr_prol%cp_to_fmt(coo_restr,info)
!!$ write(0,*)me,' ',name,' new into transposition ',coo_restr%get_nrows(),&
!!$ & coo_restr%get_ncols(),coo_restr%get_nzeros()
call coo_restr%transp()
nzl = coo_restr%get_nzeros()
call desc_cprol%l2gip(coo_restr%ia(1:nzl),info)
i=0
!
! Now we have to fix this. The only rows of the restrictor that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < coo_restr%ia(k)) .and.(coo_restr%ia(k) <= naggrp1)) then
i = i+1
coo_restr%val(i) = coo_restr%val(k)
coo_restr%ia(i) = coo_restr%ia(k)
coo_restr%ja(i) = coo_restr%ja(k)
end if
end do
call coo_restr%set_nzeros(i)
call coo_restr%fix(info)
call coo_restr%cp_to_coo(tmpcoo,info)
!!$ write(0,*)me,' ',name,' after transposition ',coo_restr%get_nrows(),coo_restr%get_ncols(),coo_restr%get_nzeros()
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv coo_restr')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
nzl = tmpcoo%get_nzeros()
call psb_glob_to_loc(tmpcoo%ia(1:nzl),desc_cprol,info,iact='I',owned=.true.)
call tmpcoo%clean_negidx(info)
nzl = tmpcoo%get_nzeros()
call tmpcoo%set_nrows(desc_cprol%get_local_rows())
call tmpcoo%set_ncols(desc_a%get_local_cols())
!!$ write(0,*)me,' ',name,' after G2L on rows ',tmpcoo%get_nrows(),tmpcoo%get_ncols(),tmpcoo%get_nzeros()
call csr_restr%mv_from_coo(tmpcoo,info)
if (debug) write(0,*) me,trim(name),' Product RESTRxAP ',&
& csr_restr%get_nrows(),csr_restr%get_ncols(), &
& desc_cprol%get_local_rows(),desc_a%get_local_cols(),&
& acsr3%get_nrows(),acsr3%get_ncols()
if (do_timings) call psb_tic(idx_spspmm)
call psb_par_spspmm(csr_restr,desc_a,acsr3,ac_csr,desc_cprol,info)
if (do_timings) call psb_toc(idx_spspmm)
call csr_restr%free()
call ac_csr%mv_to_coo(ac_coo,info)
nza = ac_coo%get_nzeros()
if (debug) write(0,*) me,trim(name),' Fixing ac ',&
& ac_coo%get_nrows(),ac_coo%get_ncols(), nza
call ac_coo%fix(info)
call desc_cprol%indxmap%l2gip(ac_coo%ia(1:nza),info)
call desc_cprol%indxmap%l2gip(ac_coo%ja(1:nza),info)
call ac_coo%set_nrows(ntaggr)
call ac_coo%set_ncols(ntaggr)
if (debug) write(0,*) me,' ',trim(name),' Before mv_from',psb_get_errstatus()
if (info == 0) call ac%mv_from(ac_coo)
if (debug) write(0,*) me,' ',trim(name),' After mv_from',psb_get_errstatus()
if (debug) write(0,*) me,' ',trim(name),' ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros(),naggr,ntaggr
! write(0,*) me,' ',trim(name),' Final AC newstyle ',ac%get_fmt(),ac%get_nrows(),ac%get_ncols(),ac%get_nzeros()
nza = coo_prol%get_nzeros()
call desc_cprol%indxmap%l2gip(coo_prol%ja(1:nza),info)
if (debug) then
write(0,*) me,' ',trim(name),' Checkpoint at exit'
call psb_barrier(ictxt)
write(0,*) me,' ',trim(name),' Checkpoint through'
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Build ac = coo_restr x am3')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done smooth_aggregate '
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_spmm_bld_inner
subroutine mld_z_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
use psb_base_mod
use mld_z_inner_mod
use mld_z_base_aggregator_mod, mld_protect_name => mld_z_old_spmm_bld_inner
implicit none
! Arguments
type(psb_lz_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
@ -218,4 +398,4 @@ subroutine mld_z_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
return
end subroutine mld_z_spmm_bld_inner
end subroutine mld_z_old_spmm_bld_inner

243
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)
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
allocate(arwsum(nrow))
call acsr%arwsum(arwsum)
anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
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
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

114
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,67 @@ 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))
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 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)
if (info /= psb_success_) goto 9999
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)
end if
call psb_erractionrestore(err_act)
return

34
mlprec/impl/aggregator/mld_zaggrmat_smth_bld.f90
Unescape Escape View File

@ -126,7 +126,7 @@ subroutine mld_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
integer(psb_ipk_) :: ictxt, np, me
character(len=20) :: name
type(psb_desc_type) :: tmp_desc
type(psb_lz_coo_sparse_mat) :: tmpcoo
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(:)
@ -176,7 +176,7 @@ subroutine mld_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
if (info == psb_success_) call a%cp_to(acsr)
call op_prol%cp_to(tmpcoo)
call op_prol%mv_to(coo_prol)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
@ -257,10 +257,10 @@ subroutine mld_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
inaggr = naggr
call psb_cdall(ictxt,tmp_desc,info,nl=inaggr)
nzlp = tmpcoo%get_nzeros()
call tmp_desc%indxmap%g2lip_ins(tmpcoo%ja(1:nzlp),info)
call tmpcoo%set_ncols(tmp_desc%get_local_cols())
call tmpcoo%mv_to_fmt(csr_prol,info)
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
@ -278,16 +278,24 @@ subroutine mld_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_rest
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done SPSPMM 1'
if (.true.) then
nzl = acsr1%get_nzeros()
call acsr1%mv_to_coo(coo_prol,info)
nzl = acsr1%get_nzeros()
call tmp_desc%l2gip(acsr1%ja(1:nzl),info)
call acsr1%set_dupl(psb_dupl_add_)
call acsr1%set_ncols(ntaggr)
call op_prol%cp_from(acsr1)
call mld_spmm_bld_inner(acsr,desc_a,nlaggr,parms,ac,&
& coo_prol,tmp_desc,coo_restr,info)
call mld_spmm_bld_inner(acsr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
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),&

72
mlprec/impl/level/mld_c_base_onelev_build.f90
Unescape Escape View File

@ -45,10 +45,8 @@ subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold)
class(psb_c_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize
integer(psb_ipk_) :: err,i,k, err_act
integer(psb_ipk_) :: ictxt, me, np
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -68,7 +66,6 @@ subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold)
goto 9999
end if
info = psb_success_
int_err(1) = 0
ictxt = lv%base_desc%get_ctxt()
call psb_info(ictxt,me,np)
@ -86,8 +83,19 @@ subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold)
end if
lv%ac_nz_loc = lv%ac%get_nzeros()
lv%ac_nz_tot = lv%ac_nz_loc
call psb_sum(ictxt,lv%ac_nz_tot)
! end do
select case(lv%parms%coarse_mat)
case(mld_distr_mat_)
call psb_sum(ictxt,lv%ac_nz_tot)
case(mld_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
@ -96,44 +104,24 @@ subroutine mld_c_base_onelev_build(lv,info,amold,vmold,imold)
call mld_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative)
if (.false.) then
call lv%sm%build(lv%base_a,lv%base_desc,&
& info,amold=amold,vmold=vmold,imold=imold)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info,&
& amold=amold,vmold=vmold,imold=imold)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
else
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
call psb_erractionrestore(err_act)
return

18
mlprec/impl/level/mld_c_base_onelev_dump.f90
Unescape Escape View File

@ -47,12 +47,12 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num
! Local variables
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: i, j, il1, iln, lname, lev, ni
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
character(len=80) :: prefix_, frmt
character(len=1024) :: fname
logical :: ac_, rp_, tprol_, global_num_
integer(psb_ipk_), allocatable :: ivr(:), ivc(:)
integer(psb_lpk_), allocatable :: ivr(:), ivc(:)
info = 0
@ -68,6 +68,7 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
else
icontxt = -1
iam = -1
np = -1
end if
if (present(ac)) then
ac_ = ac
@ -91,8 +92,13 @@ subroutine mld_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (np > 0) then
ni = floor(log10(1.0*np)) + 1
write(frmt,'(a,i3.3,a,i3.3,a)') '(a,i',ni,'.',ni,')'
write(fname(lname+1:lname+ni+2),frmt) '_p',iam
lname = lname + ni + 2
end if
if (global_num_) then
if (level >= 2) then

3
mlprec/impl/level/mld_c_base_onelev_mat_asb.f90
Unescape Escape View File

@ -105,11 +105,10 @@ subroutine mld_c_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
character(len=24) :: name
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
type(psb_lcspmat_type) :: lac, lac1, op_restr
type(psb_lcspmat_type) :: lac, op_restr
type(psb_cspmat_type) :: ac, iop_restr, iop_prol
type(psb_lc_coo_sparse_mat) :: acoo, bcoo
type(psb_lc_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit

72
mlprec/impl/level/mld_d_base_onelev_build.f90
Unescape Escape View File

@ -45,10 +45,8 @@ subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold)
class(psb_d_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize
integer(psb_ipk_) :: err,i,k, err_act
integer(psb_ipk_) :: ictxt, me, np
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -68,7 +66,6 @@ subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold)
goto 9999
end if
info = psb_success_
int_err(1) = 0
ictxt = lv%base_desc%get_ctxt()
call psb_info(ictxt,me,np)
@ -86,8 +83,19 @@ subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold)
end if
lv%ac_nz_loc = lv%ac%get_nzeros()
lv%ac_nz_tot = lv%ac_nz_loc
call psb_sum(ictxt,lv%ac_nz_tot)
! end do
select case(lv%parms%coarse_mat)
case(mld_distr_mat_)
call psb_sum(ictxt,lv%ac_nz_tot)
case(mld_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
@ -96,44 +104,24 @@ subroutine mld_d_base_onelev_build(lv,info,amold,vmold,imold)
call mld_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative)
if (.false.) then
call lv%sm%build(lv%base_a,lv%base_desc,&
& info,amold=amold,vmold=vmold,imold=imold)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info,&
& amold=amold,vmold=vmold,imold=imold)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
else
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
call psb_erractionrestore(err_act)
return

18
mlprec/impl/level/mld_d_base_onelev_dump.f90
Unescape Escape View File

@ -47,12 +47,12 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num
! Local variables
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: i, j, il1, iln, lname, lev, ni
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
character(len=80) :: prefix_, frmt
character(len=1024) :: fname
logical :: ac_, rp_, tprol_, global_num_
integer(psb_ipk_), allocatable :: ivr(:), ivc(:)
integer(psb_lpk_), allocatable :: ivr(:), ivc(:)
info = 0
@ -68,6 +68,7 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
else
icontxt = -1
iam = -1
np = -1
end if
if (present(ac)) then
ac_ = ac
@ -91,8 +92,13 @@ subroutine mld_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (np > 0) then
ni = floor(log10(1.0*np)) + 1
write(frmt,'(a,i3.3,a,i3.3,a)') '(a,i',ni,'.',ni,')'
write(fname(lname+1:lname+ni+2),frmt) '_p',iam
lname = lname + ni + 2
end if
if (global_num_) then
if (level >= 2) then

3
mlprec/impl/level/mld_d_base_onelev_mat_asb.f90
Unescape Escape View File

@ -105,11 +105,10 @@ subroutine mld_d_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
character(len=24) :: name
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
type(psb_ldspmat_type) :: lac, lac1, op_restr
type(psb_ldspmat_type) :: lac, op_restr
type(psb_dspmat_type) :: ac, iop_restr, iop_prol
type(psb_ld_coo_sparse_mat) :: acoo, bcoo
type(psb_ld_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit

72
mlprec/impl/level/mld_s_base_onelev_build.f90
Unescape Escape View File

@ -45,10 +45,8 @@ subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold)
class(psb_s_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize
integer(psb_ipk_) :: err,i,k, err_act
integer(psb_ipk_) :: ictxt, me, np
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -68,7 +66,6 @@ subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold)
goto 9999
end if
info = psb_success_
int_err(1) = 0
ictxt = lv%base_desc%get_ctxt()
call psb_info(ictxt,me,np)
@ -86,8 +83,19 @@ subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold)
end if
lv%ac_nz_loc = lv%ac%get_nzeros()
lv%ac_nz_tot = lv%ac_nz_loc
call psb_sum(ictxt,lv%ac_nz_tot)
! end do
select case(lv%parms%coarse_mat)
case(mld_distr_mat_)
call psb_sum(ictxt,lv%ac_nz_tot)
case(mld_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
@ -96,44 +104,24 @@ subroutine mld_s_base_onelev_build(lv,info,amold,vmold,imold)
call mld_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative)
if (.false.) then
call lv%sm%build(lv%base_a,lv%base_desc,&
& info,amold=amold,vmold=vmold,imold=imold)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info,&
& amold=amold,vmold=vmold,imold=imold)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
else
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
call psb_erractionrestore(err_act)
return

18
mlprec/impl/level/mld_s_base_onelev_dump.f90
Unescape Escape View File

@ -47,12 +47,12 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num
! Local variables
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: i, j, il1, iln, lname, lev, ni
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
character(len=80) :: prefix_, frmt
character(len=1024) :: fname
logical :: ac_, rp_, tprol_, global_num_
integer(psb_ipk_), allocatable :: ivr(:), ivc(:)
integer(psb_lpk_), allocatable :: ivr(:), ivc(:)
info = 0
@ -68,6 +68,7 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
else
icontxt = -1
iam = -1
np = -1
end if
if (present(ac)) then
ac_ = ac
@ -91,8 +92,13 @@ subroutine mld_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (np > 0) then
ni = floor(log10(1.0*np)) + 1
write(frmt,'(a,i3.3,a,i3.3,a)') '(a,i',ni,'.',ni,')'
write(fname(lname+1:lname+ni+2),frmt) '_p',iam
lname = lname + ni + 2
end if
if (global_num_) then
if (level >= 2) then

3
mlprec/impl/level/mld_s_base_onelev_mat_asb.f90
Unescape Escape View File

@ -105,11 +105,10 @@ subroutine mld_s_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
character(len=24) :: name
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
type(psb_lsspmat_type) :: lac, lac1, op_restr
type(psb_lsspmat_type) :: lac, op_restr
type(psb_sspmat_type) :: ac, iop_restr, iop_prol
type(psb_ls_coo_sparse_mat) :: acoo, bcoo
type(psb_ls_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit

72
mlprec/impl/level/mld_z_base_onelev_build.f90
Unescape Escape View File

@ -45,10 +45,8 @@ subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold)
class(psb_z_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize
integer(psb_ipk_) :: err,i,k, err_act
integer(psb_ipk_) :: ictxt, me, np
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -68,7 +66,6 @@ subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold)
goto 9999
end if
info = psb_success_
int_err(1) = 0
ictxt = lv%base_desc%get_ctxt()
call psb_info(ictxt,me,np)
@ -86,8 +83,19 @@ subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold)
end if
lv%ac_nz_loc = lv%ac%get_nzeros()
lv%ac_nz_tot = lv%ac_nz_loc
call psb_sum(ictxt,lv%ac_nz_tot)
! end do
select case(lv%parms%coarse_mat)
case(mld_distr_mat_)
call psb_sum(ictxt,lv%ac_nz_tot)
case(mld_repl_mat_)
! Do nothing
case default
! Should never get here
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='Wrong lv%parms')
goto 9999
end select
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Calling mlprcbld at level ',i
@ -96,44 +104,24 @@ subroutine mld_z_base_onelev_build(lv,info,amold,vmold,imold)
call mld_check_def(lv%parms%sweeps_post,&
& 'Jacobi sweeps',izero,is_int_non_negative)
if (.false.) then
call lv%sm%build(lv%base_a,lv%base_desc,&
& info,amold=amold,vmold=vmold,imold=imold)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info,&
& amold=amold,vmold=vmold,imold=imold)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
else
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
call lv%sm%build(lv%base_a,lv%base_desc,info)
if (info == 0) then
if (allocated(lv%sm2a)) then
call lv%sm2a%build(lv%base_a,lv%base_desc,info)
lv%sm2 => lv%sm2a
else
lv%sm2 => lv%sm
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
end if
if (info /=0 ) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Smoother bld error')
goto 9999
end if
if (any((/present(amold),present(vmold),present(imold)/))) &
& call lv%cnv(info,amold=amold,vmold=vmold,imold=imold)
call psb_erractionrestore(err_act)
return

18
mlprec/impl/level/mld_z_base_onelev_dump.f90
Unescape Escape View File

@ -47,12 +47,12 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: ac, rp, smoother, solver, tprol, global_num
! Local variables
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: i, j, il1, iln, lname, lev, ni
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
character(len=80) :: prefix_, frmt
character(len=1024) :: fname
logical :: ac_, rp_, tprol_, global_num_
integer(psb_ipk_), allocatable :: ivr(:), ivc(:)
integer(psb_lpk_), allocatable :: ivr(:), ivc(:)
info = 0
@ -68,6 +68,7 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
else
icontxt = -1
iam = -1
np = -1
end if
if (present(ac)) then
ac_ = ac
@ -91,8 +92,13 @@ subroutine mld_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
if (np > 0) then
ni = floor(log10(1.0*np)) + 1
write(frmt,'(a,i3.3,a,i3.3,a)') '(a,i',ni,'.',ni,')'
write(fname(lname+1:lname+ni+2),frmt) '_p',iam
lname = lname + ni + 2
end if
if (global_num_) then
if (level >= 2) then

3
mlprec/impl/level/mld_z_base_onelev_mat_asb.f90
Unescape Escape View File

@ -105,11 +105,10 @@ subroutine mld_z_base_onelev_mat_asb(lv,a,desc_a,ilaggr,nlaggr,op_prol,info)
character(len=24) :: name
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: err_act
type(psb_lzspmat_type) :: lac, lac1, op_restr
type(psb_lzspmat_type) :: lac, op_restr
type(psb_zspmat_type) :: ac, iop_restr, iop_prol
type(psb_lz_coo_sparse_mat) :: acoo, bcoo
type(psb_lz_csr_sparse_mat) :: acsr1
integer(psb_lpk_) :: ntaggr, nr, nc
integer(psb_ipk_) :: nzl, inl
integer(psb_ipk_) :: debug_level, debug_unit

8
mlprec/impl/mld_c_hierarchy_bld.f90
Unescape Escape View File

@ -90,7 +90,6 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info)
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_lcspmat_type) :: op_prol
type(mld_c_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -105,7 +104,6 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info)
name = 'mld_c_hierarchy_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
prec%ictxt = ictxt
@ -193,9 +191,9 @@ subroutine mld_c_hierarchy_bld(a,desc_a,prec,info)
! Default to the cubic root of the size at base level.
!
casize = desc_a%get_global_rows()
casize = int((sone*casize)**(sone/(sone*3)),psb_ipk_)
casize = max(casize,ione)
casize = casize*40_psb_ipk_
casize = int((sone*casize)**(sone/(sone*3)),psb_lpk_)
casize = max(casize,lone)
casize = casize*40_psb_lpk_
call psb_bcast(ictxt,casize)
if (casize > huge(prec%ag_data%min_coarse_size)) then
!

6
mlprec/impl/mld_c_smoothers_bld.f90
Unescape Escape View File

@ -97,10 +97,9 @@ subroutine mld_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! Local Variables
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, nplevs, mxplevs
real(psb_spk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val, coarse_solve_id
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: coarse_solve_id
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -115,7 +114,6 @@ subroutine mld_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
name = 'mld_c_smoothers_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

4
mlprec/impl/mld_cmlprec_aply.f90
Unescape Escape View File

@ -224,7 +224,7 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
! Local variables
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act
integer(psb_ipk_) :: nlev, nc2l, level, isweep, err_act
character(len=20) :: name
character :: trans_
complex(psb_spk_) :: beta_
@ -638,7 +638,7 @@ contains
end if
endif
!
! Compute the residual and call recursively
! Compute the residual for next level and call recursively
!
if (pre) then
call psb_geaxpby(cone,vx2l,&

2
mlprec/impl/mld_cmlprec_bld.f90
Unescape Escape View File

@ -98,7 +98,6 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
real(psb_spk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -113,7 +112,6 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
name = 'mld_cmlprec_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

8
mlprec/impl/mld_d_hierarchy_bld.f90
Unescape Escape View File

@ -90,7 +90,6 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info)
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_ldspmat_type) :: op_prol
type(mld_d_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -105,7 +104,6 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info)
name = 'mld_d_hierarchy_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
prec%ictxt = ictxt
@ -193,9 +191,9 @@ subroutine mld_d_hierarchy_bld(a,desc_a,prec,info)
! Default to the cubic root of the size at base level.
!
casize = desc_a%get_global_rows()
casize = int((done*casize)**(done/(done*3)),psb_ipk_)
casize = max(casize,ione)
casize = casize*40_psb_ipk_
casize = int((done*casize)**(done/(done*3)),psb_lpk_)
casize = max(casize,lone)
casize = casize*40_psb_lpk_
call psb_bcast(ictxt,casize)
if (casize > huge(prec%ag_data%min_coarse_size)) then
!

6
mlprec/impl/mld_d_smoothers_bld.f90
Unescape Escape View File

@ -97,10 +97,9 @@ subroutine mld_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! Local Variables
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, nplevs, mxplevs
real(psb_dpk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val, coarse_solve_id
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: coarse_solve_id
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -115,7 +114,6 @@ subroutine mld_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
name = 'mld_d_smoothers_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

4
mlprec/impl/mld_dmlprec_aply.f90
Unescape Escape View File

@ -224,7 +224,7 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
! Local variables
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act
integer(psb_ipk_) :: nlev, nc2l, level, isweep, err_act
character(len=20) :: name
character :: trans_
real(psb_dpk_) :: beta_
@ -638,7 +638,7 @@ contains
end if
endif
!
! Compute the residual and call recursively
! Compute the residual for next level and call recursively
!
if (pre) then
call psb_geaxpby(done,vx2l,&

2
mlprec/impl/mld_dmlprec_bld.f90
Unescape Escape View File

@ -98,7 +98,6 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
real(psb_dpk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -113,7 +112,6 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
name = 'mld_dmlprec_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

8
mlprec/impl/mld_s_hierarchy_bld.f90
Unescape Escape View File

@ -90,7 +90,6 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info)
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_lsspmat_type) :: op_prol
type(mld_s_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -105,7 +104,6 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info)
name = 'mld_s_hierarchy_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
prec%ictxt = ictxt
@ -193,9 +191,9 @@ subroutine mld_s_hierarchy_bld(a,desc_a,prec,info)
! Default to the cubic root of the size at base level.
!
casize = desc_a%get_global_rows()
casize = int((sone*casize)**(sone/(sone*3)),psb_ipk_)
casize = max(casize,ione)
casize = casize*40_psb_ipk_
casize = int((sone*casize)**(sone/(sone*3)),psb_lpk_)
casize = max(casize,lone)
casize = casize*40_psb_lpk_
call psb_bcast(ictxt,casize)
if (casize > huge(prec%ag_data%min_coarse_size)) then
!

6
mlprec/impl/mld_s_smoothers_bld.f90
Unescape Escape View File

@ -97,10 +97,9 @@ subroutine mld_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! Local Variables
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, nplevs, mxplevs
real(psb_spk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val, coarse_solve_id
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: coarse_solve_id
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -115,7 +114,6 @@ subroutine mld_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
name = 'mld_s_smoothers_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

4
mlprec/impl/mld_smlprec_aply.f90
Unescape Escape View File

@ -224,7 +224,7 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
! Local variables
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act
integer(psb_ipk_) :: nlev, nc2l, level, isweep, err_act
character(len=20) :: name
character :: trans_
real(psb_spk_) :: beta_
@ -638,7 +638,7 @@ contains
end if
endif
!
! Compute the residual and call recursively
! Compute the residual for next level and call recursively
!
if (pre) then
call psb_geaxpby(sone,vx2l,&

2
mlprec/impl/mld_smlprec_bld.f90
Unescape Escape View File

@ -98,7 +98,6 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
real(psb_spk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -113,7 +112,6 @@ subroutine mld_smlprec_bld(a,desc_a,p,info,amold,vmold,imold)
name = 'mld_smlprec_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

8
mlprec/impl/mld_z_hierarchy_bld.f90
Unescape Escape View File

@ -90,7 +90,6 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info)
integer(psb_lpk_), allocatable :: ilaggr(:), nlaggr(:)
type(psb_lzspmat_type) :: op_prol
type(mld_z_onelev_type), allocatable :: tprecv(:)
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -105,7 +104,6 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info)
name = 'mld_z_hierarchy_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
prec%ictxt = ictxt
@ -193,9 +191,9 @@ subroutine mld_z_hierarchy_bld(a,desc_a,prec,info)
! Default to the cubic root of the size at base level.
!
casize = desc_a%get_global_rows()
casize = int((done*casize)**(done/(done*3)),psb_ipk_)
casize = max(casize,ione)
casize = casize*40_psb_ipk_
casize = int((done*casize)**(done/(done*3)),psb_lpk_)
casize = max(casize,lone)
casize = casize*40_psb_lpk_
call psb_bcast(ictxt,casize)
if (casize > huge(prec%ag_data%min_coarse_size)) then
!

6
mlprec/impl/mld_z_smoothers_bld.f90
Unescape Escape View File

@ -97,10 +97,9 @@ subroutine mld_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
! Local Variables
integer(psb_ipk_) :: ictxt, me,np
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, nplevs, mxplevs
real(psb_dpk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val, coarse_solve_id
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: coarse_solve_id
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -115,7 +114,6 @@ subroutine mld_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
name = 'mld_z_smoothers_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

4
mlprec/impl/mld_zmlprec_aply.f90
Unescape Escape View File

@ -224,7 +224,7 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
! Local variables
integer(psb_ipk_) :: ictxt, np, me
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act
integer(psb_ipk_) :: nlev, nc2l, level, isweep, err_act
character(len=20) :: name
character :: trans_
complex(psb_dpk_) :: beta_
@ -638,7 +638,7 @@ contains
end if
endif
!
! Compute the residual and call recursively
! Compute the residual for next level and call recursively
!
if (pre) then
call psb_geaxpby(zone,vx2l,&

2
mlprec/impl/mld_zmlprec_bld.f90
Unescape Escape View File

@ -98,7 +98,6 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
integer(psb_ipk_) :: err,i,k, err_act, iszv, newsz, casize, nplevs, mxplevs
real(psb_dpk_) :: mnaggratio
integer(psb_ipk_) :: ipv(mld_ifpsz_), val
integer(psb_ipk_) :: int_err(5)
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name, ch_err
@ -113,7 +112,6 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info,amold,vmold,imold)
name = 'mld_zmlprec_bld'
info = psb_success_
int_err(1) = 0
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)

8
mlprec/impl/smoother/mld_c_as_smoother_bld.f90
Unescape Escape View File

@ -52,7 +52,7 @@ subroutine mld_c_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
! Local variables
type(psb_cspmat_type) :: blck, atmp
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_
complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='c_as_smoother_bld', ch_err
@ -168,10 +168,8 @@ subroutine mld_c_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
nzeros = sm%nd%get_nzeros()
!!$ write(0,*) me,' ND nzeors ',nzeros
call psb_sum(ictxt,nzeros)
sm%nd_nnz_tot = nzeros
sm%nd_nnz_tot = sm%nd%get_nzeros()
call psb_sum(ictxt,sm%nd_nnz_tot)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_c_jac_smoother_apply.f90
Unescape Escape View File

@ -103,7 +103,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)

2
mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -105,7 +105,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)

8
mlprec/impl/smoother/mld_c_jac_smoother_bld.f90
Unescape Escape View File

@ -75,13 +75,13 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
class is (mld_c_diag_solver_type)
call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
sm%nd_nnz_tot = nztota
class default
if (smsv%is_global()) then
! Do not put anything into SM%ND since the solver
! is acting globally.
sm%nnz_nd_tot = 0
sm%nd_nnz_tot = 0
else
call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
@ -94,7 +94,7 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_)
endif
end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
sm%nd_nnz_tot = sm%nd%get_nzeros()
end if
end select
if (info /= psb_success_) then
@ -102,7 +102,7 @@ subroutine mld_c_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='clip & psb_spcnv csr 4')
goto 9999
end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call psb_sum(ictxt,sm%nd_nnz_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)

2
mlprec/impl/smoother/mld_c_jac_smoother_clone.f90
Unescape Escape View File

@ -66,7 +66,7 @@ subroutine mld_c_jac_smoother_clone(sm,smout,info)
select type(smo => smout)
type is (mld_c_jac_smoother_type)
smo%nnz_nd_tot = sm%nnz_nd_tot
smo%nd_nnz_tot = sm%nd_nnz_tot
call sm%nd%clone(smo%nd,info)
if ((info==psb_success_).and.(allocated(sm%sv))) then
allocate(smout%sv,mold=sm%sv,stat=info)

8
mlprec/impl/smoother/mld_d_as_smoother_bld.f90
Unescape Escape View File

@ -52,7 +52,7 @@ subroutine mld_d_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
! Local variables
type(psb_dspmat_type) :: blck, atmp
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='d_as_smoother_bld', ch_err
@ -168,10 +168,8 @@ subroutine mld_d_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
nzeros = sm%nd%get_nzeros()
!!$ write(0,*) me,' ND nzeors ',nzeros
call psb_sum(ictxt,nzeros)
sm%nd_nnz_tot = nzeros
sm%nd_nnz_tot = sm%nd%get_nzeros()
call psb_sum(ictxt,sm%nd_nnz_tot)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_d_jac_smoother_apply.f90
Unescape Escape View File

@ -103,7 +103,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)

2
mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -105,7 +105,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)

8
mlprec/impl/smoother/mld_d_jac_smoother_bld.f90
Unescape Escape View File

@ -75,13 +75,13 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
class is (mld_d_diag_solver_type)
call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
sm%nd_nnz_tot = nztota
class default
if (smsv%is_global()) then
! Do not put anything into SM%ND since the solver
! is acting globally.
sm%nnz_nd_tot = 0
sm%nd_nnz_tot = 0
else
call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
@ -94,7 +94,7 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_)
endif
end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
sm%nd_nnz_tot = sm%nd%get_nzeros()
end if
end select
if (info /= psb_success_) then
@ -102,7 +102,7 @@ subroutine mld_d_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='clip & psb_spcnv csr 4')
goto 9999
end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call psb_sum(ictxt,sm%nd_nnz_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)

2
mlprec/impl/smoother/mld_d_jac_smoother_clone.f90
Unescape Escape View File

@ -66,7 +66,7 @@ subroutine mld_d_jac_smoother_clone(sm,smout,info)
select type(smo => smout)
type is (mld_d_jac_smoother_type)
smo%nnz_nd_tot = sm%nnz_nd_tot
smo%nd_nnz_tot = sm%nd_nnz_tot
call sm%nd%clone(smo%nd,info)
if ((info==psb_success_).and.(allocated(sm%sv))) then
allocate(smout%sv,mold=sm%sv,stat=info)

8
mlprec/impl/smoother/mld_s_as_smoother_bld.f90
Unescape Escape View File

@ -52,7 +52,7 @@ subroutine mld_s_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
! Local variables
type(psb_sspmat_type) :: blck, atmp
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_
real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='s_as_smoother_bld', ch_err
@ -168,10 +168,8 @@ subroutine mld_s_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
nzeros = sm%nd%get_nzeros()
!!$ write(0,*) me,' ND nzeors ',nzeros
call psb_sum(ictxt,nzeros)
sm%nd_nnz_tot = nzeros
sm%nd_nnz_tot = sm%nd%get_nzeros()
call psb_sum(ictxt,sm%nd_nnz_tot)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_s_jac_smoother_apply.f90
Unescape Escape View File

@ -103,7 +103,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)

2
mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -105,7 +105,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)

8
mlprec/impl/smoother/mld_s_jac_smoother_bld.f90
Unescape Escape View File

@ -75,13 +75,13 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
class is (mld_s_diag_solver_type)
call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
sm%nd_nnz_tot = nztota
class default
if (smsv%is_global()) then
! Do not put anything into SM%ND since the solver
! is acting globally.
sm%nnz_nd_tot = 0
sm%nd_nnz_tot = 0
else
call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
@ -94,7 +94,7 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_)
endif
end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
sm%nd_nnz_tot = sm%nd%get_nzeros()
end if
end select
if (info /= psb_success_) then
@ -102,7 +102,7 @@ subroutine mld_s_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='clip & psb_spcnv csr 4')
goto 9999
end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call psb_sum(ictxt,sm%nd_nnz_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)

2
mlprec/impl/smoother/mld_s_jac_smoother_clone.f90
Unescape Escape View File

@ -66,7 +66,7 @@ subroutine mld_s_jac_smoother_clone(sm,smout,info)
select type(smo => smout)
type is (mld_s_jac_smoother_type)
smo%nnz_nd_tot = sm%nnz_nd_tot
smo%nd_nnz_tot = sm%nd_nnz_tot
call sm%nd%clone(smo%nd,info)
if ((info==psb_success_).and.(allocated(sm%sv))) then
allocate(smout%sv,mold=sm%sv,stat=info)

8
mlprec/impl/smoother/mld_z_as_smoother_bld.f90
Unescape Escape View File

@ -52,7 +52,7 @@ subroutine mld_z_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
! Local variables
type(psb_zspmat_type) :: blck, atmp
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros
integer(psb_ipk_) :: n_row,n_col, nrow_a, nhalo, novr, data_
complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='z_as_smoother_bld', ch_err
@ -168,10 +168,8 @@ subroutine mld_z_as_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999
end if
nzeros = sm%nd%get_nzeros()
!!$ write(0,*) me,' ND nzeors ',nzeros
call psb_sum(ictxt,nzeros)
sm%nd_nnz_tot = nzeros
sm%nd_nnz_tot = sm%nd%get_nzeros()
call psb_sum(ictxt,sm%nd_nnz_tot)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'

2
mlprec/impl/smoother/mld_z_jac_smoother_apply.f90
Unescape Escape View File

@ -103,7 +103,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)

2
mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -105,7 +105,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
end if
endif
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then
if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nd_nnz_tot==0))) then
! if .not.sv%is_iterative, there's no need to pass init
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,wv,info)

8
mlprec/impl/smoother/mld_z_jac_smoother_bld.f90
Unescape Escape View File

@ -75,13 +75,13 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
class is (mld_z_diag_solver_type)
call sm%nd%free()
sm%pa => a
sm%nnz_nd_tot = nztota
sm%nd_nnz_tot = nztota
class default
if (smsv%is_global()) then
! Do not put anything into SM%ND since the solver
! is acting globally.
sm%nnz_nd_tot = 0
sm%nd_nnz_tot = 0
else
call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
@ -94,7 +94,7 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& type='csr',dupl=psb_dupl_add_)
endif
end if
sm%nnz_nd_tot = sm%nd%get_nzeros()
sm%nd_nnz_tot = sm%nd%get_nzeros()
end if
end select
if (info /= psb_success_) then
@ -102,7 +102,7 @@ subroutine mld_z_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
& a_err='clip & psb_spcnv csr 4')
goto 9999
end if
call psb_sum(ictxt,sm%nnz_nd_tot)
call psb_sum(ictxt,sm%nd_nnz_tot)
call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold)

2
mlprec/impl/smoother/mld_z_jac_smoother_clone.f90
Unescape Escape View File

@ -66,7 +66,7 @@ subroutine mld_z_jac_smoother_clone(sm,smout,info)
select type(smo => smout)
type is (mld_z_jac_smoother_type)
smo%nnz_nd_tot = sm%nnz_nd_tot
smo%nd_nnz_tot = sm%nd_nnz_tot
call sm%nd%clone(smo%nd,info)
if ((info==psb_success_).and.(allocated(sm%sv))) then
allocate(smout%sv,mold=sm%sv,stat=info)

5
mlprec/mld_c_as_smoother.f90
Unescape Escape View File

@ -67,7 +67,8 @@ module mld_c_as_smoother
!
type(psb_cspmat_type) :: nd
type(psb_desc_type) :: desc_data
integer(psb_ipk_) :: novr, restr, prol, nd_nnz_tot
integer(psb_ipk_) :: novr, restr, prol
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: check => mld_c_as_smoother_check
procedure, pass(sm) :: dump => mld_c_as_smoother_dmp
@ -321,7 +322,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = 3*psb_sizeof_ip + psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

20
mlprec/mld_c_base_aggregator_mod.f90
Unescape Escape View File

@ -44,7 +44,7 @@ module mld_c_base_aggregator_mod
use mld_base_prec_type, only : mld_sml_parms, mld_saggr_data
use psb_base_mod, only : psb_cspmat_type, psb_lcspmat_type, psb_c_vect_type, &
& psb_c_base_vect_type, psb_clinmap_type, psb_spk_, psb_lc_csr_sparse_mat, &
& psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_lc_coo_sparse_mat, psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_erractionsave, psb_error_handler, psb_success_, psb_toupper
!
!
@ -119,7 +119,21 @@ module mld_c_base_aggregator_mod
end interface
interface mld_spmm_bld_inner
subroutine mld_c_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
subroutine mld_c_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
import :: psb_lc_csr_sparse_mat, psb_lcspmat_type, psb_desc_type, &
& psb_lc_coo_sparse_mat, mld_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
implicit none
type(psb_lc_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_lc_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_lcspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_c_spmm_bld_inner
subroutine mld_c_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
import :: psb_lc_csr_sparse_mat, psb_lcspmat_type, psb_desc_type, &
& mld_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
@ -133,7 +147,7 @@ module mld_c_base_aggregator_mod
type(psb_lcspmat_type), intent(inout) :: op_prol, op_restr
type(psb_lcspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_c_spmm_bld_inner
end subroutine mld_c_old_spmm_bld_inner
end interface mld_spmm_bld_inner
contains

4
mlprec/mld_c_inner_mod.f90
Unescape Escape View File

@ -98,8 +98,8 @@ module mld_c_inner_mod
interface mld_map_to_tprol
subroutine mld_c_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type
use mld_c_prec_type, only : mld_c_onelev_type
import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lcspmat_type
import :: mld_c_onelev_type
implicit none
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)

4
mlprec/mld_c_jac_smoother.f90
Unescape Escape View File

@ -60,7 +60,7 @@ module mld_c_jac_smoother
!
type(psb_cspmat_type), pointer :: pa => null()
type(psb_cspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: dump => mld_c_jac_smoother_dmp
procedure, pass(sm) :: build => mld_c_jac_smoother_bld
@ -235,7 +235,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

3
mlprec/mld_c_onelev_mod.f90
Unescape Escape View File

@ -141,7 +141,6 @@ module mld_c_onelev_mod
type mld_cmlprec_wrk_type
complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_c_vect_type) :: vtx, vty, vx2l, vy2l
integer(psb_ipk_) :: wvsz = 0
type(psb_c_vect_type), allocatable :: wv(:)
contains
procedure, pass(wk) :: alloc => c_wrk_alloc
@ -420,7 +419,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = 0
val = psb_sizeof_ip+psb_sizeof_lp
val = val + lv%desc_ac%sizeof()
val = val + lv%ac%sizeof()
val = val + lv%tprol%sizeof()

33
mlprec/mld_c_prec_type.f90
Unescape Escape View File

@ -713,22 +713,25 @@ contains
end subroutine mld_c_apply1v
subroutine mld_c_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,&
subroutine mld_c_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,&
& global_num)
implicit none
class(mld_cprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
integer(psb_ipk_) :: i, j, il1, iln, lev
integer(psb_ipk_) :: icontxt, iam, np, iproc_
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_
info = 0
icontxt = prec%ictxt
call psb_info(icontxt,iam,np)
iln = size(prec%precv)
if (present(istart)) then
@ -739,13 +742,18 @@ contains
if (present(iend)) then
iln = min(iln, iend)
end if
iproc_ = -1
if (present(iproc)) then
iproc_ = iproc
end if
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
end if
end subroutine mld_c_dump
subroutine mld_c_cnv(prec,info,amold,vmold,imold)
@ -788,7 +796,7 @@ contains
class(psb_cprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info
! Local vars
integer(psb_ipk_) :: i, j, il1, ln, lname, lev
integer(psb_ipk_) :: i, j, ln, lev
integer(psb_ipk_) :: icontxt,iam, np
info = psb_success_
@ -892,8 +900,7 @@ contains
if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='complex(psb_spk_)')
goto 9999
end if
end do

5
mlprec/mld_d_as_smoother.f90
Unescape Escape View File

@ -67,7 +67,8 @@ module mld_d_as_smoother
!
type(psb_dspmat_type) :: nd
type(psb_desc_type) :: desc_data
integer(psb_ipk_) :: novr, restr, prol, nd_nnz_tot
integer(psb_ipk_) :: novr, restr, prol
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: check => mld_d_as_smoother_check
procedure, pass(sm) :: dump => mld_d_as_smoother_dmp
@ -321,7 +322,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = 3*psb_sizeof_ip + psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

20
mlprec/mld_d_base_aggregator_mod.f90
Unescape Escape View File

@ -44,7 +44,7 @@ module mld_d_base_aggregator_mod
use mld_base_prec_type, only : mld_dml_parms, mld_daggr_data
use psb_base_mod, only : psb_dspmat_type, psb_ldspmat_type, psb_d_vect_type, &
& psb_d_base_vect_type, psb_dlinmap_type, psb_dpk_, psb_ld_csr_sparse_mat, &
& psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_ld_coo_sparse_mat, psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_erractionsave, psb_error_handler, psb_success_, psb_toupper
!
!
@ -119,7 +119,21 @@ module mld_d_base_aggregator_mod
end interface
interface mld_spmm_bld_inner
subroutine mld_d_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
subroutine mld_d_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
import :: psb_ld_csr_sparse_mat, psb_ldspmat_type, psb_desc_type, &
& psb_ld_coo_sparse_mat, mld_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
implicit none
type(psb_ld_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_ld_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_ldspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_d_spmm_bld_inner
subroutine mld_d_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
import :: psb_ld_csr_sparse_mat, psb_ldspmat_type, psb_desc_type, &
& mld_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
@ -133,7 +147,7 @@ module mld_d_base_aggregator_mod
type(psb_ldspmat_type), intent(inout) :: op_prol, op_restr
type(psb_ldspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_d_spmm_bld_inner
end subroutine mld_d_old_spmm_bld_inner
end interface mld_spmm_bld_inner
contains

4
mlprec/mld_d_inner_mod.f90
Unescape Escape View File

@ -98,8 +98,8 @@ module mld_d_inner_mod
interface mld_map_to_tprol
subroutine mld_d_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type
use mld_d_prec_type, only : mld_d_onelev_type
import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_ldspmat_type
import :: mld_d_onelev_type
implicit none
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)

4
mlprec/mld_d_jac_smoother.f90
Unescape Escape View File

@ -60,7 +60,7 @@ module mld_d_jac_smoother
!
type(psb_dspmat_type), pointer :: pa => null()
type(psb_dspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: dump => mld_d_jac_smoother_dmp
procedure, pass(sm) :: build => mld_d_jac_smoother_bld
@ -235,7 +235,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

3
mlprec/mld_d_onelev_mod.f90
Unescape Escape View File

@ -141,7 +141,6 @@ module mld_d_onelev_mod
type mld_dmlprec_wrk_type
real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l
integer(psb_ipk_) :: wvsz = 0
type(psb_d_vect_type), allocatable :: wv(:)
contains
procedure, pass(wk) :: alloc => d_wrk_alloc
@ -420,7 +419,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = 0
val = psb_sizeof_ip+psb_sizeof_lp
val = val + lv%desc_ac%sizeof()
val = val + lv%ac%sizeof()
val = val + lv%tprol%sizeof()

33
mlprec/mld_d_prec_type.f90
Unescape Escape View File

@ -713,22 +713,25 @@ contains
end subroutine mld_d_apply1v
subroutine mld_d_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,&
subroutine mld_d_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,&
& global_num)
implicit none
class(mld_dprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
integer(psb_ipk_) :: i, j, il1, iln, lev
integer(psb_ipk_) :: icontxt, iam, np, iproc_
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_
info = 0
icontxt = prec%ictxt
call psb_info(icontxt,iam,np)
iln = size(prec%precv)
if (present(istart)) then
@ -739,13 +742,18 @@ contains
if (present(iend)) then
iln = min(iln, iend)
end if
iproc_ = -1
if (present(iproc)) then
iproc_ = iproc
end if
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
end if
end subroutine mld_d_dump
subroutine mld_d_cnv(prec,info,amold,vmold,imold)
@ -788,7 +796,7 @@ contains
class(psb_dprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info
! Local vars
integer(psb_ipk_) :: i, j, il1, ln, lname, lev
integer(psb_ipk_) :: i, j, ln, lev
integer(psb_ipk_) :: icontxt,iam, np
info = psb_success_
@ -892,8 +900,7 @@ contains
if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='real(psb_dpk_)')
goto 9999
end if
end do

5
mlprec/mld_s_as_smoother.f90
Unescape Escape View File

@ -67,7 +67,8 @@ module mld_s_as_smoother
!
type(psb_sspmat_type) :: nd
type(psb_desc_type) :: desc_data
integer(psb_ipk_) :: novr, restr, prol, nd_nnz_tot
integer(psb_ipk_) :: novr, restr, prol
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: check => mld_s_as_smoother_check
procedure, pass(sm) :: dump => mld_s_as_smoother_dmp
@ -321,7 +322,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = 3*psb_sizeof_ip + psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

20
mlprec/mld_s_base_aggregator_mod.f90
Unescape Escape View File

@ -44,7 +44,7 @@ module mld_s_base_aggregator_mod
use mld_base_prec_type, only : mld_sml_parms, mld_saggr_data
use psb_base_mod, only : psb_sspmat_type, psb_lsspmat_type, psb_s_vect_type, &
& psb_s_base_vect_type, psb_slinmap_type, psb_spk_, psb_ls_csr_sparse_mat, &
& psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_ls_coo_sparse_mat, psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_erractionsave, psb_error_handler, psb_success_, psb_toupper
!
!
@ -119,7 +119,21 @@ module mld_s_base_aggregator_mod
end interface
interface mld_spmm_bld_inner
subroutine mld_s_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
subroutine mld_s_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
import :: psb_ls_csr_sparse_mat, psb_lsspmat_type, psb_desc_type, &
& psb_ls_coo_sparse_mat, mld_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
implicit none
type(psb_ls_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_sml_parms), intent(inout) :: parms
type(psb_ls_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_lsspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_s_spmm_bld_inner
subroutine mld_s_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
import :: psb_ls_csr_sparse_mat, psb_lsspmat_type, psb_desc_type, &
& mld_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
@ -133,7 +147,7 @@ module mld_s_base_aggregator_mod
type(psb_lsspmat_type), intent(inout) :: op_prol, op_restr
type(psb_lsspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_s_spmm_bld_inner
end subroutine mld_s_old_spmm_bld_inner
end interface mld_spmm_bld_inner
contains

4
mlprec/mld_s_inner_mod.f90
Unescape Escape View File

@ -98,8 +98,8 @@ module mld_s_inner_mod
interface mld_map_to_tprol
subroutine mld_s_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type
use mld_s_prec_type, only : mld_s_onelev_type
import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_, psb_lsspmat_type
import :: mld_s_onelev_type
implicit none
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)

4
mlprec/mld_s_jac_smoother.f90
Unescape Escape View File

@ -60,7 +60,7 @@ module mld_s_jac_smoother
!
type(psb_sspmat_type), pointer :: pa => null()
type(psb_sspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: dump => mld_s_jac_smoother_dmp
procedure, pass(sm) :: build => mld_s_jac_smoother_bld
@ -235,7 +235,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

3
mlprec/mld_s_onelev_mod.f90
Unescape Escape View File

@ -141,7 +141,6 @@ module mld_s_onelev_mod
type mld_smlprec_wrk_type
real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_s_vect_type) :: vtx, vty, vx2l, vy2l
integer(psb_ipk_) :: wvsz = 0
type(psb_s_vect_type), allocatable :: wv(:)
contains
procedure, pass(wk) :: alloc => s_wrk_alloc
@ -420,7 +419,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = 0
val = psb_sizeof_ip+psb_sizeof_lp
val = val + lv%desc_ac%sizeof()
val = val + lv%ac%sizeof()
val = val + lv%tprol%sizeof()

33
mlprec/mld_s_prec_type.f90
Unescape Escape View File

@ -713,22 +713,25 @@ contains
end subroutine mld_s_apply1v
subroutine mld_s_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,&
subroutine mld_s_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,&
& global_num)
implicit none
class(mld_sprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
integer(psb_ipk_) :: i, j, il1, iln, lev
integer(psb_ipk_) :: icontxt, iam, np, iproc_
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_
info = 0
icontxt = prec%ictxt
call psb_info(icontxt,iam,np)
iln = size(prec%precv)
if (present(istart)) then
@ -739,13 +742,18 @@ contains
if (present(iend)) then
iln = min(iln, iend)
end if
iproc_ = -1
if (present(iproc)) then
iproc_ = iproc
end if
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
end if
end subroutine mld_s_dump
subroutine mld_s_cnv(prec,info,amold,vmold,imold)
@ -788,7 +796,7 @@ contains
class(psb_sprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info
! Local vars
integer(psb_ipk_) :: i, j, il1, ln, lname, lev
integer(psb_ipk_) :: i, j, ln, lev
integer(psb_ipk_) :: icontxt,iam, np
info = psb_success_
@ -892,8 +900,7 @@ contains
if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='real(psb_spk_)')
goto 9999
end if
end do

5
mlprec/mld_z_as_smoother.f90
Unescape Escape View File

@ -67,7 +67,8 @@ module mld_z_as_smoother
!
type(psb_zspmat_type) :: nd
type(psb_desc_type) :: desc_data
integer(psb_ipk_) :: novr, restr, prol, nd_nnz_tot
integer(psb_ipk_) :: novr, restr, prol
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: check => mld_z_as_smoother_check
procedure, pass(sm) :: dump => mld_z_as_smoother_dmp
@ -321,7 +322,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = 3*psb_sizeof_ip + psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

20
mlprec/mld_z_base_aggregator_mod.f90
Unescape Escape View File

@ -44,7 +44,7 @@ module mld_z_base_aggregator_mod
use mld_base_prec_type, only : mld_dml_parms, mld_daggr_data
use psb_base_mod, only : psb_zspmat_type, psb_lzspmat_type, psb_z_vect_type, &
& psb_z_base_vect_type, psb_zlinmap_type, psb_dpk_, psb_lz_csr_sparse_mat, &
& psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_lz_coo_sparse_mat, psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
& psb_erractionsave, psb_error_handler, psb_success_, psb_toupper
!
!
@ -119,7 +119,21 @@ module mld_z_base_aggregator_mod
end interface
interface mld_spmm_bld_inner
subroutine mld_z_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
subroutine mld_z_spmm_bld_inner(a_csr,desc_a,nlaggr,parms,ac,&
& coo_prol,desc_cprol,coo_restr,info)
import :: psb_lz_csr_sparse_mat, psb_lzspmat_type, psb_desc_type, &
& psb_lz_coo_sparse_mat, mld_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
implicit none
type(psb_lz_csr_sparse_mat), intent(inout) :: a_csr
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), intent(inout) :: nlaggr(:)
type(mld_dml_parms), intent(inout) :: parms
type(psb_lz_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
type(psb_desc_type), intent(inout) :: desc_cprol
type(psb_lzspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_z_spmm_bld_inner
subroutine mld_z_old_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,ac,&
& op_prol,op_restr,info)
import :: psb_lz_csr_sparse_mat, psb_lzspmat_type, psb_desc_type, &
& mld_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
@ -133,7 +147,7 @@ module mld_z_base_aggregator_mod
type(psb_lzspmat_type), intent(inout) :: op_prol, op_restr
type(psb_lzspmat_type), intent(out) :: ac
integer(psb_ipk_), intent(out) :: info
end subroutine mld_z_spmm_bld_inner
end subroutine mld_z_old_spmm_bld_inner
end interface mld_spmm_bld_inner
contains

4
mlprec/mld_z_inner_mod.f90
Unescape Escape View File

@ -98,8 +98,8 @@ module mld_z_inner_mod
interface mld_map_to_tprol
subroutine mld_z_map_to_tprol(desc_a,ilaggr,nlaggr,op_prol,info)
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type
use mld_z_prec_type, only : mld_z_onelev_type
import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_, psb_lzspmat_type
import :: mld_z_onelev_type
implicit none
type(psb_desc_type), intent(in) :: desc_a
integer(psb_lpk_), allocatable, intent(inout) :: ilaggr(:),nlaggr(:)

4
mlprec/mld_z_jac_smoother.f90
Unescape Escape View File

@ -60,7 +60,7 @@ module mld_z_jac_smoother
!
type(psb_zspmat_type), pointer :: pa => null()
type(psb_zspmat_type) :: nd
integer(psb_ipk_) :: nnz_nd_tot
integer(psb_lpk_) :: nd_nnz_tot
contains
procedure, pass(sm) :: dump => mld_z_jac_smoother_dmp
procedure, pass(sm) :: build => mld_z_jac_smoother_bld
@ -235,7 +235,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = psb_sizeof_ip
val = psb_sizeof_lp
if (allocated(sm%sv)) val = val + sm%sv%sizeof()
val = val + sm%nd%sizeof()

3
mlprec/mld_z_onelev_mod.f90
Unescape Escape View File

@ -141,7 +141,6 @@ module mld_z_onelev_mod
type mld_zmlprec_wrk_type
complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_z_vect_type) :: vtx, vty, vx2l, vy2l
integer(psb_ipk_) :: wvsz = 0
type(psb_z_vect_type), allocatable :: wv(:)
contains
procedure, pass(wk) :: alloc => z_wrk_alloc
@ -420,7 +419,7 @@ contains
integer(psb_epk_) :: val
integer(psb_ipk_) :: i
val = 0
val = psb_sizeof_ip+psb_sizeof_lp
val = val + lv%desc_ac%sizeof()
val = val + lv%ac%sizeof()
val = val + lv%tprol%sizeof()

33
mlprec/mld_z_prec_type.f90
Unescape Escape View File

@ -713,22 +713,25 @@ contains
end subroutine mld_z_apply1v
subroutine mld_z_dump(prec,info,istart,iend,prefix,head,ac,rp,smoother,solver,tprol,&
subroutine mld_z_dump(prec,info,istart,iend,iproc,prefix,head,&
& ac,rp,smoother,solver,tprol,&
& global_num)
implicit none
class(mld_zprec_type), intent(in) :: prec
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: istart, iend
integer(psb_ipk_), intent(in), optional :: istart, iend, iproc
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: smoother, solver,ac, rp, tprol, global_num
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
integer(psb_ipk_) :: i, j, il1, iln, lev
integer(psb_ipk_) :: icontxt, iam, np, iproc_
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
! len of prefix_
info = 0
icontxt = prec%ictxt
call psb_info(icontxt,iam,np)
iln = size(prec%precv)
if (present(istart)) then
@ -739,13 +742,18 @@ contains
if (present(iend)) then
iln = min(iln, iend)
end if
iproc_ = -1
if (present(iproc)) then
iproc_ = iproc
end if
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
if ((iproc_ == -1).or.(iproc_==iam)) then
do lev=il1, iln
call prec%precv(lev)%dump(lev,info,prefix=prefix,head=head,&
& ac=ac,smoother=smoother,solver=solver,rp=rp,tprol=tprol, &
& global_num=global_num)
end do
end if
end subroutine mld_z_dump
subroutine mld_z_cnv(prec,info,amold,vmold,imold)
@ -788,7 +796,7 @@ contains
class(psb_zprec_type), target, intent(inout) :: precout
integer(psb_ipk_), intent(out) :: info
! Local vars
integer(psb_ipk_) :: i, j, il1, ln, lname, lev
integer(psb_ipk_) :: i, j, ln, lev
integer(psb_ipk_) :: icontxt,iam, np
info = psb_success_
@ -892,8 +900,7 @@ contains
if (psb_errstatus_fatal()) then
nc2l = prec%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
call psb_errpush(info,name,i_err=(/2*nc2l/), a_err='complex(psb_dpk_)')
goto 9999
end if
end do

Write Preview
Loading…
Cancel
Save