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

mld2p4:

Browse Source 
mlprec/mld_caggrmat_smth_asb.F90
 mlprec/mld_d_as_smoother.f90
 mlprec/mld_d_ilu_solver.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_d_jac_smoother.f90
 mlprec/mld_daggrmat_minnrg_asb.F90
 mlprec/mld_daggrmat_smth_asb.F90
 mlprec/mld_dmlprec_aply.f90
 mlprec/mld_dmlprec_bld.f90
 mlprec/mld_dprecaply.f90
 mlprec/mld_saggrmat_smth_asb.F90
 mlprec/mld_zaggrmat_smth_asb.F90
 tests/pdegen/ppde.f90
 tests/pdegen/runs/ppde.inp

D: work on VECT/MAT version. Updated test program.
Initial version of MINNRGY, to be tested.
stopcriterion
Salvatore Filippone 13 years ago
parent 859ba30c4c
commit 08f013aac9
14 changed files with 1718 additions and 866 deletions
Show Stats Download Patch File Download Diff File

74
mlprec/mld_caggrmat_smth_asb.F90
Unescape Escape View File

@ -111,11 +111,11 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
#endif
! Arguments
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(psb_cspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_conelev_type), intent(inout), target :: p
integer, intent(out) :: info
integer, intent(out) :: info
! Local variables
type(psb_cspmat_type) :: b
@ -125,8 +125,8 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_cspmat_type) :: am1,am2, am3, am4
type(psb_c_coo_sparse_mat) :: acoo1, acoo2, acoof, acoo3,acoo4, bcoo, cootmp
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsrf, acsr3,acsr4, bcsr
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(:)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
@ -205,25 +205,25 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo4%allocate(ncol,ntaggr,ncol)
call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
acoo4%val(i) = cone
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = cone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(ncol)
call acoo%set_nzeros(ncol)
else
call acoo4%allocate(ncol,naggr,ncol)
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo4%val(i) = cone
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = cone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(nrow)
call acoo%set_nzeros(nrow)
endif
call acoo4%set_dupl(psb_dupl_add_)
call acoo%set_dupl(psb_dupl_add_)
call acsr4%mv_from_coo(acoo4,info)
call ptilde%mv_from_coo(acoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) &
@ -359,13 +359,13 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,acsr4,acsr1,info)
call psb_symbmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsrf,acsr4,acsr1)
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -395,20 +395,20 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,acsr4,acsr1,info)
call psb_symbmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsr3,acsr4,acsr1)
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
end if
call acsr4%free()
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1)
@ -442,24 +442,24 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_kind == mld_smooth_prol_) then
call am2%transp(am1)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of B that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < acoo2%ia(k)) .and.(acoo2%ia(k) <= naggrp1)) then
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then
i = i+1
acoo2%val(i) = acoo2%val(k)
acoo2%ia(i) = acoo2%ia(k)
acoo2%ja(i) = acoo2%ja(k)
acoo%val(i) = acoo%val(k)
acoo%ia(i) = acoo%ia(k)
acoo%ja(i) = acoo%ja(k)
end if
end do
call acoo2%set_nzeros(i)
call acoo2%trim()
call am2%mv_from(acoo2)
call acoo%set_nzeros(i)
call acoo%trim()
call am2%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2')
@ -551,11 +551,11 @@ subroutine mld_caggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo2%ia(1:nzl),p%desc_ac,info,'I')
call acoo2%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local')

666
mlprec/mld_d_as_smoother.f90
Unescape Escape View File

@ -166,6 +166,8 @@ contains
integer :: n_row,n_col, nrow_d, i
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_dpk_), allocatable :: vx(:)
type(psb_d_vect_type) :: vtx, vty, vww
integer :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_
character(len=20) :: name='d_as_smoother_apply', ch_err
@ -193,10 +195,11 @@ contains
end if
n_row = sm%desc_data%get_local_rows()
n_col = sm%desc_data%get_local_cols()
n_row = sm%desc_data%get_local_rows()
n_col = sm%desc_data%get_local_cols()
nrow_d = desc_data%get_local_rows()
isz=max(n_row,N_COL)
isz = max(n_row,N_COL)
if ((6*isz) <= size(work)) then
ww => work(1:isz)
tx => work(isz+1:2*isz)
@ -246,322 +249,344 @@ contains
else
!!$
!!$ tx(1:nrow_d) = x(1:nrow_d)
!!$ tx(nrow_d+1:isz) = dzero
!!$
!!$
!!$ if (sweeps == 1) then
!!$
!!$ select case(trans_)
!!$ case('N')
!!$ !
!!$ ! Get the overlap entries of tx (tx == x)
!!$ !
!!$ if (sm%restr == psb_halo_) then
!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_halo'
!!$ goto 9999
!!$ end if
!!$ else if (sm%restr /= psb_none_) then
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_')
!!$ goto 9999
!!$ end if
!!$
!!$
!!$ case('T','C')
!!$ !
!!$ ! With transpose, we have to do it here
!!$ !
!!$
!!$ select case (sm%prol)
!!$
!!$ case(psb_none_)
!!$ !
!!$ ! Do nothing
!!$
!!$ case(psb_sum_)
!!$ !
!!$ ! The transpose of sum is halo
!!$ !
!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_halo'
!!$ goto 9999
!!$ end if
!!$
!!$ case(psb_avg_)
!!$ !
!!$ ! Tricky one: first we have to scale the overlap entries,
!!$ ! which we can do by assignind mode=0, i.e. no communication
!!$ ! (hence only scaling), then we do the halo
!!$ !
!!$ call psb_ovrl(tx,sm%desc_data,info,&
!!$ & update=psb_avg_,work=aux,mode=0)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_ovrl'
!!$ goto 9999
!!$ end if
!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_halo'
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_')
!!$ goto 9999
!!$ end select
!!$
!!$
!!$ case default
!!$ info=psb_err_iarg_invalid_i_
!!$ int_err(1)=6
!!$ ch_err(2:2)=trans
!!$ goto 9999
!!$ end select
!!$
!!$ call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_internal_error_,name,&
!!$ & a_err='Error in sub_aply Jacobi Sweeps = 1')
!!$ goto 9999
!!$ endif
!!$
!!$ select case(trans_)
!!$ case('N')
!!$
!!$ select case (sm%prol)
!!$
!!$ case(psb_none_)
!!$ !
!!$ ! Would work anyway, but since it is supposed to do nothing ...
!!$ ! call psb_ovrl(ty,sm%desc_data,info,&
!!$ ! & update=sm%prol,work=aux)
!!$
!!$
!!$ case(psb_sum_,psb_avg_)
!!$ !
!!$ ! Update the overlap of ty
!!$ !
!!$ call psb_ovrl(ty,sm%desc_data,info,&
!!$ & update=sm%prol,work=aux)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_ovrl'
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_')
!!$ goto 9999
!!$ end select
!!$
!!$ case('T','C')
!!$ !
!!$ ! With transpose, we have to do it here
!!$ !
!!$ if (sm%restr == psb_halo_) then
!!$ call psb_ovrl(ty,sm%desc_data,info,&
!!$ & update=psb_sum_,work=aux)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_ovrl'
!!$ goto 9999
!!$ end if
!!$ else if (sm%restr /= psb_none_) then
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_')
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ info=psb_err_iarg_invalid_i_
!!$ int_err(1)=6
!!$ ch_err(2:2)=trans
!!$ goto 9999
!!$ end select
!!$
!!$
!!$
!!$ else if (sweeps > 1) then
!!$
!!$ !
!!$ !
!!$ ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver
!!$ ! to compute an approximate solution of a linear system.
!!$ !
!!$ !
!!$ ty = dzero
!!$ do i=1, sweeps
!!$ select case(trans_)
!!$ case('N')
!!$ !
!!$ ! Get the overlap entries of tx (tx == x)
!!$ !
!!$ if (sm%restr == psb_halo_) then
!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_halo'
!!$ goto 9999
!!$ end if
!!$ else if (sm%restr /= psb_none_) then
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_')
!!$ goto 9999
!!$ end if
!!$
!!$
!!$ case('T','C')
!!$ !
!!$ ! With transpose, we have to do it here
!!$ !
!!$
!!$ select case (sm%prol)
!!$
!!$ case(psb_none_)
!!$ !
!!$ ! Do nothing
!!$
!!$ case(psb_sum_)
!!$ !
!!$ ! The transpose of sum is halo
!!$ !
!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_halo'
!!$ goto 9999
!!$ end if
!!$
!!$ case(psb_avg_)
!!$ !
!!$ ! Tricky one: first we have to scale the overlap entries,
!!$ ! which we can do by assignind mode=0, i.e. no communication
!!$ ! (hence only scaling), then we do the halo
!!$ !
!!$ call psb_ovrl(tx,sm%desc_data,info,&
!!$ & update=psb_avg_,work=aux,mode=0)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_ovrl'
!!$ goto 9999
!!$ end if
!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_halo'
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_')
!!$ goto 9999
!!$ end select
!!$
!!$
!!$ case default
!!$ info=psb_err_iarg_invalid_i_
!!$ int_err(1)=6
!!$ ch_err(2:2)=trans
!!$ goto 9999
!!$ end select
!!$ !
!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
!!$ ! block diagonal part and the remaining part of the local matrix
!!$ ! and Y(j) is the approximate solution at sweep j.
!!$ !
!!$ ww(1:n_row) = tx(1:n_row)
!!$ call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,work=aux,trans=trans_)
!!$
!!$ if (info /= psb_success_) exit
!!$
!!$ call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info)
!!$
!!$ if (info /= psb_success_) exit
!!$
!!$
!!$ select case(trans_)
!!$ case('N')
!!$
!!$ select case (sm%prol)
!!$
!!$ case(psb_none_)
!!$ !
!!$ ! Would work anyway, but since it is supposed to do nothing ...
!!$ ! call psb_ovrl(ty,sm%desc_data,info,&
!!$ ! & update=sm%prol,work=aux)
!!$
!!$
!!$ case(psb_sum_,psb_avg_)
!!$ !
!!$ ! Update the overlap of ty
!!$ !
!!$ call psb_ovrl(ty,sm%desc_data,info,&
!!$ & update=sm%prol,work=aux)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_ovrl'
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_')
!!$ goto 9999
!!$ end select
!!$
!!$ case('T','C')
!!$ !
!!$ ! With transpose, we have to do it here
!!$ !
!!$ if (sm%restr == psb_halo_) then
!!$ call psb_ovrl(ty,sm%desc_data,info,&
!!$ & update=psb_sum_,work=aux)
!!$ if(info /= psb_success_) then
!!$ info=psb_err_from_subroutine_
!!$ ch_err='psb_ovrl'
!!$ goto 9999
!!$ end if
!!$ else if (sm%restr /= psb_none_) then
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_')
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ info=psb_err_iarg_invalid_i_
!!$ int_err(1)=6
!!$ ch_err(2:2)=trans
!!$ goto 9999
!!$ end select
!!$ end do
!!$
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
!!$
!!$
!!$ else
!!$
!!$ info = psb_err_iarg_neg_
!!$ call psb_errpush(info,name,&
!!$ & i_err=(/2,sweeps,0,0,0/))
!!$ goto 9999
!!$
!!$
!!$ end if
!!$
!!$ !
!!$ ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!!$ !
!!$ call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
!!$
vx = x%getCopy()
call psb_geall(vtx,sm%desc_data,info)
call psb_geasb(vtx,sm%desc_data,info,mold=x%v)
call psb_geall(vty,sm%desc_data,info)
call psb_geasb(vty,sm%desc_data,info,mold=x%v)
call psb_geall(vww,sm%desc_data,info)
call psb_geasb(vww,sm%desc_data,info,mold=x%v)
call vtx%set(dzero)
call vty%set(dzero)
call vww%set(dzero)
call vtx%set(vx(1:nrow_d))
if (sweeps == 1) then
select case(trans_)
case('N')
!
! Get the overlap entries of tx (tx == x)
!
if (sm%restr == psb_halo_) then
call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_halo'
goto 9999
end if
else if (sm%restr /= psb_none_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_restr_')
goto 9999
end if
case('T','C')
!
! With transpose, we have to do it here
!
select case (sm%prol)
case(psb_none_)
!
! Do nothing
case(psb_sum_)
!
! The transpose of sum is halo
!
call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_halo'
goto 9999
end if
case(psb_avg_)
!
! Tricky one: first we have to scale the overlap entries,
! which we can do by assignind mode=0, i.e. no communication
! (hence only scaling), then we do the halo
!
call psb_ovrl(vtx,sm%desc_data,info,&
& update=psb_avg_,work=aux,mode=0)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_halo'
goto 9999
end if
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_prol_')
goto 9999
end select
case default
info=psb_err_iarg_invalid_i_
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
call sm%sv%apply(done,vtx,dzero,vty,sm%desc_data,trans_,aux,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999
endif
select case(trans_)
case('N')
select case (sm%prol)
case(psb_none_)
!
! Would work anyway, but since it is supposed to do nothing ...
! call psb_ovrl(ty,sm%desc_data,info,&
! & update=sm%prol,work=aux)
case(psb_sum_,psb_avg_)
!
! Update the overlap of ty
!
call psb_ovrl(vty,sm%desc_data,info,&
& update=sm%prol,work=aux)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_prol_')
goto 9999
end select
case('T','C')
!
! With transpose, we have to do it here
!
if (sm%restr == psb_halo_) then
call psb_ovrl(vty,sm%desc_data,info,&
& update=psb_sum_,work=aux)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
else if (sm%restr /= psb_none_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_restr_')
goto 9999
end if
case default
info=psb_err_iarg_invalid_i_
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
else if (sweeps > 1) then
!
!
! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system.
!
!
call vty%set(dzero)
do i=1, sweeps
select case(trans_)
case('N')
!
! Get the overlap entries of tx (tx == x)
!
if (sm%restr == psb_halo_) then
call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_halo'
goto 9999
end if
else if (sm%restr /= psb_none_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_restr_')
goto 9999
end if
case('T','C')
!
! With transpose, we have to do it here
!
select case (sm%prol)
case(psb_none_)
!
! Do nothing
case(psb_sum_)
!
! The transpose of sum is halo
!
call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_halo'
goto 9999
end if
case(psb_avg_)
!
! Tricky one: first we have to scale the overlap entries,
! which we can do by assignind mode=0, i.e. no communication
! (hence only scaling), then we do the halo
!
call psb_ovrl(vtx,sm%desc_data,info,&
& update=psb_avg_,work=aux,mode=0)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(vtx,sm%desc_data,info,work=aux,data=psb_comm_ext_)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_halo'
goto 9999
end if
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_prol_')
goto 9999
end select
case default
info=psb_err_iarg_invalid_i_
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
!
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! block diagonal part and the remaining part of the local matrix
! and Y(j) is the approximate solution at sweep j.
!
call psb_geaxpby(done,vtx,dzero,vww,sm%desc_data,info)
call psb_spmm(-done,sm%nd,vty,done,vww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(done,vww,dzero,vty,sm%desc_data,trans_,aux,info)
if (info /= psb_success_) exit
select case(trans_)
case('N')
select case (sm%prol)
case(psb_none_)
!
! Would work anyway, but since it is supposed to do nothing ...
! call psb_ovrl(ty,sm%desc_data,info,&
! & update=sm%prol,work=aux)
case(psb_sum_,psb_avg_)
!
! Update the overlap of ty
!
call psb_ovrl(vty,sm%desc_data,info,&
& update=sm%prol,work=aux)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_prol_')
goto 9999
end select
case('T','C')
!
! With transpose, we have to do it here
!
if (sm%restr == psb_halo_) then
call psb_ovrl(vty,sm%desc_data,info,&
& update=psb_sum_,work=aux)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
else if (sm%restr /= psb_none_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid mld_sub_restr_')
goto 9999
end if
case default
info=psb_err_iarg_invalid_i_
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
else
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/2,sweeps,0,0,0/))
goto 9999
end if
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,vty,beta,y,desc_data,info)
end if
@ -573,6 +598,9 @@ contains
else
deallocate(ww,aux,tx,ty)
endif
call vww%free(info)
call vtx%free(info)
call vty%free(info)
call psb_erractionrestore(err_act)
return

4
mlprec/mld_d_ilu_solver.f90
Unescape Escape View File

@ -243,6 +243,8 @@ contains
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
! For real this is the same of T but in the future we will move to
! a preprocessed version, so we need a placeholder
call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,&
& trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,&
@ -354,7 +356,7 @@ contains
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
! For real this is the same of T but in the future we will move to
! a preprocessed version
! a preprocessed version, so we need a placeholder
call psb_spsm(done,sv%u,x,dzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=sv%d,choice=psb_none_,work=aux)
if (info == psb_success_) call psb_spsm(alpha,sv%l,ww,beta,y,desc_data,info,&

13
mlprec/mld_d_inner_mod.f90
Unescape Escape View File

@ -79,6 +79,19 @@ module mld_d_inner_mod
real(psb_dpk_),target :: work(:)
integer, intent(out) :: info
end subroutine mld_dmlprec_aply
subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, &
& psb_dpk_, psb_d_vect_type
use mld_d_prec_type, only : mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(inout) :: p
real(psb_dpk_),intent(in) :: alpha,beta
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
character,intent(in) :: trans
real(psb_dpk_),target :: work(:)
integer, intent(out) :: info
end subroutine mld_dmlprec_aply_vect
end interface mld_mlprec_aply

2
mlprec/mld_d_jac_smoother.f90
Unescape Escape View File

@ -170,7 +170,9 @@ contains
goto 9999
end if
call tx%bld(x%get_nrows())
call tx%set(dzero)
call ty%bld(x%get_nrows())
call ty%set(dzero)
do i=1, sweeps
!

958
mlprec/mld_daggrmat_minnrg_asb.F90
Unescape Escape View File

@ -111,7 +111,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
#endif
! Arguments
type(psb_dspmat_type), intent(in) :: a
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer, intent(inout) :: ilaggr(:), nlaggr(:)
type(mld_donelev_type), intent(inout), target :: p
@ -123,9 +123,12 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
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
character(len=20) :: name
type(psb_dspmat_type) :: am1,am2, am3, am4, atmp, atmp2, atran
type(psb_d_coo_sparse_mat) :: acoo, acoof, bcoo, tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, bcsr, acsr, acsrf, ptilde, rtilde
type(psb_d_csr_sparse_mat) :: ra, rada, arp, ardap, artp, artdatp, acrtran
type(psb_d_csc_sparse_mat) :: ap, adap, atp, atdatp, acsc
real(psb_dpk_), allocatable :: adiag(:), pj(:), xj(:), yj(:), omf(:),omp(:),omi(:),&
& oden(:), adinv(:)
logical :: filter_mat
@ -145,299 +148,537 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
!!$
!!$
!!$ call psb_nullify_sp(b)
!!$ call psb_nullify_sp(am3)
!!$ call psb_nullify_sp(am4)
!!$ call psb_nullify_sp(am1)
!!$ call psb_nullify_sp(am2)
!!$ call psb_nullify_sp(Ap)
!!$ call psb_nullify_sp(Adap)
!!$ call psb_nullify_sp(Atmp)
!!$ call psb_nullify_sp(Atmp2)
!!$ call psb_nullify_sp(Atran)
!!$ call psb_nullify_sp(Atp)
!!$ call psb_nullify_sp(atdatp)
!!$ call psb_nullify_sp(AF)
!!$ call psb_nullify_sp(ra)
!!$ call psb_nullify_sp(rada)
!!$ call psb_nullify_sp(ptilde)
!!$ call psb_nullify_sp(rtilde)
!!$
!!$ nglob = desc_a%get_global_rows()
!!$ nrow = desc_a%get_local_rows()
!!$ ncol = desc_a%get_local_cols()
!!$
!!$ theta = p%rprcparm(mld_aggr_thresh_)
!!$
!!$ naggr = nlaggr(me+1)
!!$ 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')
!!$ goto 9999
!!$ end if
!!$
!!$ naggrm1 = sum(nlaggr(1:me))
!!$ naggrp1 = sum(nlaggr(1:me+1))
!!$
!!$ filter_mat = (p%parms%aggr_filter == mld_filter_mat_)
!!$
!!$ ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
!!$ call psb_halo(ilaggr,desc_a,info)
!!$
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
!!$ goto 9999
!!$ end if
!!$
!!$ ! naggr: number of local aggregates
!!$ ! nrow: local rows.
!!$ !
!!$ allocate(adiag(ncol),adinv(ncol),xj(ncol),&
!!$ & yj(ncol),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_)')
!!$ goto 9999
!!$ end if
!!$
!!$ ! Get the diagonal D
!!$ call psb_sp_getdiag(a,adiag,info)
!!$ if (info == psb_success_) &
!!$ & call psb_halo(adiag,desc_a,info)
!!$
!!$ 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
!!$ ptilde%fida='COO'
!!$ ptilde%m=ncol
!!$ ptilde%k=ntaggr
!!$ call psb_sp_all(ncol,ntaggr,ptilde,ncol,info)
!!$
!!$
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spall')
!!$ goto 9999
!!$ end if
!!$
!!$ do i=1,ncol
!!$ ptilde%aspk(i) = done
!!$ ptilde%ia1(i) = i
!!$ ptilde%ia2(i) = ilaggr(i)
!!$ end do
!!$ ptilde%infoa(psb_nnz_) = ncol
!!$
!!$ call psb_spcnv(ptilde,info,afmt='csr',dupl=psb_dupl_add_)
!!$ if (info == psb_success_) call psb_spcnv(a,am3,info,afmt='csr',dupl=psb_dupl_add_)
!!$ if (info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
!!$ goto 9999
!!$ end if
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),&
!!$ & ' Initial copies done.'
!!$
!!$ call psb_symbmm(am3,ptilde,ap,info)
!!$ if (info == psb_success_) call psb_numbmm(am3,ptilde,ap)
!!$
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_sp_clone(ap,atmp,info)
!!$
!!$
!!$ do i=1,size(adiag)
!!$ if (adiag(i) /= dzero) then
!!$ adinv(i) = done / adiag(i)
!!$ else
!!$ adinv(i) = done
!!$ end if
!!$ end do
!!$ call psb_sp_scal(adinv,atmp,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 psb_sp_free(am4,info)
!!$
!!$ call psb_symbmm(am3,atmp,adap,info)
!!$ call psb_numbmm(am3,atmp,adap)
!!$ call psb_sp_free(atmp,info)
!!$
!!$! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
!!$! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
!!$ call psb_spcnv(ap,info,afmt='coo')
!!$ if (info == psb_success_) call psb_spcnv(ap,info,afmt='csc')
!!$ if (info == psb_success_) call psb_spcnv(adap,info,afmt='coo')
!!$ if (info == psb_success_) call psb_spcnv(adap,info,afmt='csc')
!!$ if (info /= psb_success_) then
!!$ write(0,*) 'Failed conversion to CSC'
!!$ end if
!!$
!!$ call csc_mat_col_prod(ap,adap,omp,info)
!!$ call csc_mat_col_prod(adap,adap,oden,info)
!!$ call psb_sum(ictxt,omp)
!!$ call psb_sum(ictxt,oden)
!!$! !$ write(debug_unit,*) trim(name),' OMP :',omp
!!$! !$ write(debug_unit,*) trim(name),' ODEN:',oden
!!$ omp = omp/oden
!!$! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),&
!!$ & 'Done NUMBMM 1'
!!$
!!$ ! Compute omega_int
!!$ ommx = -1d300
!!$ do i=1, ncol
!!$ omi(i) = omp(ilaggr(i))
!!$ ommx = max(ommx,omi(i))
!!$ end do
!!$ ! Compute omega_fine
!!$ do i=1, nrow
!!$ omf(i) = ommx
!!$ do j=am3%ia2(i),am3%ia2(i+1)-1
!!$ omf(i) = min(omf(i),omi(am3%ia1(j)))
!!$ end do
!!$ omf(i) = max(dzero,omf(i))
!!$ end do
!!$
!!$
!!$ if (filter_mat) then
!!$ !
!!$ ! Build the filtered matrix Af from A
!!$ !
!!$ call psb_spcnv(a,af,info,afmt='csr',dupl=psb_dupl_add_)
!!$
!!$ do i=1,nrow
!!$ tmp = dzero
!!$ jd = -1
!!$ do j=af%ia2(i),af%ia2(i+1)-1
!!$ if (af%ia1(j) == i) jd = j
!!$ if (abs(af%aspk(j)) < theta*sqrt(abs(adiag(i)*adiag(af%ia1(j))))) then
!!$ tmp=tmp+af%aspk(j)
!!$ af%aspk(j)=dzero
!!$ endif
!!$ enddo
!!$ if (jd == -1) then
!!$ write(0,*) 'Wrong input: we need the diagonal!!!!', i
!!$ else
!!$ af%aspk(jd)=af%aspk(jd)-tmp
!!$ end if
!!$ enddo
!!$ ! Take out zeroed terms
!!$ call psb_spcnv(af,info,afmt='coo')
!!$ k = 0
!!$ do j=1,psb_sp_get_nnzeros(af)
!!$ if ((af%aspk(j) /= dzero) .or. (af%ia1(j) == af%ia2(j))) then
!!$ k = k + 1
!!$ af%aspk(k) = af%aspk(j)
!!$ af%ia1(k) = af%ia1(j)
!!$ af%ia2(k) = af%ia2(j)
!!$ end if
!!$ end do
!!$! !$ write(debug_unit,*) me,' ',trim(name),' Non zeros from filtered matrix:',k,af%m,af%k
!!$ call psb_sp_setifld(k,psb_nnz_,af,info)
!!$ call psb_spcnv(af,info,afmt='csr')
!!$ end if
!!$
!!$ omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
!!$
!!$ if (filter_mat) then
!!$ !
!!$ ! Build the smoothed prolongator using the filtered matrix
!!$ !
!!$ if (psb_toupper(af%fida) == 'CSR') then
!!$ do i=1,af%m
!!$ do j=af%ia2(i),af%ia2(i+1)-1
!!$ if (af%ia1(j) == i) then
!!$ af%aspk(j) = done - omf(i)*af%aspk(j)
!!$ else
!!$ af%aspk(j) = - omf(i)*af%aspk(j)
!!$ end if
!!$ end do
!!$ end do
!!$ else
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AF storage format')
!!$ goto 9999
!!$ end if
!!$
!!$ 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.
!!$ !
!!$ ! am1 = (I-w*D*Af)Ptilde
!!$ ! Doing it this way means to consider diag(Af_i)
!!$ !
!!$ !
!!$ call psb_symbmm(af,ptilde,am1,info)
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_numbmm(af,ptilde,am1)
!!$
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),&
!!$ & 'Done NUMBMM 1'
!!$ else
!!$ !
!!$ ! Build the smoothed prolongator using the original matrix
!!$ !
!!$ if (psb_toupper(am3%fida) == 'CSR') then
!!$ do i=1,am3%m
!!$ do j=am3%ia2(i),am3%ia2(i+1)-1
!!$ if (am3%ia1(j) == i) then
!!$ am3%aspk(j) = done - omf(i)*am3%aspk(j)
!!$ else
!!$ am3%aspk(j) = - omf(i)*am3%aspk(j)
!!$ end if
!!$ end do
!!$ end do
!!$ else
!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid AM3 storage format')
!!$ goto 9999
!!$ end if
!!$
!!$ 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.
!!$ !
!!$ ! am1 = (I-w*D*A)Ptilde
!!$ !
!!$ !
!!$ call psb_symbmm(am3,ptilde,am1,info)
!!$ if(info /= psb_success_) then
!!$ call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
!!$ goto 9999
!!$ end if
!!$
!!$ call psb_numbmm(am3,ptilde,am1)
!!$
!!$ if (debug_level >= psb_debug_outer_) &
!!$ & write(debug_unit,*) me,' ',trim(name),&
!!$ & 'Done NUMBMM 1'
!!$
!!$ end if
!!$
!!$ !
!!$ ! Ok, let's start over with the restrictor
!!$ !
nglob = desc_a%get_global_rows()
nrow = desc_a%get_local_rows()
ncol = desc_a%get_local_cols()
theta = p%parms%aggr_thresh
naggr = nlaggr(me+1)
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')
goto 9999
end if
naggrm1 = sum(nlaggr(1:me))
naggrp1 = sum(nlaggr(1:me+1))
filter_mat = (p%parms%aggr_filter == mld_filter_mat_)
ilaggr(1:nrow) = ilaggr(1:nrow) + naggrm1
call psb_halo(ilaggr,desc_a,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='psb_halo')
goto 9999
end if
! naggr: number of local aggregates
! nrow: local rows.
!
allocate(adiag(ncol),adinv(ncol),xj(ncol),&
& yj(ncol),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_)')
goto 9999
end if
! Get the diagonal D
! Get the diagonal D
call a%get_diag(adiag,info)
if (info == psb_success_) &
& call psb_halo(adiag,desc_a,info)
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 acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
acoo%val(i) = done
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo%set_nzeros(ncol)
call acoo%set_dupl(psb_dupl_add_)
call ptilde%mv_from_coo(acoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Initial copies done.'
call psb_symbmm(acsr3,ptilde,arp,info)
if (info == psb_success_) call psb_numbmm(acsr3,ptilde,arp)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call atmp%cp_from(arp)
do i=1,size(adiag)
if (adiag(i) /= dzero) then
adinv(i) = done / adiag(i)
else
adinv(i) = done
end if
end do
call atmp%scal(adinv,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 psb_sp_free(am4,info)
call atmp%mv_to(acsr1)
call psb_symbmm(acsr3,acsr1,ardap,info)
call psb_numbmm(acsr3,acsr1,ardap)
call acsr1%free()
! !$ write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
! !$ write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
call ap%mv_from_fmt(arp,info)
call adap%mv_from_fmt(ardap,info)
call csc_mat_col_prod(ap,adap,omp,info)
call csc_mat_col_prod(adap,adap,oden,info)
call psb_sum(ictxt,omp)
call psb_sum(ictxt,oden)
! !$ write(0,*) trim(name),' OMP :',omp
! !$ write(0,*) trim(name),' ODEN:',oden
omp = omp/oden
! !$ write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
! Compute omega_int
ommx = -1d300
do i=1, ncol
omi(i) = omp(ilaggr(i))
ommx = max(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)))
end do
omf(i) = max(dzero,omf(i))
end do
omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
if (filter_mat) then
!
! Build the filtered matrix Af from A
!
call a%cscnv(acsrf,info,dupl=psb_dupl_add_)
do i=1,nrow
tmp = dzero
jd = -1
do j=acsrf%irp(i),acsrf%irp(i+1)-1
if (acsrf%ja(j) == i) jd = j
if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
tmp=tmp+acsrf%val(j)
acsrf%val(j)=dzero
endif
enddo
if (jd == -1) then
write(0,*) 'Wrong input: we need the diagonal!!!!', i
else
acsrf%val(jd)=acsrf%val(jd)-tmp
end if
enddo
! Take out zeroed terms
call acsrf%mv_to_coo(tmpcoo,info)
k = 0
do j=1,tmpcoo%get_nzeros()
if ((tmpcoo%val(j) /= dzero) .or. (tmpcoo%ia(j) == tmpcoo%ja(j))) then
k = k + 1
tmpcoo%val(k) = tmpcoo%val(j)
tmpcoo%ia(k) = tmpcoo%ia(j)
tmpcoo%ja(k) = tmpcoo%ja(j)
end if
end do
! !$ write(debug_unit,*) me,' ',trim(name),' Non zeros from filtered matrix:',k,acsrf%m,acsrf%k
call tmpcoo%set_nzeros(k)
call acsrf%mv_from_coo(tmpcoo,info)
!
! 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 - omf(i)*acsrf%val(j)
else
acsrf%val(j) = - omf(i)*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.
!
! am1 = (I-w*D*Af)Ptilde
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 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 - omf(i)*acsr3%val(j)
else
acsr3%val(j) = - omf(i)*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.
!
! am1 = (I-w*D*A)Ptilde
!
!
call psb_symbmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
end if
!
! Now encapsulate in AM1
!
call am1%mv_from(acsr1)
!
! Ok, let's start over with the restrictor
!
call ptilde%transp(rtilde)
call atmp%mv_from(acsr3)
call psb_sphalo(atmp,desc_a,am4,info,&
& colcnv=.true.,rowscale=.true.)
nrt = am4%get_nrows()
call am4%csclip(atmp2,info,1,nrt,1,ncol)
call atmp2%cscnv(info,type='CSR')
if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)
call am4%free()
call atmp2%free()
! This is to compute the transpose. It ONLY works if the
! original A has a symmetric pattern.
call atmp%transp(atmp2)
call atmp2%csclip(atran,info,1,nrow,1,ncol)
call atmp2%free()
call atran%mv_to(acrtran)
! Now for the product.
call psb_symbmm(acrtran,ptilde,artp,info)
if (info == psb_success_) call psb_numbmm(acrtran,ptilde,artp)
call atmp2%cp_from(artp)
call atmp2%scal(adinv,info)
call psb_sphalo(atmp2,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.,outfmt='CSR ')
if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)
if (info == psb_success_) call am4%free()
call atmp2%mv_to(acsr2)
call psb_symbmm(acrtran,acsr2,artdatp,info)
call psb_numbmm(acrtran,acsr2,artdatp)
call acsr2%free()
call artp%mv_to_fmt(atp,info)
call artdatp%mv_to_fmt(atdatp,info)
call csc_mat_col_prod(atp,atdatp,omp,info)
call csc_mat_col_prod(atdatp,atdatp,oden,info)
call psb_sum(ictxt,omp)
call psb_sum(ictxt,oden)
! !$ write(debug_unit,*) trim(name),' OMP_R :',omp
! ! $ write(debug_unit,*) trim(name),' ODEN_R:',oden
omp = omp/oden
! !$ write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
! Compute omega_int
ommx = -1d300
do i=1, ncol
omi(i) = omp(ilaggr(i))
ommx = max(ommx,omi(i))
end do
! Compute omega_fine
! Going over the columns of atmp means going over the rows
! of A^T. Hopefully ;-)
call atmp%cp_to(acsc)
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)))
end do
omf(i) = max(dzero,omf(i))
end do
omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
call psb_halo(omf,desc_a,info)
call acsc%free()
call atmp%mv_to(acsr1)
do i=1,acsr1%get_nrows()
do j=acsr1%irp(i),acsr1%irp(i+1)-1
if (acsr1%ja(j) == i) then
acsr1%val(j) = done - acsr1%val(j)*omf(acsr1%ja(j))
else
acsr1%val(j) = - acsr1%val(j)*omf(acsr1%ja(j))
end if
end do
end do
call psb_symbmm(rtilde,acsr1,acsr2,info)
call psb_numbmm(rtilde,acsr1,acsr2)
!
! Now we have to fix this. The only rows of B that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
call acsr2%mv_to_coo(tmpcoo,info)
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 acsr2%mv_from_coo(tmpcoo,info)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd'
!
! Final steps: build the product.
! Now we have to gather the halo of am1, and add it to itself
! to multiply it by A,
!
call psb_sphalo(am1,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am1,info,b=am4)
if (info == psb_success_) call am4%free()
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of am1')
goto 9999
end if
call am1%cp_to(acsr1)
call a%cp_to(acsr)
call psb_symbmm(acsr,acsr1,acsr3,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999
end if
call psb_numbmm(acsr,acsr1,acsr3)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2'
call am3%mv_from(acsr3)
! am2 = ((i-wDA)Ptilde)^T
call psb_sphalo(am3,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)
if (info == psb_success_) call am4%free()
call am3%mv_to(acsr3)
call psb_symbmm(acsr2,acsr3,bcsr,info)
if (info == psb_success_) call psb_numbmm(acsr2,acsr3,bcsr)
call acsr3%free()
call bcsr%mv_to_coo(bcoo,info)
call am2%mv_from(acsr2)
select case(p%parms%coarse_mat)
case(mld_distr_mat_)
nzl = bcoo%get_nrows()
if (info == psb_success_) call psb_cdall(ictxt,p%desc_ac,info,nl=nlaggr(me+1))
if (info == psb_success_) call psb_cdins(nzl,bcoo%ia,bcoo%ja,p%desc_ac,info)
if (info == psb_success_) call psb_cdasb(p%desc_ac,info)
if (info == psb_success_) call psb_glob_to_loc(bcoo%ia(1:nzl),p%desc_ac,info,iact='I')
if (info == psb_success_) call psb_glob_to_loc(bcoo%ja(1:nzl),p%desc_ac,info,iact='I')
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Creating p%desc_ac and converting ac')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Assembld aux descr. distr.'
call bcoo%set_nrows(p%desc_ac%get_local_rows())
call bcoo%set_ncols(p%desc_ac%get_local_cols())
call p%ac%mv_from(bcoo)
if (np>1) then
call am1%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros()
call psb_glob_to_loc(tmpcoo%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(tmpcoo)
endif
call am1%set_ncols(p%desc_ac%get_local_cols())
if (np>1) then
call am2%mv_to(tmpcoo)
nzl = tmpcoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(tmpcoo%ia(1:nzl),p%desc_ac,info,'I')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local')
goto 9999
end if
call am2%mv_from(tmpcoo)
end if
call am2%set_nrows(p%desc_ac%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done ac '
case(mld_repl_mat_)
!
!
call psb_cdall(ictxt,p%desc_ac,info,mg=ntaggr,repl=.true.)
nzbr(:) = 0
nzbr(me+1) = bcoo%get_nzeros()
call psb_sum(ictxt,nzbr(1:np))
nzac = sum(nzbr)
if (info == psb_success_) call tmpcoo%allocate(ntaggr,ntaggr,nzac)
if (info /= psb_success_) goto 9999
do ip=1,np
idisp(ip) = sum(nzbr(1:ip-1))
enddo
ndx = nzbr(me+1)
call mpi_allgatherv(bcoo%val,ndx,mpi_double_precision,tmpcoo%val,nzbr,idisp,&
& mpi_double_precision,icomm,info)
if (info == psb_success_) call mpi_allgatherv(bcoo%ia,ndx,mpi_integer,tmpcoo%ia,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info == psb_success_) call mpi_allgatherv(bcoo%ja,ndx,mpi_integer,tmpcoo%ja,nzbr,idisp,&
& mpi_integer,icomm,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
goto 9999
end if
call tmpcoo%fix(info)
call p%ac%mv_from(tmpcoo)
call bcoo%free()
if(info /= psb_success_) goto 9999
deallocate(nzbr,idisp,stat=info)
if (info /= psb_success_) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid mld_coarse_mat_')
goto 9999
end select
call p%ac%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
goto 9999
end if
!
! Copy the prolongation/restriction matrices into the descriptor map.
! am2 => R i.e. restriction operator
! am1 => P i.e. prolongation operator
!
p%map = psb_linmap(psb_map_aggr_,desc_a,&
& p%desc_ac,am2,am1,ilaggr,nlaggr)
if (info == psb_success_) call am1%free()
if (info == psb_success_) call am2%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_Free')
goto 9999
end if
!!$ if (.false.) then
!!$ i = 4
!!$ select case (i)
@ -850,76 +1091,59 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,p,info)
contains
!!$ subroutine csc_mat_col_prod(a,b,v,info)
!!$ type(psb_dspmat_type), intent(in) :: a, b
!!$ real(psb_dpk_), intent(out) :: v(:)
!!$ integer, intent(out) :: info
!!$
!!$ integer :: i,j,k, nr, nc,iap,nra,ibp,nrb
!!$ logical :: csca, cscb
!!$
!!$ info = psb_success_
!!$ nc = psb_sp_get_ncols(a)
!!$ if (nc /= psb_sp_get_ncols(b)) then
!!$ write(0,*) 'Matrices A and B should have same columns'
!!$ info = -1
!!$ return
!!$ end if
!!$ csca = (psb_toupper(a%fida(1:3)) == 'CSC')
!!$ cscb = (psb_toupper(b%fida(1:3)) == 'CSC')
!!$
!!$ if (.not.(csca.and.cscb)) then
!!$ write(0,*) 'Matrices A and B should be in CSC'
!!$ info = -2
!!$ return
!!$ end if
!!$
!!$ do j=1, nc
!!$ iap = a%ia2(j)
!!$ nra = a%ia2(j+1)-iap
!!$ ibp = b%ia2(j)
!!$ nrb = b%ia2(j+1)-ibp
!!$ v(j) = sparse_srtd_dot(nra,a%ia1(iap:iap+nra-1),a%aspk(iap:iap+nra-1),&
!!$ & nrb,b%ia1(ibp:ibp+nrb-1),b%aspk(ibp:ibp+nrb-1))
!!$ end do
!!$
!!$ end subroutine csc_mat_col_prod
!!$
!!$
!!$ subroutine csr_mat_row_prod(a,b,v,info)
!!$ type(psb_dspmat_type), intent(in) :: a, b
!!$ real(psb_dpk_), intent(out) :: v(:)
!!$ integer, intent(out) :: info
!!$
!!$ integer :: i,j,k, nr, nc,iap,nca,ibp,ncb
!!$ logical :: csra, csrb
!!$
!!$ info = psb_success_
!!$ nr = psb_sp_get_nrows(a)
!!$ if (nr /= psb_sp_get_nrows(b)) then
!!$ write(0,*) 'Matrices A and B should have same rows'
!!$ info = -1
!!$ return
!!$ end if
!!$ csra = (psb_toupper(a%fida(1:3)) == 'CSR')
!!$ csrb = (psb_toupper(b%fida(1:3)) == 'CSR')
!!$
!!$ if (.not.(csra.and.csrb)) then
!!$ write(0,*) 'Matrices A and B should be in CSR'
!!$ info = -2
!!$ return
!!$ end if
!!$
!!$ do j=1, nr
!!$ iap = a%ia2(j)
!!$ nca = a%ia2(j+1)-iap
!!$ ibp = b%ia2(j)
!!$ ncb = b%ia2(j+1)-ibp
!!$ v(j) = sparse_srtd_dot(nca,a%ia1(iap:iap+nca-1),a%aspk(iap:iap+nca-1),&
!!$ & ncb,b%ia1(ibp:ibp+ncb-1),b%aspk(ibp:ibp+ncb-1))
!!$ end do
!!$
!!$ end subroutine csr_mat_row_prod
subroutine csc_mat_col_prod(a,b,v,info)
type(psb_d_csc_sparse_mat), intent(in) :: a, b
real(psb_dpk_), intent(out) :: v(:)
integer, intent(out) :: info
integer :: i,j,k, nr, nc,iap,nra,ibp,nrb
info = psb_success_
nc = a%get_ncols()
if (nc /= b%get_ncols()) then
write(0,*) 'Matrices A and B should have same columns'
info = -1
return
end if
do j=1, nc
iap = a%icp(j)
nra = a%icp(j+1)-iap
ibp = b%icp(j)
nrb = b%icp(j+1)-ibp
v(j) = sparse_srtd_dot(nra,a%ia(iap:iap+nra-1),a%val(iap:iap+nra-1),&
& nrb,b%ia(ibp:ibp+nrb-1),b%val(ibp:ibp+nrb-1))
end do
end subroutine csc_mat_col_prod
subroutine csr_mat_row_prod(a,b,v,info)
type(psb_d_csr_sparse_mat), intent(in) :: a, b
real(psb_dpk_), intent(out) :: v(:)
integer, intent(out) :: info
integer :: i,j,k, nr, nc,iap,nca,ibp,ncb
info = psb_success_
nr = a%get_nrows()
if (nr /= b%get_nrows()) then
write(0,*) 'Matrices A and B should have same rows'
info = -1
return
end if
do j=1, nr
iap = a%irp(j)
nca = a%irp(j+1)-iap
ibp = b%irp(j)
ncb = b%irp(j+1)-ibp
v(j) = sparse_srtd_dot(nca,a%ja(iap:iap+nca-1),a%val(iap:iap+nca-1),&
& ncb,b%ja(ibp:ibp+ncb-1),b%val(ibp:ibp+ncb-1))
end do
end subroutine csr_mat_row_prod
function sparse_srtd_dot(nv1,iv1,v1,nv2,iv2,v2) result(dot)
integer, intent(in) :: nv1,nv2

66
mlprec/mld_daggrmat_smth_asb.F90
Unescape Escape View File

@ -125,8 +125,8 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_dspmat_type) :: am1,am2, am3, am4
type(psb_d_coo_sparse_mat) :: acoo1, acoo2, acoof, acoo3,acoo4, bcoo, cootmp
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsrf, acsr3,acsr4, bcsr
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(:)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
@ -205,25 +205,25 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo4%allocate(ncol,ntaggr,ncol)
call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
acoo4%val(i) = done
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = done
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(ncol)
call acoo%set_nzeros(ncol)
else
call acoo4%allocate(ncol,naggr,ncol)
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo4%val(i) = done
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = done
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(nrow)
call acoo%set_nzeros(nrow)
endif
call acoo4%set_dupl(psb_dupl_add_)
call acoo%set_dupl(psb_dupl_add_)
call acsr4%mv_from_coo(acoo4,info)
call ptilde%mv_from_coo(acoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) &
@ -359,13 +359,13 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,acsr4,acsr1,info)
call psb_symbmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsrf,acsr4,acsr1)
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -395,20 +395,20 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,acsr4,acsr1,info)
call psb_symbmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsr3,acsr4,acsr1)
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
end if
call acsr4%free()
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1)
@ -442,24 +442,24 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_kind == mld_smooth_prol_) then
call am2%transp(am1)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of B that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < acoo2%ia(k)) .and.(acoo2%ia(k) <= naggrp1)) then
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then
i = i+1
acoo2%val(i) = acoo2%val(k)
acoo2%ia(i) = acoo2%ia(k)
acoo2%ja(i) = acoo2%ja(k)
acoo%val(i) = acoo%val(k)
acoo%ia(i) = acoo%ia(k)
acoo%ja(i) = acoo%ja(k)
end if
end do
call acoo2%set_nzeros(i)
call acoo2%trim()
call am2%mv_from(acoo2)
call acoo%set_nzeros(i)
call acoo%trim()
call am2%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2')
@ -551,11 +551,11 @@ subroutine mld_daggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo2%ia(1:nzl),p%desc_ac,info,'I')
call acoo2%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local')

573
mlprec/mld_dmlprec_aply.f90
Unescape Escape View File

@ -440,7 +440,8 @@ contains
if (level > 1) then
nc2l = psb_cd_get_local_cols(p%precv(level)%base_desc)
nr2l = psb_cd_get_local_rows(p%precv(level)%base_desc)
allocate(mlprec_wrk(level)%x2l(nc2l),mlprec_wrk(level)%y2l(nc2l),&
allocate(mlprec_wrk(level)%x2l(nc2l),&
& mlprec_wrk(level)%y2l(nc2l),&
& stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
@ -857,3 +858,573 @@ contains
end subroutine mld_dmlprec_aply
subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_dmlprec_aply_vect
implicit none
! Arguments
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(inout) :: p
real(psb_dpk_),intent(in) :: alpha,beta
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: ictxt, np, me, err_act
integer :: debug_level, debug_unit, nlev,nc2l,nr2l,level
character(len=20) :: name
character :: trans_
type psb_mlprec_wrk_type
real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l
end type psb_mlprec_wrk_type
type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:)
name='mld_dmlprec_aply'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_data%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name),&
& ' Entry ', size(p%precv)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
call psb_geall(mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,mold=x%v)
call psb_geall(mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info,mold=x%v)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call mlprec_wrk(level)%vy2l%set(dzero)
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,&
& p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
contains
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
implicit none
! Arguments
integer :: level
type(mld_dprec_type), intent(inout) :: p
type(psb_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: ictxt,np,me,i, nr2l,nc2l,err_act
integer :: debug_level, debug_unit
integer :: nlev, ilev, sweeps
character(len=20) :: name
name = 'inner_ml_aply'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_ml')
goto 9999
end if
ictxt = psb_cd_get_context(p%precv(level)%base_desc)
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if (level > 1) then
call psb_geall(mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info, &
& mold=mlprec_wrk(1)%vx2l%v)
call psb_geall(mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info, &
& mold=mlprec_wrk(1)%vx2l%v)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
end if
select case(p%precv(level)%parms%ml_type)
case(mld_no_ml_)
!
! No preconditioning, should not really get here
!
call psb_errpush(psb_err_internal_error_,name,&
& a_err='mld_no_ml_ in mlprc_aply?')
goto 9999
case(mld_add_ml_)
!
! Additive multilevel
!
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) goto 9999
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) goto 9999
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) goto 9999
end if
case(mld_mult_ml_)
!
! Multiplicative multilevel (multiplicative among the levels, additive inside
! each level)
!
! Pre/post-smoothing versions.
! Note that the transpose switches pre <-> post.
!
select case(p%precv(level)%parms%smoother_pos)
case(mld_post_smooth_)
select case (trans_)
case('N')
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
! This is one step of post-smoothing
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) goto 9999
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) goto 9999
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) goto 9999
sweeps = p%precv(level)%parms%sweeps_post
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
case('T','C')
! Post-smoothing transpose is pre-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
!
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_post
else
sweeps = p%precv(level)%parms%sweeps
end if
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) goto 9999
!
! Compute the residual (at all levels but the coarsest one)
!
if (level < nlev) then
call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) goto 9999
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) goto 9999
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) goto 9999
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans')
goto 9999
end select
case(mld_pre_smooth_)
select case (trans_)
case('N')
! One step of pre-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
!
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_pre
else
sweeps = p%precv(level)%parms%sweeps
end if
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) goto 9999
!
! Compute the residual (at all levels but the coarsest one)
!
if (level < nlev) then
call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) goto 9999
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) goto 9999
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) goto 9999
end if
case('T','C')
! pre-smooth transpose is post-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) goto 9999
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) goto 9999
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) goto 9999
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans')
goto 9999
end select
case(mld_twoside_smooth_)
call psb_geall(mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info, &
& mold=mlprec_wrk(1)%vx2l%v)
call psb_geall(mlprec_wrk(level)%vty,p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%vty,p%precv(level)%base_desc,info, &
& mold=mlprec_wrk(1)%vx2l%v)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vty,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
else
sweeps = p%precv(level)%parms%sweeps_post
end if
else
sweeps = p%precv(level)%parms%sweeps
end if
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_ ) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
else
sweeps = p%precv(level)%parms%sweeps_pre
end if
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error: residual/baseprec_aply')
goto 9999
end if
endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid smooth_pos',&
& i_Err=(/p%precv(level)%parms%smoother_pos,0,0,0,0/))
goto 9999
end select
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,0,0,0,0/))
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine inner_ml_aply
end subroutine mld_dmlprec_aply_vect

1
mlprec/mld_dmlprec_bld.f90
Unescape Escape View File

@ -314,6 +314,7 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold)
if ((info == psb_success_).and.(i>1).and.(present(amold))) then
call psb_map_cscnv(p%precv(i)%map,info,mold=amold)
call p%precv(i)%ac%cscnv(info,mold=amold)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

9
mlprec/mld_dprecaply.f90
Unescape Escape View File

@ -125,6 +125,7 @@ subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
call psb_errpush(info,name)
goto 9999
end if
if (size(prec%precv) >1) then
!
! Number of levels > 1: apply the multilevel preconditioner
@ -324,8 +325,8 @@ subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work)
!
! Number of levels > 1: apply the multilevel preconditioner
!
!!$ call mld_mlprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info)
info = psb_err_missing_override_method_
call mld_mlprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply')
goto 9999
@ -337,7 +338,9 @@ subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work)
!
call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,&
& prec%precv(1)%parms%sweeps, work_,info)
else
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Invalid size of precv',&
& i_Err=(/size(prec%precv),0,0,0,0/))
@ -345,7 +348,7 @@ subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work)
endif
! If the original distribution has an overlap we should fix that.
!!$ call psb_halo(y,desc_data,info,data=psb_comm_mov_)
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (present(work)) then

66
mlprec/mld_saggrmat_smth_asb.F90
Unescape Escape View File

@ -125,8 +125,8 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_sspmat_type) :: am1,am2, am3, am4
type(psb_s_coo_sparse_mat) :: acoo1, acoo2, acoof, acoo3,acoo4, bcoo, cootmp
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsrf, acsr3,acsr4, bcsr
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(:)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
@ -205,25 +205,25 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo4%allocate(ncol,ntaggr,ncol)
call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
acoo4%val(i) = sone
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = sone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(ncol)
call acoo%set_nzeros(ncol)
else
call acoo4%allocate(ncol,naggr,ncol)
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo4%val(i) = sone
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = sone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(nrow)
call acoo%set_nzeros(nrow)
endif
call acoo4%set_dupl(psb_dupl_add_)
call acoo%set_dupl(psb_dupl_add_)
call acsr4%mv_from_coo(acoo4,info)
call ptilde%mv_from_coo(acoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) &
@ -359,13 +359,13 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,acsr4,acsr1,info)
call psb_symbmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsrf,acsr4,acsr1)
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -395,20 +395,20 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,acsr4,acsr1,info)
call psb_symbmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsr3,acsr4,acsr1)
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
end if
call acsr4%free()
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1)
@ -442,24 +442,24 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_kind == mld_smooth_prol_) then
call am2%transp(am1)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of B that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < acoo2%ia(k)) .and.(acoo2%ia(k) <= naggrp1)) then
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then
i = i+1
acoo2%val(i) = acoo2%val(k)
acoo2%ia(i) = acoo2%ia(k)
acoo2%ja(i) = acoo2%ja(k)
acoo%val(i) = acoo%val(k)
acoo%ia(i) = acoo%ia(k)
acoo%ja(i) = acoo%ja(k)
end if
end do
call acoo2%set_nzeros(i)
call acoo2%trim()
call am2%mv_from(acoo2)
call acoo%set_nzeros(i)
call acoo%trim()
call am2%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2')
@ -551,11 +551,11 @@ subroutine mld_saggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo2%ia(1:nzl),p%desc_ac,info,'I')
call acoo2%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local')

66
mlprec/mld_zaggrmat_smth_asb.F90
Unescape Escape View File

@ -125,8 +125,8 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
integer ::ictxt,np,me, err_act, icomm
character(len=20) :: name
type(psb_zspmat_type) :: am1,am2, am3, am4
type(psb_z_coo_sparse_mat) :: acoo1, acoo2, acoof, acoo3,acoo4, bcoo, cootmp
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsrf, acsr3,acsr4, bcsr
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(:)
logical :: ml_global_nmb, filter_mat
integer :: debug_level, debug_unit
@ -205,25 +205,25 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! 1. Allocate Ptilde in sparse matrix form
if (ml_global_nmb) then
call acoo4%allocate(ncol,ntaggr,ncol)
call acoo%allocate(ncol,ntaggr,ncol)
do i=1,ncol
acoo4%val(i) = zone
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = zone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(ncol)
call acoo%set_nzeros(ncol)
else
call acoo4%allocate(ncol,naggr,ncol)
call acoo%allocate(ncol,naggr,ncol)
do i=1,nrow
acoo4%val(i) = zone
acoo4%ia(i) = i
acoo4%ja(i) = ilaggr(i)
acoo%val(i) = zone
acoo%ia(i) = i
acoo%ja(i) = ilaggr(i)
end do
call acoo4%set_nzeros(nrow)
call acoo%set_nzeros(nrow)
endif
call acoo4%set_dupl(psb_dupl_add_)
call acoo%set_dupl(psb_dupl_add_)
call acsr4%mv_from_coo(acoo4,info)
call ptilde%mv_from_coo(acoo,info)
if (info == psb_success_) call a%cscnv(acsr3,info,dupl=psb_dupl_add_)
if (debug_level >= psb_debug_outer_) &
@ -359,13 +359,13 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(Af_i)
!
!
call psb_symbmm(acsrf,acsr4,acsr1,info)
call psb_symbmm(acsrf,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsrf,acsr4,acsr1)
call psb_numbmm(acsrf,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
@ -395,20 +395,20 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
! Doing it this way means to consider diag(A_i)
!
!
call psb_symbmm(acsr3,acsr4,acsr1,info)
call psb_symbmm(acsr3,ptilde,acsr1,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 1')
goto 9999
end if
call psb_numbmm(acsr3,acsr4,acsr1)
call psb_numbmm(acsr3,ptilde,acsr1)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 1'
end if
call acsr4%free()
call ptilde%free()
call acsr1%set_dupl(psb_dupl_add_)
call am1%mv_from(acsr1)
@ -442,24 +442,24 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (p%parms%aggr_kind == mld_smooth_prol_) then
call am2%transp(am1)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
i=0
!
! Now we have to fix this. The only rows of B that are correct
! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
!
do k=1, nzl
if ((naggrm1 < acoo2%ia(k)) .and.(acoo2%ia(k) <= naggrp1)) then
if ((naggrm1 < acoo%ia(k)) .and.(acoo%ia(k) <= naggrp1)) then
i = i+1
acoo2%val(i) = acoo2%val(k)
acoo2%ia(i) = acoo2%ia(k)
acoo2%ja(i) = acoo2%ja(k)
acoo%val(i) = acoo%val(k)
acoo%ia(i) = acoo%ia(k)
acoo%ja(i) = acoo%ja(k)
end if
end do
call acoo2%set_nzeros(i)
call acoo2%trim()
call am2%mv_from(acoo2)
call acoo%set_nzeros(i)
call acoo%trim()
call am2%mv_from(acoo)
call am2%cscnv(info,type='csr',dupl=psb_dupl_add_)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv am2')
@ -554,11 +554,11 @@ subroutine mld_zaggrmat_smth_asb(a,desc_a,ilaggr,nlaggr,p,info)
if (np>1) then
call am2%cscnv(info,type='coo',dupl=psb_dupl_add_)
call am2%mv_to(acoo2)
nzl = acoo2%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo2%ia(1:nzl),p%desc_ac,info,'I')
call acoo2%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo2)
call am2%mv_to(acoo)
nzl = acoo%get_nzeros()
if (info == psb_success_) call psb_glob_to_loc(acoo%ia(1:nzl),p%desc_ac,info,'I')
call acoo%set_dupl(psb_dupl_add_)
if (info == psb_success_) call am2%mv_from(acoo)
if (info == psb_success_) call am2%cscnv(info,type='csr')
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Converting am2 to local')

76
tests/pdegen/ppde.f90
Unescape Escape View File

@ -91,7 +91,8 @@ program ppde
! descriptor
type(psb_desc_type) :: desc_a
! dense matrices
real(psb_dpk_), allocatable :: b(:), x(:)
type(psb_d_vect_type) :: x,b, vtst
real(psb_dpk_), allocatable :: tst(:)
! blacs parameters
integer :: ictxt, iam, np
@ -174,8 +175,10 @@ program ppde
end if
if (iam == psb_root_) write(*,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) write(*,'(" ")')
if (iam == psb_root_) &
& write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) &
& write(psb_out_unit,'(" ")')
!
! prepare the preconditioner.
!
@ -229,14 +232,17 @@ program ppde
call psb_amx(ictxt,tprec)
if (iam == psb_root_) write(*,'("Preconditioner time : ",es12.5)')tprec
if (iam == psb_root_) &
& write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec
if (iam == psb_root_) call mld_precdescr(prec,info)
if (iam == psb_root_) write(*,'(" ")')
if (iam == psb_root_) &
& write(psb_out_unit,'(" ")')
!
! iterative method parameters
!
if(iam == psb_root_) write(*,'("Calling iterative method ",a)')kmethd
if(iam == psb_root_) &
& write(psb_out_unit,'("Calling iterative method ",a)')kmethd
call psb_barrier(ictxt)
t1 = psb_wtime()
call psb_krylov(kmethd,a,prec,b,x,eps,desc_a,info,&
@ -254,21 +260,22 @@ program ppde
call psb_amx(ictxt,t2)
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
!!$ descsize = psb_sizeof(desc_a)
descsize = desc_a%sizeof()
precsize = mld_sizeof(prec)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
call psb_sum(ictxt,precsize)
if (iam == psb_root_) then
write(*,'(" ")')
write(*,'("Time to solve matrix : ",es12.5)')t2
write(*,'("Time per iteration : ",es12.5)')t2/iter
write(*,'("Number of iterations : ",i0)')iter
write(*,'("Convergence indicator on exit : ",es12.5)')err
write(*,'("Info on exit : ",i0)')info
write(*,'("Total memory occupation for A: ",i12)')amatsize
write(*,'("Total memory occupation for DESC_A: ",i12)')descsize
write(*,'("Total memory occupation for PREC: ",i12)')precsize
write(psb_out_unit,'(" ")')
write(psb_out_unit,'("Time to solve matrix : ",es12.5)') t2
write(psb_out_unit,'("Time per iteration : ",es12.5)') t2/iter
write(psb_out_unit,'("Number of iterations : ",i0)') iter
write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)') err
write(psb_out_unit,'("Info on exit : ",i0)') info
write(psb_out_unit,'("Total memory occupation for A: ",i12)') amatsize
write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)') descsize
write(psb_out_unit,'("Total memory occupation for PREC: ",i12)') precsize
end if
!
@ -380,13 +387,13 @@ contains
end if
if (iam == psb_root_) then
write(*,'("Solving matrix : ell1")')
write(*,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
write(*,'("Number of processors : ",i0)') np
write(*,'("Data distribution : BLOCK")')
write(*,'("Preconditioner : ",a)') prectype%descr
write(*,'("Iterative method : ",a)') kmethd
write(*,'(" ")')
write(psb_out_unit,'("Solving matrix : ell1")')
write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
write(psb_out_unit,'("Number of processors : ",i0)') np
write(psb_out_unit,'("Data distribution : BLOCK")')
write(psb_out_unit,'("Preconditioner : ",a)') prectype%descr
write(psb_out_unit,'("Iterative method : ",a)') kmethd
write(psb_out_unit,'(" ")')
endif
return
@ -438,7 +445,7 @@ contains
implicit none
integer :: idim
integer, parameter :: nb=20
real(psb_dpk_), allocatable :: b(:),xv(:)
type(psb_d_vect_type) :: b,xv
type(psb_desc_type) :: desc_a
integer :: ictxt, info
character :: afmt*5
@ -452,12 +459,12 @@ contains
real(psb_dpk_), allocatable :: val(:)
! deltah dimension of each grid cell
! deltat discretization time
real(psb_dpk_) :: deltah, deltah2
real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0
real(psb_dpk_) :: deltah, deltah2
real(psb_dpk_), parameter :: rhs=0.d0,one=1.d0,zero=0.d0
real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen
real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3
external :: a1, a2, a3, a4, b1, b2, b3
integer :: err_act
external :: a1, a2, a3, a4, b1, b2, b3
integer :: err_act
character(len=20) :: name, ch_err
@ -494,7 +501,7 @@ contains
! define rhs from boundary conditions; also build initial guess
if (info == psb_success_) call psb_geall(b,desc_a,info)
if (info == psb_success_) call psb_geall(xv,desc_a,info)
nlr = psb_cd_get_local_rows(desc_a)
nlr = desc_a%get_local_rows()
call psb_barrier(ictxt)
talc = psb_wtime()-t0
@ -675,12 +682,13 @@ contains
call psb_amx(ictxt,ttot)
if(iam == psb_root_) then
ch_err = a%get_fmt()
write(*,'("The matrix has been generated and assembled in ",a3," format.")')&
write(psb_out_unit,&
& '("The matrix has been generated and assembled in ",a3," format.")')&
& ch_err(1:3)
write(*,'("-allocation time : ",es12.5)') talc
write(*,'("-coeff. gen. time : ",es12.5)') tgen
write(*,'("-assembly time : ",es12.5)') tasb
write(*,'("-total time : ",es12.5)') ttot
write(psb_out_unit,'("-allocation time : ",es12.5)') talc
write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
write(psb_out_unit,'("-assembly time : ",es12.5)') tasb
write(psb_out_unit,'("-total time : ",es12.5)') ttot
end if
call psb_erractionrestore(err_act)

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

@ -8,11 +8,11 @@ CSR ! Storage format CSR COO JAD
1.d-9 ! EPS
3L-M-RAS-I-D4 ! Longer descriptive name for preconditioner (up to 20 chars)
ML ! Preconditioner NONE JACOBI BJAC AS ML
1 ! Number of overlap layers for AS preconditioner at finest level
0 ! 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
01 ! Level-set N for ILU(N)
1 ! Level-set N for ILU(N)
1.d-5 ! Threshold T for ILU(T,P)
1 ! Smoother/Jacobi sweeps
JACOBI ! Smoother type JACOBI BJAC AS; ignored for non-ML
@ -20,10 +20,10 @@ JACOBI ! Smoother type JACOBI BJAC AS; ignored for non-ML
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED
DEC ! Type of aggregation DEC SYMDEC GLB
MULT ! Type of multilevel correction: ADD MULT
TWOSIDE ! Side of correction PRE POST TWOSIDE (ignored for ADD)
POST ! Side of correction PRE POST TWOSIDE (ignored for ADD)
DIST ! Coarse level: matrix distribution DIST REPL
JACOBI ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST
DSCALE ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
ILU ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST
1 ! Coarse level: Level-set N for ILU(N)
1.d-4 ! Coarse level: Threshold T for ILU(T,P)
4 ! Coarse level: Number of Jacobi sweeps

Write Preview
Loading…
Cancel
Save