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

mld2p4-2:

Browse Source 
mlprec/impl/mld_caggrmat_asb.f90
 mlprec/impl/mld_caggrmat_minnrg_asb.F90
 mlprec/impl/mld_caggrmat_nosmth_asb.F90
 mlprec/impl/mld_caggrmat_smth_asb.F90
 mlprec/impl/mld_daggrmat_asb.f90
 mlprec/impl/mld_daggrmat_minnrg_asb.F90
 mlprec/impl/mld_daggrmat_nosmth_asb.F90
 mlprec/impl/mld_daggrmat_smth_asb.F90
 mlprec/impl/mld_saggrmat_asb.f90
 mlprec/impl/mld_saggrmat_minnrg_asb.F90
 mlprec/impl/mld_saggrmat_nosmth_asb.F90
 mlprec/impl/mld_saggrmat_smth_asb.F90
 mlprec/impl/mld_zaggrmat_asb.f90
 mlprec/impl/mld_zaggrmat_minnrg_asb.F90
 mlprec/impl/mld_zaggrmat_nosmth_asb.F90
 mlprec/impl/mld_zaggrmat_smth_asb.F90
 tests/pdegen/runs/ppde.inp

Unified matrix aggregation routines.
stopcriterion
Salvatore Filippone 13 years ago
parent ad88b4298d
commit 34581e785e
17 changed files with 222 additions and 224 deletions
Show Stats Download Patch File Download Diff File

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

@ -113,7 +113,7 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer :: ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
character(len=20) :: name
name='mld_aggrmat_asb'
@ -122,7 +122,6 @@ subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)

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

@ -59,7 +59,7 @@
! of A, and omega is a suitable smoothing parameter. An estimate of the spectral
! radius of D^(-1)A, to be used in the computation of omega, is provided,
! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_cprecinit and mld_zprecset.
! through mld_cprecinit and mld_cprecset.
!
! This routine can also build A_C according to a "bizarre" aggregation algorithm,
! using a "naive" prolongator proposed by the authors of MLD2P4. However, this
@ -119,7 +119,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Local variables
type(psb_cspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm
@ -133,6 +133,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
complex(psb_spk_), allocatable :: adiag(:), adinv(:)
complex(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(psb_spk_) :: anorm, theta
@ -162,9 +163,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -188,9 +188,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),&
& a_err='real(psb_spk_)')
info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr;
call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_spk_)')
goto 9999
end if
@ -236,7 +235,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies cone.'
& ' Initial copies done.'
call da%scal(adinv,info)
@ -292,7 +291,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
end do
if(min(real(omf(i)),aimag(omf(i))) .lt. dzero) omf(i) = zzero
!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero
if(psb_minreal(omf(i)) < szero) omf(i) = czero
end do
omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
@ -411,7 +411,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Ok, let's start over with the restrictor
!
call ptilde%transp(rtilde)
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
@ -424,7 +424,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transp(atmp2)
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,1,nrow,1,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
@ -473,7 +473,8 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
do j= acsc%icp(i),acsc%icp(i+1)-1
if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
end do
if(min(real(omf(i)),aimag(omf(i))) .lt. dzero) omf(i) = zzero
!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero
if(psb_minreal(omf(i)) < szero) omf(i) = czero
end do
omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
call psb_halo(omf,desc_a,info)
@ -671,14 +672,14 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,tmpcoo%val,nzbr,idisp,&
& mpi_double_complex,icomm,info)
call mpi_allgatherv(bcoo%val,ndx,mpi_complex,tmpcoo%val,nzbr,idisp,&
& mpi_complex,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

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

@ -101,13 +101,15 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer ::ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name
type(psb_cspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5)
type(psb_cspmat_type) :: am1,am2
type(psb_c_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, &
& naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb'
@ -128,9 +130,8 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -201,12 +202,12 @@ subroutine mld_caggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,ac_coo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
call mpi_allgatherv(bcoo%val,ndx,mpi_complex,ac_coo%val,nzbr,idisp,&
& mpi_complex,icomm,minfo)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
if(info /= psb_success_) then
info=-1
call psb_errpush(info,name)

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

@ -122,12 +122,13 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt,np,me, err_act, icomm
integer ::ictxt, np, me, err_act
character(len=20) :: name
type(psb_cspmat_type) :: am1,am2, am3, am4
type(psb_c_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
@ -141,7 +142,6 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -157,9 +157,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -187,9 +186,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(adiag(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),&
& a_err='complex(psb_spk_)')
info=psb_err_alloc_request_; ierr(1)=nrow;
call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_spk_)')
goto 9999
end if
@ -519,7 +517,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac')
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &

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

@ -52,7 +52,7 @@
! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine.
! The prolongator P_C is built here from this mapping, according to the
! value of p%iprcparm(mld_aggr_kind_), specified by the user through
! mld_dprecinit and mld_dprecset.
! mld_dprecinit and mld_zprecset.
!
! Currently three different prolongators are implemented, corresponding to
! three aggregation algorithms:
@ -78,7 +78,7 @@
!
!
! Arguments:
! a - type(psb_dspmat_type), input.
! a - type(psb_cspmat_type), input.
! The sparse matrix structure containing the local part of
! the fine-level matrix.
! desc_a - type(psb_desc_type), input.
@ -113,7 +113,7 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer :: ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
character(len=20) :: name
name='mld_aggrmat_asb'
@ -122,7 +122,6 @@ subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)

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

@ -68,7 +68,7 @@
!
! The coarse-level matrix A_C is distributed among the parallel processes or
! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_dprecinit and mld_dprecset.
! specified by the user through mld_dprecinit and mld_zprecset.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
@ -119,23 +119,25 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Local variables
type(psb_dspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_dspmat_type) :: am1, am2, af, ptilde, rtilde, atran, atp, atdatp
type(psb_dspmat_type) :: am3, am4, ap, adap, atmp, rada, ra, atmp2, dap, dadap, da
type(psb_dspmat_type) :: am1,am2, af, ptilde, rtilde, atran, atp, atdatp
type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_dspmat_type) :: dat, datp, datdatp, atmp3
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf
type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
real(psb_dpk_), allocatable :: adiag(:), omf(:),omp(:),omi(:),&
& oden(:), adinv(:)
real(psb_dpk_), allocatable :: adiag(:), adinv(:)
real(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(psb_dpk_) :: omega, anorm, tmp, dg, theta, alpha,beta, ommx
real(psb_dpk_) :: anorm, theta
real(psb_dpk_) :: tmp, alpha, beta, ommx
name='mld_aggrmat_minnrg'
if(psb_get_errstatus().ne.0) return
@ -161,9 +163,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -187,9 +188,8 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),&
& a_err='real(psb_dpk_)')
info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr;
call psb_errpush(info,name,i_err=ierr,a_err='real(psb_dpk_)')
goto 9999
end if
@ -249,14 +249,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call dap%clone(atmp,info)
if (info == psb_success_) call psb_sphalo(atmp,desc_a,am4,info,&
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.,outfmt='CSR ')
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)
if (info == psb_success_) call am4%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sphalo 2')
goto 9999
end if
call psb_symbmm(da,atmp,dadap,info)
call psb_numbmm(da,atmp,dadap)
@ -284,18 +280,19 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call am3%mv_to(acsr3)
! Compute omega_int
ommx = -huge(done)
ommx = cmplx(szero,szero)
do i=1, ncol
omi(i) = omp(ilaggr(i))
ommx = max(ommx,omi(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
end do
! Compute omega_fine
do i=1, nrow
omf(i) = ommx
do j=acsr3%irp(i),acsr3%irp(i+1)-1
omf(i) = min(omf(i),omi(acsr3%ja(j)))
if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
end do
omf(i) = max(dzero,omf(i))
!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero
if(psb_minreal(omf(i)) < dzero) omf(i) = dzero
end do
omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
@ -414,7 +411,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Ok, let's start over with the restrictor
!
call ptilde%transp(rtilde)
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
@ -427,7 +424,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transp(atmp2)
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,1,nrow,1,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
@ -461,10 +458,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = -huge(done)
ommx = cmplx(szero,szero)
do i=1, ncol
omi(i) = omp(ilaggr(i))
ommx = max(ommx,omi(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
end do
! Compute omega_fine
! Going over the columns of atmp means going over the rows
@ -474,9 +471,10 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
do i=1, nrow
omf(i) = ommx
do j= acsc%icp(i),acsc%icp(i+1)-1
omf(i) = min(omf(i),omi(acsc%ia(j)))
if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
end do
omf(i) = max(dzero,omf(i))
!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero
if(psb_minreal(omf(i)) < dzero) omf(i) = dzero
end do
omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
call psb_halo(omf,desc_a,info)
@ -629,14 +627,14 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call p%ac%cscnv(info,type='csr')
if (np>1) then
call am1%mv_to(acsr1)
nzl = acsr1%get_nzeros()
call psb_glob_to_loc(acsr1%ja(1:nzl),p%desc_ac,info,'I')
call am1%mv_to(acsr)
nzl = acsr%get_nzeros()
call psb_glob_to_loc(acsr%ja(1:nzl),p%desc_ac,info,'I')
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_glob_to_loc')
goto 9999
end if
call am1%mv_from(acsr1)
call am1%mv_from(acsr)
endif
call am1%set_ncols(p%desc_ac%get_local_cols())
@ -677,11 +675,11 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,tmpcoo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -817,7 +815,7 @@ contains
if (ip1 > nv1) exit
if (ip2 > nv2) exit
if (iv1(ip1) == iv2(ip2)) then
dot = dot + v1(ip1)*v2(ip2)
dot = dot + (v1(ip1))*v2(ip2)
ip1 = ip1 + 1
ip2 = ip2 + 1
else if (iv1(ip1) < iv2(ip2)) then

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

@ -51,7 +51,7 @@
!
! The coarse-level matrix A_C is distributed among the parallel processes or
! replicated on each of them, according to the value of p%parms%coarse_mat
! specified by the user through mld_dprecinit and mld_dprecset.
! specified by the user through mld_dprecinit and mld_zprecset.
!
! For details see
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
@ -101,13 +101,15 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer ::ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name
type(psb_dspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5)
type(psb_dspmat_type) :: am1,am2
type(psb_d_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, &
& naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb'
@ -128,9 +130,8 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -202,11 +203,11 @@ subroutine mld_daggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,ac_coo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& mpi_double_precision,icomm,minfo)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
if(info /= psb_success_) then
info=-1
call psb_errpush(info,name)

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

@ -59,7 +59,7 @@
! of A, and omega is a suitable smoothing parameter. An estimate of the spectral
! radius of D^(-1)A, to be used in the computation of omega, is provided,
! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_dprecinit and mld_dprecset.
! through mld_dprecinit and mld_zprecset.
!
! This routine can also build A_C according to a "bizarre" aggregation algorithm,
! using a "naive" prolongator proposed by the authors of MLD2P4. However, this
@ -68,7 +68,7 @@
!
! The coarse-level matrix A_C is distributed among the parallel processes or
! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_dprecinit and mld_dprecset.
! specified by the user through mld_dprecinit and mld_zprecset.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
@ -122,16 +122,17 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt,np,me, err_act, icomm
integer ::ictxt, np, me, err_act
character(len=20) :: name
type(psb_dspmat_type) :: am1,am2, am3, am4
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(psb_dpk_) :: omega, anorm, tmp, dg, theta
real(psb_dpk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_smth_asb'
if(psb_get_errstatus().ne.0) return
@ -141,7 +142,6 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -157,9 +157,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -187,9 +186,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(adiag(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),&
& a_err='real(psb_dpk_)')
info=psb_err_alloc_request_; ierr(1)=nrow;
call psb_errpush(info,name,i_err=ierr,a_err='real(psb_dpk_)')
goto 9999
end if
@ -228,7 +226,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies done.'
& ' Initial copies sone.'
if (filter_mat) then
!
@ -296,7 +294,7 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = dzero
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
@ -519,7 +517,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac')
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &

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

@ -52,7 +52,7 @@
! adjacency graph of A_C has been computed by the mld_aggrmap_bld subroutine.
! The prolongator P_C is built here from this mapping, according to the
! value of p%iprcparm(mld_aggr_kind_), specified by the user through
! mld_sprecinit and mld_dprecset.
! mld_sprecinit and mld_zprecset.
!
! Currently three different prolongators are implemented, corresponding to
! three aggregation algorithms:
@ -78,7 +78,7 @@
!
!
! Arguments:
! a - type(psb_sspmat_type), input.
! a - type(psb_cspmat_type), input.
! The sparse matrix structure containing the local part of
! the fine-level matrix.
! desc_a - type(psb_desc_type), input.
@ -113,7 +113,7 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer :: ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
character(len=20) :: name
name='mld_aggrmat_asb'
@ -122,7 +122,6 @@ subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)

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

@ -59,7 +59,7 @@
! of A, and omega is a suitable smoothing parameter. An estimate of the spectral
! radius of D^(-1)A, to be used in the computation of omega, is provided,
! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_sprecinit and mld_dprecset.
! through mld_sprecinit and mld_sprecset.
!
! This routine can also build A_C according to a "bizarre" aggregation algorithm,
! using a "naive" prolongator proposed by the authors of MLD2P4. However, this
@ -68,7 +68,7 @@
!
! The coarse-level matrix A_C is distributed among the parallel processes or
! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_sprecinit and mld_dprecset.
! specified by the user through mld_sprecinit and mld_zprecset.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
@ -119,7 +119,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Local variables
type(psb_sspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm
@ -130,12 +130,14 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf
type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
real(psb_spk_), allocatable :: adiag(:), omf(:),omp(:),omi(:),&
& oden(:), adinv(:)
real(psb_spk_), allocatable :: adiag(:), adinv(:)
real(psb_spk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(psb_spk_) :: omega, anorm, tmp, dg, theta, alpha,beta, ommx
real(psb_spk_) :: anorm, theta
real(psb_spk_) :: tmp, alpha, beta, ommx
name='mld_aggrmat_minnrg'
if(psb_get_errstatus().ne.0) return
@ -161,9 +163,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -187,9 +188,8 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),&
& a_err='real(psb_spk_)')
info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr;
call psb_errpush(info,name,i_err=ierr,a_err='real(psb_spk_)')
goto 9999
end if
@ -235,7 +235,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies sone.'
& ' Initial copies done.'
call da%scal(adinv,info)
@ -280,18 +280,19 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call am3%mv_to(acsr3)
! Compute omega_int
ommx = -huge(sone)
ommx = cmplx(szero,szero)
do i=1, ncol
omi(i) = omp(ilaggr(i))
ommx = max(ommx,omi(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
end do
! Compute omega_fine
do i=1, nrow
omf(i) = ommx
do j=acsr3%irp(i),acsr3%irp(i+1)-1
omf(i) = min(omf(i),omi(acsr3%ja(j)))
if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
end do
omf(i) = max(szero,omf(i))
!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero
if(psb_minreal(omf(i)) < szero) omf(i) = szero
end do
omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
@ -410,7 +411,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Ok, let's start over with the restrictor
!
call ptilde%transp(rtilde)
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
@ -423,7 +424,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transp(atmp2)
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,1,nrow,1,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
@ -457,10 +458,10 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = -huge(sone)
ommx = cmplx(szero,szero)
do i=1, ncol
omi(i) = omp(ilaggr(i))
ommx = max(ommx,omi(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
end do
! Compute omega_fine
! Going over the columns of atmp means going over the rows
@ -470,9 +471,10 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
do i=1, nrow
omf(i) = ommx
do j= acsc%icp(i),acsc%icp(i+1)-1
omf(i) = min(omf(i),omi(acsc%ia(j)))
if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
end do
omf(i) = max(szero,omf(i))
!!$ if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero
if(psb_minreal(omf(i)) < szero) omf(i) = szero
end do
omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
call psb_halo(omf,desc_a,info)
@ -670,14 +672,14 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,tmpcoo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(bcoo%val,ndx,mpi_real,tmpcoo%val,nzbr,idisp,&
& mpi_real,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@ -813,7 +815,7 @@ contains
if (ip1 > nv1) exit
if (ip2 > nv2) exit
if (iv1(ip1) == iv2(ip2)) then
dot = dot + v1(ip1)*v2(ip2)
dot = dot + (v1(ip1))*v2(ip2)
ip1 = ip1 + 1
ip2 = ip2 + 1
else if (iv1(ip1) < iv2(ip2)) then

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

@ -51,7 +51,7 @@
!
! The coarse-level matrix A_C is distributed among the parallel processes or
! replicated on each of them, according to the value of p%parms%coarse_mat
! specified by the user through mld_sprecinit and mld_dprecset.
! specified by the user through mld_sprecinit and mld_zprecset.
!
! For details see
! P. D'Ambra, D. di Serafino and S. Filippone, On the development of
@ -101,13 +101,15 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer ::ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name
type(psb_sspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5)
type(psb_sspmat_type) :: am1,am2
type(psb_s_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, &
& naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb'
@ -128,9 +130,8 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -201,12 +202,12 @@ subroutine mld_saggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,ac_coo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
call mpi_allgatherv(bcoo%val,ndx,mpi_real,ac_coo%val,nzbr,idisp,&
& mpi_real,icomm,minfo)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
if(info /= psb_success_) then
info=-1
call psb_errpush(info,name)

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

@ -59,7 +59,7 @@
! of A, and omega is a suitable smoothing parameter. An estimate of the spectral
! radius of D^(-1)A, to be used in the computation of omega, is provided,
! according to the value of p%parms%aggr_omega_alg, specified by the user
! through mld_sprecinit and mld_dprecset.
! through mld_sprecinit and mld_zprecset.
!
! This routine can also build A_C according to a "bizarre" aggregation algorithm,
! using a "naive" prolongator proposed by the authors of MLD2P4. However, this
@ -68,7 +68,7 @@
!
! The coarse-level matrix A_C is distributed among the parallel processes or
! replicated on each of them, according to the value of p%parms%coarse_mat,
! specified by the user through mld_sprecinit and mld_dprecset.
! specified by the user through mld_sprecinit and mld_zprecset.
!
! For more details see
! M. Brezina and P. Vanek, A black-box iterative solver based on a
@ -122,16 +122,17 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt,np,me, err_act, icomm
integer ::ictxt, np, me, err_act
character(len=20) :: name
type(psb_sspmat_type) :: am1,am2, am3, am4
type(psb_s_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_spk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(psb_spk_) :: omega, anorm, tmp, dg, theta
real(psb_spk_) :: anorm, omega, tmp, dg, theta
name='mld_aggrmat_smth_asb'
if(psb_get_errstatus().ne.0) return
@ -141,7 +142,6 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -157,9 +157,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -187,9 +186,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(adiag(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),&
& a_err='real(psb_spk_)')
info=psb_err_alloc_request_; ierr(1)=nrow;
call psb_errpush(info,name,i_err=ierr,a_err='real(psb_spk_)')
goto 9999
end if
@ -519,7 +517,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac')
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &

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

@ -78,7 +78,7 @@
!
!
! Arguments:
! a - type(psb_zspmat_type), input.
! a - type(psb_cspmat_type), input.
! The sparse matrix structure containing the local part of
! the fine-level matrix.
! desc_a - type(psb_desc_type), input.
@ -113,7 +113,7 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer :: ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
character(len=20) :: name
name='mld_aggrmat_asb'
@ -122,7 +122,6 @@ subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info)
call psb_erractionsave(err_act)
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
call psb_info(ictxt, me, np)

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

@ -119,7 +119,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Local variables
type(psb_zspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrt
integer :: ictxt,np,me, err_act, icomm
@ -133,6 +133,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
complex(psb_dpk_), allocatable :: adiag(:), adinv(:)
complex(psb_dpk_), allocatable :: omf(:), omp(:), omi(:), oden(:)
logical :: filter_mat
integer(psb_ipk_) :: ierr(5)
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
real(psb_dpk_) :: anorm, theta
@ -162,9 +163,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -188,9 +188,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
& omf(ncol),omp(ntaggr),oden(ntaggr),omi(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/6*ncol+ntaggr,0,0,0,0/),&
& a_err='real(psb_dpk_)')
info=psb_err_alloc_request_; ierr(1)=6*ncol+ntaggr;
call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_dpk_)')
goto 9999
end if
@ -236,7 +235,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies zone.'
& ' Initial copies done.'
call da%scal(adinv,info)
@ -281,7 +280,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call am3%mv_to(acsr3)
! Compute omega_int
ommx = cmplx(dzero,dzero)
ommx = cmplx(szero,szero)
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@ -292,7 +291,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
end do
if(min(real(omf(i)),aimag(omf(i))) .lt. dzero) omf(i) = zzero
!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero
if(psb_minreal(omf(i)) < dzero) omf(i) = zzero
end do
omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
@ -411,7 +411,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
!
! Ok, let's start over with the restrictor
!
call ptilde%transp(rtilde)
call ptilde%transc(rtilde)
call a%cscnv(atmp,info,type='csr')
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
@ -424,7 +424,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transp(atmp2)
call atmp%transc(atmp2)
call atmp2%csclip(dat,info,1,nrow,1,ncol)
call dat%cscnv(info,type='csr')
call dat%scal(adinv,info)
@ -458,7 +458,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = cmplx(dzero,dzero)
ommx = cmplx(szero,szero)
do i=1, ncol
omi(i) = omp(ilaggr(i))
if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
@ -473,7 +473,8 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
do j= acsc%icp(i),acsc%icp(i+1)-1
if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
end do
if(min(real(omf(i)),aimag(omf(i))) .lt. dzero) omf(i) = zzero
!!$ if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero
if(psb_minreal(omf(i)) < dzero) omf(i) = zzero
end do
omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
call psb_halo(omf,desc_a,info)
@ -674,11 +675,11 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,tmpcoo%val,nzbr,idisp,&
& mpi_double_complex,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,tmpcoo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info == psb_success_)&
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,tmpcoo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

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

@ -101,13 +101,15 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, intent(out) :: info
! Local variables
integer ::ictxt,np,me, err_act, icomm
integer :: ictxt,np,me, err_act
integer(psb_mpik_) :: icomm, ndx, minfo
character(len=20) :: name
type(psb_zspmat_type) :: b
integer, allocatable :: nzbr(:), idisp(:)
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5)
type(psb_zspmat_type) :: am1,am2
type(psb_z_coo_sparse_mat) :: acoo1, acoo2, bcoo, ac_coo
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, &
& naggr, nzt, naggrm1, i
name='mld_aggrmat_nosmth_asb'
@ -128,9 +130,8 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
ntaggr = sum(nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -202,11 +203,11 @@ subroutine mld_zaggrmat_nosmth_asb(a,desc_a,ilaggr,nlaggr,p,info)
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_complex,ac_coo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_integer,icomm,info)
& mpi_double_complex,icomm,minfo)
call mpi_allgatherv(bcoo%ia,ndx,psb_mpi_ipk_integer,ac_coo%ia,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
call mpi_allgatherv(bcoo%ja,ndx,psb_mpi_ipk_integer,ac_coo%ja,nzbr,idisp,&
& psb_mpi_ipk_integer,icomm,minfo)
if(info /= psb_success_) then
info=-1
call psb_errpush(info,name)

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

@ -122,12 +122,13 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer, allocatable :: nzbr(:), idisp(:)
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
integer ::ictxt,np,me, err_act, icomm
integer ::ictxt, np, me, err_act
character(len=20) :: name
type(psb_zspmat_type) :: am1,am2, am3, am4
type(psb_z_coo_sparse_mat) :: acoo, acoof, bcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_dpk_), allocatable :: adiag(:)
integer(psb_ipk_) :: ierr(5)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
integer, parameter :: ncmax=16
@ -141,7 +142,6 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
icomm = desc_a%get_mpic()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
@ -157,9 +157,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
& a_err='integer')
info=psb_err_alloc_request_; ierr(1)=2*np;
call psb_errpush(info,name,i_err=ierr,a_err='integer')
goto 9999
end if
@ -187,9 +186,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
allocate(adiag(ncol),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/nrow,0,0,0,0/),&
& a_err='complex(psb_dpk_)')
info=psb_err_alloc_request_; ierr(1)=nrow;
call psb_errpush(info,name,i_err=ierr,a_err='complex(psb_dpk_)')
goto 9999
end if
@ -228,7 +226,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies done.'
& ' Initial copies sone.'
if (filter_mat) then
!
@ -296,7 +294,7 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
dg = done
nrw = acsr3%get_nrows()
do i=1, nrw
tmp = dzero
tmp = szero
do j=acsr3%irp(i),acsr3%irp(i+1)-1
if (acsr3%ja(j) <= nrw) then
tmp = tmp + abs(acsr3%val(j))
@ -519,7 +517,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac')
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &

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

@ -5,16 +5,16 @@ CSR ! Storage format CSR COO JAD
0100 ! ITMAX
-1 ! ITRACE
30 ! IRST (restart for RGMRES and BiCGSTABL)
1.d-9 ! EPS
3L-M-RAS-I-D4 ! Longer descriptive name for preconditioner (up to 20 chars)
1.d-6 ! EPS
3L-MUL-RAS-BJAC4-ILU ! Descriptive name for preconditioner (up to 40 chars)
ML ! Preconditioner NONE JACOBI BJAC AS ML
0 ! Number of overlap layers for AS preconditioner at finest level
1 ! Number of overlap layers for AS preconditioner at finest level
HALO ! Restriction operator NONE HALO
NONE ! Prolongation operator NONE SUM AVG
ILU ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU
0 ! Level-set N for ILU(N)
1.d-5 ! Threshold T for ILU(T,P)
1 ! Smoother/Jacobi sweeps
0 ! Level-set N for ILU(N), and P for ILUT
1.d-4 ! Threshold T for ILU(T,P)
4 ! Smoother/Jacobi sweeps
BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML
3 ! Number of levels in a multilevel preconditioner
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED

Write Preview
Loading…
Cancel
Save