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/mld_d_as_smoother.f90
 mlprec/mld_d_diag_solver.f90
 mlprec/mld_d_id_solver.f90
 mlprec/mld_d_ilu_solver.f90
 mlprec/mld_d_inner_mod.f90
 mlprec/mld_d_jac_smoother.f90
 mlprec/mld_d_prec_mod.f90
 mlprec/mld_d_prec_type.f90
 mlprec/mld_d_slu_solver.f90
 mlprec/mld_d_sludist_solver.f90
 mlprec/mld_d_umf_solver.f90
 mlprec/mld_dmlprec_bld.f90
 mlprec/mld_dprecaply.f90
 mlprec/mld_dprecbld.f90

Started inclusion of _vect methods.
stopcriterion
Salvatore Filippone 13 years ago
parent 05a910d78e
commit 859ba30c4c
14 changed files with 1353 additions and 167 deletions
Show Stats Download Patch File Download Diff File

464
mlprec/mld_d_as_smoother.f90
Unescape Escape View File

@ -58,7 +58,8 @@ module mld_d_as_smoother
procedure, pass(sm) :: check => d_as_smoother_check
procedure, pass(sm) :: dump => d_as_smoother_dmp
procedure, pass(sm) :: build => d_as_smoother_bld
procedure, pass(sm) :: apply => d_as_smoother_apply
procedure, pass(sm) :: apply_v => d_as_smoother_apply_vect
procedure, pass(sm) :: apply_a => d_as_smoother_apply
procedure, pass(sm) :: free => d_as_smoother_free
procedure, pass(sm) :: seti => d_as_smoother_seti
procedure, pass(sm) :: setc => d_as_smoother_setc
@ -74,7 +75,7 @@ module mld_d_as_smoother
& d_as_smoother_setc, d_as_smoother_setr,&
& d_as_smoother_descr, d_as_smoother_sizeof, &
& d_as_smoother_check, d_as_smoother_default,&
& d_as_smoother_dmp
& d_as_smoother_dmp, d_as_smoother_apply_vect
character(len=6), parameter, private :: &
& restrict_names(0:4)=(/'none ','halo ',' ',' ',' '/)
@ -151,6 +152,442 @@ contains
return
end subroutine d_as_smoother_check
subroutine d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_as_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
integer, intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
integer :: n_row,n_col, nrow_d, i
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_
character(len=20) :: name='d_as_smoother_apply', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = desc_data%get_context()
call psb_info(ictxt,me,np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T')
case('C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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)
if ((6*isz) <= size(work)) then
ww => work(1:isz)
tx => work(isz+1:2*isz)
ty => work(2*isz+1:3*isz)
aux => work(3*isz+1:)
else if ((4*isz) <= size(work)) then
aux => work(1:)
allocate(ww(isz),tx(isz),ty(isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,i_err=(/3*isz,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
else if ((3*isz) <= size(work)) then
ww => work(1:isz)
tx => work(isz+1:2*isz)
ty => work(2*isz+1:3*isz)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,i_err=(/4*isz,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
else
allocate(ww(isz),tx(isz),ty(isz),&
&aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,i_err=(/4*isz,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
if ((sm%novr == 0).and.(sweeps == 1)) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi.
!
call sm%sv%apply(alpha,x,beta,y,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
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)
!!$
end if
if ((6*isz) <= size(work)) then
else if ((4*isz) <= size(work)) then
deallocate(ww,tx,ty)
else if ((3*isz) <= size(work)) then
deallocate(aux)
else
deallocate(ww,aux,tx,ty)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_as_smoother_apply_vect
subroutine d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
@ -178,7 +615,8 @@ contains
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case('T')
case('C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
@ -585,7 +1023,7 @@ contains
end subroutine d_as_smoother_apply
subroutine d_as_smoother_bld(a,desc_a,sm,upd,info,mold)
subroutine d_as_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
use psb_base_mod
@ -597,7 +1035,9 @@ contains
class(mld_d_as_smoother_type), intent(inout) :: sm
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
type(psb_dspmat_type) :: blck, atmp
integer :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros
@ -688,7 +1128,8 @@ contains
End if
if (info == psb_success_) &
& call sm%sv%build(a,sm%desc_data,upd,info,blck,mold=mold)
& call sm%sv%build(a,sm%desc_data,upd,info,&
& blck,amold=amold,vmold=vmold)
nrow_a = a%get_nrows()
n_row = sm%desc_data%get_local_rows()
@ -699,9 +1140,16 @@ contains
if (info == psb_success_) call blck%csclip(atmp,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
if (info == psb_success_) call psb_rwextd(n_row,sm%nd,info,b=atmp)
if (info == psb_success_) call sm%nd%cscnv(info,&
& type='csr',dupl=psb_dupl_add_)
if (info == psb_success_) then
if (present(amold)) then
call sm%nd%cscnv(info,&
& mold=amold,dupl=psb_dupl_add_)
else
call sm%nd%cscnv(info,&
& type='csr',dupl=psb_dupl_add_)
end if
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4')
goto 9999

95
mlprec/mld_d_diag_solver.f90
Unescape Escape View File

@ -48,9 +48,11 @@ module mld_d_diag_solver
use mld_d_prec_type
type, extends(mld_d_base_solver_type) :: mld_d_diag_solver_type
type(psb_d_vect_type), allocatable :: dv
real(psb_dpk_), allocatable :: d(:)
contains
procedure, pass(sv) :: build => d_diag_solver_bld
procedure, pass(sv) :: apply_v => d_diag_solver_apply_vect
procedure, pass(sv) :: apply_a => d_diag_solver_apply
procedure, pass(sv) :: free => d_diag_solver_free
procedure, pass(sv) :: seti => d_diag_solver_seti
@ -69,6 +71,79 @@ module mld_d_diag_solver
contains
subroutine d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_diag_solver_type), intent(inout) :: sv
type(psb_d_vect_type), intent(inout) :: x
type(psb_d_vect_type), intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
integer :: n_row,n_col
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me,i, err_act
character :: trans_
character(len=20) :: name='d_diag_solver_apply'
call psb_erractionsave(err_act)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols()
if (x%get_nrows() < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (y%get_nrows() < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
if (.not.allocated(sv%dv)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (sv%dv%get_nrows() < n_row) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
call y%mlt(alpha,sv%dv,x,beta,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_solver_apply_vect
subroutine d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
@ -189,7 +264,7 @@ contains
end subroutine d_diag_solver_apply
subroutine d_diag_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -202,7 +277,8 @@ contains
character, intent(in) :: upd
integer, intent(out) :: info
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
integer :: n_row,n_col, nrow_a, nztota
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@ -251,7 +327,20 @@ contains
sv%d(i) = done/sv%d(i)
end if
end do
allocate(sv%dv,stat=info)
if (info == psb_success_) then
if (present(vmold)) then
allocate(sv%dv%v,mold=vmold,stat=info)
else
allocate(psb_d_base_vect_type :: sv%dv%v,stat=info)
end if
end if
if (info == psb_success_) then
call sv%dv%bld(sv%d)
else
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'

56
mlprec/mld_d_id_solver.f90
Unescape Escape View File

@ -50,6 +50,7 @@ module mld_d_id_solver
type, extends(mld_d_base_solver_type) :: mld_d_id_solver_type
contains
procedure, pass(sv) :: build => d_id_solver_bld
procedure, pass(sv) :: apply_v => d_id_solver_apply_vect
procedure, pass(sv) :: apply_a => d_id_solver_apply
procedure, pass(sv) :: free => d_id_solver_free
procedure, pass(sv) :: seti => d_id_solver_seti
@ -68,6 +69,53 @@ module mld_d_id_solver
contains
subroutine d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_id_solver_type), intent(inout) :: sv
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
integer :: n_row,n_col
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me,i, err_act
character :: trans_
character(len=20) :: name='d_id_solver_apply'
call psb_erractionsave(err_act)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T')
case('C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_id_solver_apply_vect
subroutine d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
@ -92,7 +140,8 @@ contains
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case('T')
case('C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
@ -113,7 +162,7 @@ contains
end subroutine d_id_solver_apply
subroutine d_id_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_id_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -125,8 +174,9 @@ contains
class(mld_d_id_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
integer :: n_row,n_col, nrow_a, nztota
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)

192
mlprec/mld_d_ilu_solver.f90
Unescape Escape View File

@ -51,11 +51,13 @@ module mld_d_ilu_solver
type, extends(mld_d_base_solver_type) :: mld_d_ilu_solver_type
type(psb_dspmat_type) :: l, u
real(psb_dpk_), allocatable :: d(:)
type(psb_d_vect_type), allocatable :: dv
integer :: fact_type, fill_in
real(psb_dpk_) :: thresh
contains
procedure, pass(sv) :: dump => d_ilu_solver_dmp
procedure, pass(sv) :: build => d_ilu_solver_bld
procedure, pass(sv) :: apply_v => d_ilu_solver_apply_vect
procedure, pass(sv) :: apply_a => d_ilu_solver_apply
procedure, pass(sv) :: free => d_ilu_solver_free
procedure, pass(sv) :: seti => d_ilu_solver_seti
@ -63,6 +65,7 @@ module mld_d_ilu_solver
procedure, pass(sv) :: setr => d_ilu_solver_setr
procedure, pass(sv) :: descr => d_ilu_solver_descr
procedure, pass(sv) :: sizeof => d_ilu_solver_sizeof
procedure, pass(sv) :: get_nzeros => d_ilu_solver_get_nzeros
procedure, pass(sv) :: default => d_ilu_solver_default
end type mld_d_ilu_solver_type
@ -71,7 +74,8 @@ module mld_d_ilu_solver
& d_ilu_solver_free, d_ilu_solver_seti, &
& d_ilu_solver_setc, d_ilu_solver_setr,&
& d_ilu_solver_descr, d_ilu_solver_sizeof, &
& d_ilu_solver_default, d_ilu_solver_dmp
& d_ilu_solver_default, d_ilu_solver_dmp, &
& d_ilu_solver_apply_vect, d_ilu_solver_get_nzeros
character(len=15), parameter, private :: &
@ -143,6 +147,142 @@ contains
end subroutine d_ilu_solver_check
subroutine d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_ilu_solver_type), intent(inout) :: sv
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
integer :: n_row,n_col
type(psb_d_vect_type) :: wv
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
integer :: ictxt,np,me,i, err_act
character :: trans_
character(len=20) :: name='d_ilu_solver_apply'
call psb_erractionsave(err_act)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T')
case('C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols()
if (x%get_nrows() < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/))
goto 9999
end if
if (y%get_nrows() < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/))
goto 9999
end if
if (.not.allocated(sv%dv)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: DV")
goto 9999
end if
if (sv%dv%get_nrows() < n_row) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: DV")
goto 9999
end if
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
endif
if (info == psb_success_) allocate(wv%v,mold=x%v)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
call wv%bld(n_col)
select case(trans_)
case('N')
call psb_spsm(done,sv%l,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%u,wv,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
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,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
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,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
goto 9999
endif
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_ilu_solver_apply_vect
subroutine d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
@ -207,7 +347,14 @@ contains
if (info == psb_success_) call psb_spsm(alpha,sv%u,ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T','C')
case('T')
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,&
& 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
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,&
@ -246,7 +393,7 @@ contains
end subroutine d_ilu_solver_apply
subroutine d_ilu_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_ilu_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -258,8 +405,9 @@ contains
class(mld_d_ilu_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
integer :: n_row,n_col, nrow_a, nztota
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@ -301,13 +449,22 @@ contains
endif
if (.not.allocated(sv%d)) then
allocate(sv%d(n_row),stat=info)
endif
if (info == psb_success_) then
allocate(sv%dv, stat=info)
if (info == 0) then
if (present(vmold)) then
allocate(sv%dv%v,mold=vmold,stat=info)
else
allocate(psb_d_base_vect_type :: sv%dv%v,stat=info)
end if
end if
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
endif
endif
select case(sv%fact_type)
@ -403,10 +560,11 @@ contains
call sv%l%trim()
call sv%u%set_asb()
call sv%u%trim()
call sv%dv%bld(sv%d)
if (present(mold)) then
call sv%l%cscnv(info,mold=mold)
call sv%u%cscnv(info,mold=mold)
if (present(amold)) then
call sv%l%cscnv(info,mold=amold)
call sv%u%cscnv(info,mold=amold)
end if
if (debug_level >= psb_debug_outer_) &
@ -626,6 +784,22 @@ contains
return
end subroutine d_ilu_solver_descr
function d_ilu_solver_get_nzeros(sv) result(val)
use psb_base_mod
implicit none
! Arguments
class(mld_d_ilu_solver_type), intent(in) :: sv
integer(psb_long_int_k_) :: val
integer :: i
val = 0
if (allocated(sv%dv)) val = val + sv%dv%get_nrows()
val = val + sv%l%get_nzeros()
val = val + sv%u%get_nzeros()
return
end function d_ilu_solver_get_nzeros
function d_ilu_solver_sizeof(sv) result(val)
use psb_base_mod
implicit none

7
mlprec/mld_d_inner_mod.f90
Unescape Escape View File

@ -50,16 +50,17 @@ module mld_d_inner_mod
interface mld_mlprec_bld
subroutine mld_dmlprec_bld(a,desc_a,prec,info, mold)
subroutine mld_dmlprec_bld(a,desc_a,prec,info, amold, vmold)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, &
& psb_dpk_, psb_d_base_sparse_mat
& psb_dpk_, psb_d_base_sparse_mat, psb_d_base_vect_type
use mld_d_prec_type, only : mld_dprec_type
implicit none
type(psb_dspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type), intent(inout), target :: prec
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
!!$ character, intent(in),optional :: upd
end subroutine mld_dmlprec_bld
end interface mld_mlprec_bld

182
mlprec/mld_d_jac_smoother.f90
Unescape Escape View File

@ -55,7 +55,8 @@ module mld_d_jac_smoother
integer :: nnz_nd_tot
contains
procedure, pass(sm) :: build => d_jac_smoother_bld
procedure, pass(sm) :: apply => d_jac_smoother_apply
procedure, pass(sm) :: apply_v => d_jac_smoother_apply_vect
procedure, pass(sm) :: apply_a => d_jac_smoother_apply
procedure, pass(sm) :: free => d_jac_smoother_free
procedure, pass(sm) :: seti => d_jac_smoother_seti
procedure, pass(sm) :: setc => d_jac_smoother_setc
@ -68,12 +69,173 @@ module mld_d_jac_smoother
private :: d_jac_smoother_bld, d_jac_smoother_apply, &
& d_jac_smoother_free, d_jac_smoother_seti, &
& d_jac_smoother_setc, d_jac_smoother_setr,&
& d_jac_smoother_descr, d_jac_smoother_sizeof
& d_jac_smoother_descr, d_jac_smoother_sizeof, &
& d_jac_smoother_apply_vect
contains
subroutine d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_jac_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
integer, intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
integer :: n_row,n_col
type(psb_d_vect_type) :: tx, ty
real(psb_dpk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me,i, err_act
character :: trans_
character(len=20) :: name='d_jac_smoother_apply'
call psb_erractionsave(err_act)
info = psb_success_
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T')
case('C')
case default
call psb_errpush(psb_err_iarg_invalid_i_,name)
goto 9999
end select
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
n_row = desc_data%get_local_rows()
n_col = desc_data%get_local_cols()
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then
call sm%sv%apply(alpha,x,beta,y,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
else if (sweeps > 1) then
!
!
! Apply multiple sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system.
!
!
allocate(tx%v,mold=x%v,stat=info)
if (info == psb_success_) allocate(ty%v,mold=x%v,stat=info)
if (info /= psb_success_) then
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*n_col,0,0,0,0/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
call tx%bld(x%get_nrows())
call ty%bld(x%get_nrows())
do i=1, sweeps
!
! 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,x,dzero,ty,desc_data,info)
call psb_spmm(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(done,ty,dzero,tx,desc_data,trans_,aux,info)
if (info /= psb_success_) exit
end do
if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
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
call tx%free(info)
if (info == psb_success_) call ty%free(info)
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,a_err='final cleanup 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
endif
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_jac_smoother_apply_vect
subroutine d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
@ -228,7 +390,7 @@ contains
end subroutine d_jac_smoother_apply
subroutine d_jac_smoother_bld(a,desc_a,sm,upd,info,mold)
subroutine d_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
use psb_base_mod
use mld_d_diag_solver
@ -240,7 +402,8 @@ contains
class(mld_d_jac_smoother_type), intent(inout) :: sm
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
integer :: n_row,n_col, nrow_a, nztota, nzeros
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@ -268,15 +431,22 @@ contains
call a%csclip(sm%nd,info,&
& jmin=nrow_a+1,rscale=.false.,cscale=.false.)
end select
if (info == psb_success_) call sm%nd%cscnv(info,&
if (info == psb_success_) then
if (present(amold)) then
call sm%nd%cscnv(info,&
& mold=amold,dupl=psb_dupl_add_)
else
call sm%nd%cscnv(info,&
& type='csr',dupl=psb_dupl_add_)
endif
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='clip & psb_spcnv csr 4')
goto 9999
end if
call sm%sv%build(a,desc_a,upd,info,mold=mold)
call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='solver build')

7
mlprec/mld_d_prec_mod.f90
Unescape Escape View File

@ -111,16 +111,17 @@ module mld_d_prec_mod
end interface
interface mld_precbld
subroutine mld_dprecbld(a,desc_a,prec,info,mold)
subroutine mld_dprecbld(a,desc_a,prec,info,amold,vmold)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, &
& psb_dpk_, psb_d_base_sparse_mat
& psb_dpk_, psb_d_base_sparse_mat, psb_d_base_vect_type
use mld_d_prec_type, only : mld_dprec_type
implicit none
type(psb_dspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type), intent(inout), target :: prec
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
!!$ character, intent(in),optional :: upd
end subroutine mld_dprecbld
end interface

170
mlprec/mld_d_prec_type.f90
Unescape Escape View File

@ -60,6 +60,7 @@
module mld_d_prec_type
use mld_base_prec_type
use psb_base_mod, only : psb_d_vect_type, psb_d_base_vect_type
!
! Type: mld_Tprec_type.
!
@ -191,6 +192,7 @@ module mld_d_prec_type
procedure, pass(sv) :: default => d_base_solver_default
procedure, pass(sv) :: descr => d_base_solver_descr
procedure, pass(sv) :: sizeof => d_base_solver_sizeof
procedure, pass(sv) :: get_nzeros => d_base_solver_get_nzeros
end type mld_d_base_solver_type
type mld_d_base_smoother_type
@ -199,7 +201,9 @@ module mld_d_prec_type
procedure, pass(sm) :: check => d_base_smoother_check
procedure, pass(sm) :: dump => d_base_smoother_dmp
procedure, pass(sm) :: build => d_base_smoother_bld
procedure, pass(sm) :: apply => d_base_smoother_apply
procedure, pass(sm) :: apply_v => d_base_smoother_apply_vect
procedure, pass(sm) :: apply_a => d_base_smoother_apply
generic, public :: apply => apply_a, apply_v
procedure, pass(sm) :: free => d_base_smoother_free
procedure, pass(sm) :: seti => d_base_smoother_seti
procedure, pass(sm) :: setc => d_base_smoother_setc
@ -208,6 +212,7 @@ module mld_d_prec_type
procedure, pass(sm) :: default => d_base_smoother_default
procedure, pass(sm) :: descr => d_base_smoother_descr
procedure, pass(sm) :: sizeof => d_base_smoother_sizeof
procedure, pass(sm) :: get_nzeros => d_base_smoother_get_nzeros
end type mld_d_base_smoother_type
type mld_donelev_type
@ -227,6 +232,7 @@ module mld_d_prec_type
procedure, pass(lv) :: setr => d_base_onelev_setr
procedure, pass(lv) :: setc => d_base_onelev_setc
generic, public :: set => seti, setr, setc
procedure, pass(lv) :: get_nzeros => d_base_onelev_get_nzeros
end type mld_donelev_type
type, extends(psb_dprec_type) :: mld_dprec_type
@ -234,11 +240,13 @@ module mld_d_prec_type
real(psb_dpk_) :: op_complexity=-done
type(mld_donelev_type), allocatable :: precv(:)
contains
procedure, pass(prec) :: d_apply2_vect => mld_d_apply2_vect
procedure, pass(prec) :: d_apply2v => mld_d_apply2v
procedure, pass(prec) :: d_apply1v => mld_d_apply1v
procedure, pass(prec) :: dump => mld_d_dump
procedure, pass(prec) :: get_complexity => mld_d_get_compl
procedure, pass(prec) :: cmp_complexity => mld_d_cmp_compl
procedure, pass(prec) :: get_nzeros => mld_d_get_nzeros
end type mld_dprec_type
private :: d_base_solver_bld, d_base_solver_apply, &
@ -246,18 +254,20 @@ module mld_d_prec_type
& d_base_solver_setc, d_base_solver_setr, &
& d_base_solver_descr, d_base_solver_sizeof, &
& d_base_solver_default, d_base_solver_check,&
& d_base_solver_dmp, &
& d_base_solver_dmp, d_base_solver_apply_vect, &
& d_base_smoother_bld, d_base_smoother_apply, &
& d_base_smoother_free, d_base_smoother_seti, &
& d_base_smoother_setc, d_base_smoother_setr,&
& d_base_smoother_descr, d_base_smoother_sizeof, &
& d_base_smoother_default, d_base_smoother_check, &
& d_base_smoother_dmp, &
& d_base_smoother_dmp, d_base_smoother_apply_vect, &
& d_base_onelev_seti, d_base_onelev_setc, &
& d_base_onelev_setr, d_base_onelev_check, &
& d_base_onelev_default, d_base_onelev_dump, &
& d_base_onelev_descr, mld_d_dump, &
& mld_d_get_compl, mld_d_cmp_compl
& mld_d_get_compl, mld_d_cmp_compl,&
& mld_d_get_nzeros, d_base_onelev_get_nzeros, &
& d_base_smoother_get_nzeros, d_base_solver_get_nzeros
!
@ -282,6 +292,18 @@ module mld_d_prec_type
end interface
interface mld_precaply
subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, &
& psb_dpk_, psb_d_vect_type
import mld_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(inout) :: prec
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
end subroutine mld_dprecaply_vect
subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work)
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_
import mld_dprec_type
@ -309,6 +331,48 @@ contains
! Function returning the size of the mld_prec_type data structure
!
function d_base_solver_get_nzeros(sv) result(val)
implicit none
class(mld_d_base_solver_type), intent(in) :: sv
integer(psb_long_int_k_) :: val
integer :: i
val = 0
end function d_base_solver_get_nzeros
function d_base_smoother_get_nzeros(sm) result(val)
implicit none
class(mld_d_base_smoother_type), intent(in) :: sm
integer(psb_long_int_k_) :: val
integer :: i
val = 0
if (allocated(sm%sv)) &
& val = sm%sv%get_nzeros()
end function d_base_smoother_get_nzeros
function d_base_onelev_get_nzeros(lv) result(val)
implicit none
class(mld_donelev_type), intent(in) :: lv
integer(psb_long_int_k_) :: val
integer :: i
val = 0
if (allocated(lv%sm)) &
& val = lv%sm%get_nzeros()
end function d_base_onelev_get_nzeros
function mld_d_get_nzeros(prec) result(val)
implicit none
class(mld_dprec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
integer :: i
val = 0
if (allocated(prec%precv)) then
do i=1, size(prec%precv)
val = val + prec%precv(i)%get_nzeros()
end do
end if
end function mld_d_get_nzeros
function mld_dprec_sizeof(prec) result(val)
implicit none
type(mld_dprec_type), intent(in) :: prec
@ -451,6 +515,10 @@ contains
endif
call p%precv(1)%sm%descr(info,iout=iout_)
if (nlev == 1) then
if (p%precv(1)%parms%sweeps > 1) then
write(iout_,*) ' Number of sweeps : ',&
& p%precv(1)%parms%sweeps
end if
write(iout_,*)
return
end if
@ -688,6 +756,47 @@ contains
end subroutine d_base_smoother_apply
subroutine d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_base_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1),intent(in) :: trans
integer, intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:)
integer, intent(out) :: info
Integer :: err_act
character(len=20) :: name='d_base_smoother_apply'
call psb_erractionsave(err_act)
info = psb_success_
if (allocated(sm%sv)) then
call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
else
info = 1121
endif
if (info /= psb_success_) then
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_smoother_apply_vect
subroutine d_base_smoother_check(sm,info)
use psb_base_mod
@ -830,7 +939,7 @@ contains
return
end subroutine d_base_smoother_setr
subroutine d_base_smoother_bld(a,desc_a,sm,upd,info,mold)
subroutine d_base_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
use psb_base_mod
@ -842,7 +951,8 @@ contains
class(mld_d_base_smoother_type), intent(inout) :: sm
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
Integer :: err_act
character(len=20) :: name='d_base_smoother_bld'
@ -850,7 +960,7 @@ contains
info = psb_success_
if (allocated(sm%sv)) then
call sm%sv%build(a,desc_a,upd,info,mold=mold)
call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold)
else
info = 1121
call psb_errpush(info,name)
@ -989,7 +1099,6 @@ contains
end subroutine d_base_smoother_default
subroutine d_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
@ -1023,11 +1132,10 @@ contains
end subroutine d_base_solver_apply
subroutine d_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
use psb_base_mod
type(psb_desc_type), intent(in) :: desc_data
class(mld_d_base_solver_type), intent(in) :: sv
class(mld_d_base_solver_type), intent(inout) :: sv
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
real(psb_dpk_),intent(in) :: alpha,beta
@ -1057,7 +1165,7 @@ contains
end subroutine d_base_solver_apply_vect
subroutine d_base_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_base_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -1070,7 +1178,8 @@ contains
character, intent(in) :: upd
integer, intent(out) :: info
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
Integer :: err_act
character(len=20) :: name='d_base_solver_bld'
@ -1288,6 +1397,43 @@ contains
end subroutine d_base_solver_default
subroutine mld_d_apply2_vect(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(mld_dprec_type), intent(inout) :: prec
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act
character(len=20) :: name='d_prec_apply'
call psb_erractionsave(err_act)
select type(prec)
type is (mld_dprec_type)
call mld_precaply(prec,x,y,desc_data,info,trans,work)
class default
info = psb_err_missing_override_method_
call psb_errpush(info,name)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_d_apply2_vect
subroutine mld_d_apply2v(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data

5
mlprec/mld_d_slu_solver.f90
Unescape Escape View File

@ -185,7 +185,7 @@ contains
end subroutine d_slu_solver_apply
subroutine d_slu_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_slu_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -197,8 +197,9 @@ contains
class(mld_d_slu_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
type(psb_dspmat_type) :: atmp
type(psb_d_csr_sparse_mat) :: acsr

5
mlprec/mld_d_sludist_solver.f90
Unescape Escape View File

@ -185,7 +185,7 @@ contains
end subroutine d_sludist_solver_apply
subroutine d_sludist_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_sludist_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -197,8 +197,9 @@ contains
class(mld_d_sludist_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
type(psb_dspmat_type) :: atmp
type(psb_d_csr_sparse_mat) :: acsr

5
mlprec/mld_d_umf_solver.f90
Unescape Escape View File

@ -185,7 +185,7 @@ contains
end subroutine d_umf_solver_apply
subroutine d_umf_solver_bld(a,desc_a,sv,upd,info,b,mold)
subroutine d_umf_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
use psb_base_mod
@ -197,8 +197,9 @@ contains
class(mld_d_umf_solver_type), intent(inout) :: sv
character, intent(in) :: upd
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
! Local variables
type(psb_dspmat_type) :: atmp
type(psb_d_csc_sparse_mat) :: acsc

11
mlprec/mld_dmlprec_bld.f90
Unescape Escape View File

@ -63,7 +63,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_dmlprec_bld(a,desc_a,p,info,mold)
subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold)
use psb_base_mod
use mld_d_inner_mod, mld_protect_name => mld_dmlprec_bld
@ -76,7 +76,8 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,mold)
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type),intent(inout),target :: p
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
!!$ character, intent(in), optional :: upd
! Local Variables
@ -309,10 +310,10 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,mold)
! Test version for beginning of OO stuff.
!
call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,&
& 'F',info,mold=mold)
& 'F',info,amold=amold,vmold=vmold)
if ((info == psb_success_).and.(i>1).and.(present(mold))) then
call psb_map_cscnv(p%precv(i)%map,info,mold=mold)
if ((info == psb_success_).and.(i>1).and.(present(amold))) then
call psb_map_cscnv(p%precv(i)%map,info,mold=amold)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

102
mlprec/mld_dprecaply.f90
Unescape Escape View File

@ -263,3 +263,105 @@ subroutine mld_dprecaply1(prec,x,desc_data,info,trans)
end if
return
end subroutine mld_dprecaply1
subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work)
use psb_base_mod
use mld_d_inner_mod
implicit none
! Arguments
type(psb_desc_type),intent(in) :: desc_data
type(mld_dprec_type), intent(inout) :: prec
type(psb_d_vect_type),intent(inout) :: x
type(psb_d_vect_type),intent(inout) :: y
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
! Local variables
character :: trans_
real(psb_dpk_), pointer :: work_(:)
integer :: ictxt,np,me,err_act,iwsz
character(len=20) :: name
name='mld_dprecaply'
info = psb_success_
call psb_erractionsave(err_act)
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
if (present(trans)) then
trans_=psb_toupper(trans)
else
trans_='N'
end if
if (present(work)) then
work_ => work
else
iwsz = max(1,4*psb_cd_get_local_cols(desc_data))
allocate(work_(iwsz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,i_err=(/iwsz,0,0,0,0/),&
&a_err='real(psb_dpk_)')
goto 9999
end if
end if
if (.not.(allocated(prec%precv))) then
!! Error 1: should call mld_dprecbld
info=3112
call psb_errpush(info,name)
goto 9999
end if
if (size(prec%precv) >1) then
!
! 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_
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply')
goto 9999
end if
else if (size(prec%precv) == 1) then
!
! Number of levels = 1: apply the base preconditioner
!
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/))
goto 9999
endif
! If the original distribution has an overlap we should fix that.
!!$ call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (present(work)) then
else
deallocate(work_)
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
end subroutine mld_dprecaply_vect

9
mlprec/mld_dprecbld.f90
Unescape Escape View File

@ -58,7 +58,7 @@
! info - integer, output.
! Error code.
!
subroutine mld_dprecbld(a,desc_a,p,info,mold)
subroutine mld_dprecbld(a,desc_a,p,info,amold,vmold)
use psb_base_mod
use mld_d_inner_mod
@ -71,7 +71,8 @@ subroutine mld_dprecbld(a,desc_a,p,info,mold)
type(psb_desc_type), intent(in), target :: desc_a
type(mld_dprec_type),intent(inout), target :: p
integer, intent(out) :: info
class(psb_d_base_sparse_mat), intent(in), optional :: mold
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
!!$ character, intent(in), optional :: upd
! Local Variables
@ -171,7 +172,7 @@ subroutine mld_dprecbld(a,desc_a,p,info,mold)
goto 9999
endif
call p%precv(1)%sm%build(a,desc_a,upd_,info,mold=mold)
call p%precv(1)%sm%build(a,desc_a,upd_,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='One level preconditioner build.')
@ -185,7 +186,7 @@ subroutine mld_dprecbld(a,desc_a,p,info,mold)
!
! Build the multilevel preconditioner
!
call mld_mlprec_bld(a,desc_a,p,info,mold=mold)
call mld_mlprec_bld(a,desc_a,p,info,amold=amold,vmold=vmold)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&

Write Preview
Loading…
Cancel
Save