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/impl/smoother/Makefile
 mlprec/impl/smoother/mld_c_as_smoother_apply.f90
 mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_c_as_smoother_prol_a.f90
 mlprec/impl/smoother/mld_c_as_smoother_prol_v.f90
 mlprec/impl/smoother/mld_c_as_smoother_restr_a.f90
 mlprec/impl/smoother/mld_c_as_smoother_restr_v.f90
 mlprec/impl/smoother/mld_d_as_smoother_apply.f90
 mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_d_as_smoother_prol_a.f90
 mlprec/impl/smoother/mld_d_as_smoother_prol_v.f90
 mlprec/impl/smoother/mld_d_as_smoother_restr_a.f90
 mlprec/impl/smoother/mld_d_as_smoother_restr_v.f90
 mlprec/impl/smoother/mld_s_as_smoother_apply.f90
 mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_s_as_smoother_prol_a.f90
 mlprec/impl/smoother/mld_s_as_smoother_prol_v.f90
 mlprec/impl/smoother/mld_s_as_smoother_restr_a.f90
 mlprec/impl/smoother/mld_s_as_smoother_restr_v.f90
 mlprec/impl/smoother/mld_z_as_smoother_apply.f90
 mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_z_as_smoother_prol_a.f90
 mlprec/impl/smoother/mld_z_as_smoother_prol_v.f90
 mlprec/impl/smoother/mld_z_as_smoother_restr_a.f90
 mlprec/impl/smoother/mld_z_as_smoother_restr_v.f90
 mlprec/mld_c_as_smoother.f90
 mlprec/mld_d_as_smoother.f90
 mlprec/mld_s_as_smoother.f90
 mlprec/mld_z_as_smoother.f90


Refactored application of AS smoothers.
stopcriterion
Salvatore Filippone 9 years ago
parent da6bde2130
commit 999c7eaf97
29 changed files with 3552 additions and 3104 deletions
Show Stats Download Patch File Download Diff File

16
mlprec/impl/smoother/Makefile
Unescape Escape View File

@ -22,6 +22,10 @@ mld_c_as_smoother_free.o \
mld_c_as_smoother_setc.o \
mld_c_as_smoother_seti.o \
mld_c_as_smoother_setr.o \
mld_c_as_smoother_prol_a.o \
mld_c_as_smoother_prol_v.o \
mld_c_as_smoother_restr_a.o \
mld_c_as_smoother_restr_v.o \
mld_c_base_smoother_apply.o \
mld_c_base_smoother_apply_vect.o \
mld_c_base_smoother_bld.o \
@ -57,6 +61,10 @@ mld_d_as_smoother_free.o \
mld_d_as_smoother_setc.o \
mld_d_as_smoother_seti.o \
mld_d_as_smoother_setr.o \
mld_d_as_smoother_prol_a.o \
mld_d_as_smoother_prol_v.o \
mld_d_as_smoother_restr_a.o \
mld_d_as_smoother_restr_v.o \
mld_d_base_smoother_apply.o \
mld_d_base_smoother_apply_vect.o \
mld_d_base_smoother_bld.o \
@ -92,6 +100,10 @@ mld_s_as_smoother_free.o \
mld_s_as_smoother_setc.o \
mld_s_as_smoother_seti.o \
mld_s_as_smoother_setr.o \
mld_s_as_smoother_prol_a.o \
mld_s_as_smoother_prol_v.o \
mld_s_as_smoother_restr_a.o \
mld_s_as_smoother_restr_v.o \
mld_s_base_smoother_apply.o \
mld_s_base_smoother_apply_vect.o \
mld_s_base_smoother_bld.o \
@ -127,6 +139,10 @@ mld_z_as_smoother_free.o \
mld_z_as_smoother_setc.o \
mld_z_as_smoother_seti.o \
mld_z_as_smoother_setr.o \
mld_z_as_smoother_prol_a.o \
mld_z_as_smoother_prol_v.o \
mld_z_as_smoother_restr_a.o \
mld_z_as_smoother_restr_v.o \
mld_z_base_smoother_apply.o \
mld_z_base_smoother_apply_vect.o \
mld_z_base_smoother_bld.o \

457
mlprec/impl/smoother/mld_c_as_smoother_apply.f90
Unescape Escape View File

@ -54,15 +54,15 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
complex(psb_spk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
complex(psb_spk_), pointer :: ww(:), aux(:)
complex(psb_spk_), allocatable :: tx(:),ty(:)
complex(psb_spk_), pointer :: aux(:)
complex(psb_spk_), allocatable :: tx(:),ty(:), ww(:)
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='c_as_smoother_apply', ch_err
call psb_erractionsave(err_act)
info = psb_success_
info = psb_success_
ictxt = desc_data%get_context()
call psb_info(ictxt,me,np)
@ -89,55 +89,26 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
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)
if ((6*isz) <= size(work)) then
ww => work(1:isz)
aux => work(3*isz+1:)
else if ((4*isz) <= size(work)) then
if ((4*isz) <= size(work)) then
aux => work(1:)
allocate(ww(isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/3*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
else if ((3*isz) <= size(work)) then
ww => work(1:isz)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
else
allocate(ww(isz), aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
endif
if (sweeps == 0) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -147,361 +118,86 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info)
call psb_geasb(ty,sm%desc_data,info)
call psb_geasb(ww,sm%desc_data,info)
call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
select case (init_)
case('Z')
tx(:) = czero
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(cone,y,czero,tx,desc_data,info)
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,initu,czero,tx,desc_data,info)
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,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
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=izero)
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(cone,tx,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
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.
do i=1, sweeps-1
!
! 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.
!
select case (init_)
case('Z')
ty = czero
case('Y')
call psb_geaxpby(cone,y,czero,ty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,initu,czero,ty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
@ -509,17 +205,26 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
!
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
if (info ==0) deallocate(ww,tx,ty,stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

468
mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90
Unescape Escape View File

@ -54,12 +54,11 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
type(psb_c_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_spk_), allocatable :: vx(:)
type(psb_c_vect_type) :: vtx, vty, vww
complex(psb_spk_), pointer :: aux(:)
type(psb_c_vect_type) :: tx, ty, ww
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='c_as_smoother_apply', ch_err
character(len=20) :: name='c_as_smoother_apply_v', ch_err
call psb_erractionsave(err_act)
@ -95,55 +94,21 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
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,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
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,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
if (4*isz <= size(work)) then
aux => work(:)
else
allocate(ww(isz),tx(isz),ty(isz),&
&aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
endif
if (sweeps == 0) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -153,382 +118,115 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call vtx%zero()
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(cone,y,czero,vtx,desc_data,info)
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,initu,czero,vtx,desc_data,info)
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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=izero)
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(cone,vtx,czero,vty,sm%desc_data,trans_,aux,info,init='Y')
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
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
do i=1, sweeps-1
!
! 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.
!
! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system.
!
!
select case (init_)
case('Z')
call vty%zero()
case('Y')
call psb_geaxpby(cone,y,czero,vty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,initu,czero,vty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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(cone,vtx,czero,vww,sm%desc_data,info)
call psb_spmm(-cone,sm%nd,vty,cone,vww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(cone,vww,czero,vty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,vty,beta,y,desc_data,info)
end if
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
call vww%free(info)
call vtx%free(info)
call vty%free(info)
if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

150
mlprec/impl/smoother/mld_c_as_smoother_prol_a.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_c_as_smoother_prol_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_prol_a
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
complex(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='c_as_smther_prol_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_as_smoother_prol_a

150
mlprec/impl/smoother/mld_c_as_smoother_prol_v.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_c_as_smoother_prol_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_prol_v
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
type(psb_c_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='c_as_smther_prol_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_as_smoother_prol_v

169
mlprec/impl/smoother/mld_c_as_smoother_restr_a.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_c_as_smoother_restr_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_restr_a
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
complex(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='c_as_smther_restr_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_as_smoother_restr_a

169
mlprec/impl/smoother/mld_c_as_smoother_restr_v.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_c_as_smoother_restr_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_restr_v
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
type(psb_c_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='c_as_smther_restr_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_as_smoother_restr_v

457
mlprec/impl/smoother/mld_d_as_smoother_apply.f90
Unescape Escape View File

@ -54,15 +54,15 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
real(psb_dpk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
real(psb_dpk_), pointer :: ww(:), aux(:)
real(psb_dpk_), allocatable :: tx(:),ty(:)
real(psb_dpk_), pointer :: aux(:)
real(psb_dpk_), allocatable :: tx(:),ty(:), ww(:)
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='d_as_smoother_apply', ch_err
call psb_erractionsave(err_act)
info = psb_success_
info = psb_success_
ictxt = desc_data%get_context()
call psb_info(ictxt,me,np)
@ -89,55 +89,26 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
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)
if ((6*isz) <= size(work)) then
ww => work(1:isz)
aux => work(3*isz+1:)
else if ((4*isz) <= size(work)) then
if ((4*isz) <= size(work)) then
aux => work(1:)
allocate(ww(isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/3*isz,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
else if ((3*isz) <= size(work)) then
ww => work(1:isz)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
else
allocate(ww(isz), aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
if (sweeps == 0) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -147,361 +118,86 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info)
call psb_geasb(ty,sm%desc_data,info)
call psb_geasb(ww,sm%desc_data,info)
call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
select case (init_)
case('Z')
tx(:) = dzero
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(done,y,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,initu,dzero,tx,desc_data,info)
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,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
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=izero)
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,init='Y')
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.
do i=1, sweeps-1
!
! 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.
!
select case (init_)
case('Z')
ty = dzero
case('Y')
call psb_geaxpby(done,y,dzero,ty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,initu,dzero,ty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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,init='Y')
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
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
@ -509,17 +205,26 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
!
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
if (info ==0) deallocate(ww,tx,ty,stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

468
mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90
Unescape Escape View File

@ -54,12 +54,11 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
type(psb_d_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: 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
real(psb_dpk_), pointer :: aux(:)
type(psb_d_vect_type) :: tx, ty, ww
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='d_as_smoother_apply', ch_err
character(len=20) :: name='d_as_smoother_apply_v', ch_err
call psb_erractionsave(err_act)
@ -95,55 +94,21 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
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,izero,izero,izero,izero/),&
& 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,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
if (4*isz <= size(work)) then
aux => work(:)
else
allocate(ww(isz),tx(isz),ty(isz),&
&aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
endif
if (sweeps == 0) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -153,382 +118,115 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call vtx%zero()
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(done,y,dzero,vtx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,initu,dzero,vtx,desc_data,info)
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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=izero)
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,init='Y')
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
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
do i=1, sweeps-1
!
! 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.
!
! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system.
!
!
select case (init_)
case('Z')
call vty%zero()
case('Y')
call psb_geaxpby(done,y,dzero,vty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,initu,dzero,vty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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,init='Y')
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
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,vty,beta,y,desc_data,info)
end if
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
call vww%free(info)
call vtx%free(info)
call vty%free(info)
if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

150
mlprec/impl/smoother/mld_d_as_smoother_prol_a.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_d_as_smoother_prol_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_prol_a
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
real(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='d_as_smther_prol_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_as_smoother_prol_a

150
mlprec/impl/smoother/mld_d_as_smoother_prol_v.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_d_as_smoother_prol_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_prol_v
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='d_as_smther_prol_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_as_smoother_prol_v

169
mlprec/impl/smoother/mld_d_as_smoother_restr_a.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_d_as_smoother_restr_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_restr_a
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
real(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='d_as_smther_restr_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_as_smoother_restr_a

169
mlprec/impl/smoother/mld_d_as_smoother_restr_v.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_d_as_smoother_restr_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_restr_v
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='d_as_smther_restr_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_as_smoother_restr_v

457
mlprec/impl/smoother/mld_s_as_smoother_apply.f90
Unescape Escape View File

@ -54,15 +54,15 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
real(psb_spk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
real(psb_spk_), pointer :: ww(:), aux(:)
real(psb_spk_), allocatable :: tx(:),ty(:)
real(psb_spk_), pointer :: aux(:)
real(psb_spk_), allocatable :: tx(:),ty(:), ww(:)
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='s_as_smoother_apply', ch_err
call psb_erractionsave(err_act)
info = psb_success_
info = psb_success_
ictxt = desc_data%get_context()
call psb_info(ictxt,me,np)
@ -89,55 +89,26 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
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)
if ((6*isz) <= size(work)) then
ww => work(1:isz)
aux => work(3*isz+1:)
else if ((4*isz) <= size(work)) then
if ((4*isz) <= size(work)) then
aux => work(1:)
allocate(ww(isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/3*isz,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
else if ((3*isz) <= size(work)) then
ww => work(1:isz)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
else
allocate(ww(isz), aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
endif
if (sweeps == 0) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -147,361 +118,86 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info)
call psb_geasb(ty,sm%desc_data,info)
call psb_geasb(ww,sm%desc_data,info)
call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
select case (init_)
case('Z')
tx(:) = szero
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(sone,y,szero,tx,desc_data,info)
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,initu,szero,tx,desc_data,info)
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,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
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=izero)
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(sone,tx,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
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.
do i=1, sweeps-1
!
! 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.
!
select case (init_)
case('Z')
ty = szero
case('Y')
call psb_geaxpby(sone,y,szero,ty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,initu,szero,ty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
@ -509,17 +205,26 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
!
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
if (info ==0) deallocate(ww,tx,ty,stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

468
mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90
Unescape Escape View File

@ -54,12 +54,11 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
type(psb_s_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_spk_), allocatable :: vx(:)
type(psb_s_vect_type) :: vtx, vty, vww
real(psb_spk_), pointer :: aux(:)
type(psb_s_vect_type) :: tx, ty, ww
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='s_as_smoother_apply', ch_err
character(len=20) :: name='s_as_smoother_apply_v', ch_err
call psb_erractionsave(err_act)
@ -95,55 +94,21 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
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,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
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,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
if (4*isz <= size(work)) then
aux => work(:)
else
allocate(ww(isz),tx(isz),ty(isz),&
&aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
endif
if (sweeps == 0) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -153,382 +118,115 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call vtx%zero()
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(sone,y,szero,vtx,desc_data,info)
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,initu,szero,vtx,desc_data,info)
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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=izero)
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(sone,vtx,szero,vty,sm%desc_data,trans_,aux,info,init='Y')
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
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
do i=1, sweeps-1
!
! 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.
!
! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system.
!
!
select case (init_)
case('Z')
call vty%zero()
case('Y')
call psb_geaxpby(sone,y,szero,vty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,initu,szero,vty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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(sone,vtx,szero,vww,sm%desc_data,info)
call psb_spmm(-sone,sm%nd,vty,sone,vww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(sone,vww,szero,vty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,vty,beta,y,desc_data,info)
end if
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
call vww%free(info)
call vtx%free(info)
call vty%free(info)
if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

150
mlprec/impl/smoother/mld_s_as_smoother_prol_a.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_s_as_smoother_prol_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_prol_a
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
real(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='s_as_smther_prol_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_as_smoother_prol_a

150
mlprec/impl/smoother/mld_s_as_smoother_prol_v.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_s_as_smoother_prol_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_prol_v
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
type(psb_s_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='s_as_smther_prol_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_as_smoother_prol_v

169
mlprec/impl/smoother/mld_s_as_smoother_restr_a.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_s_as_smoother_restr_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_restr_a
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
real(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='s_as_smther_restr_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_as_smoother_restr_a

169
mlprec/impl/smoother/mld_s_as_smoother_restr_v.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_s_as_smoother_restr_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_restr_v
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
type(psb_s_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='s_as_smther_restr_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_as_smoother_restr_v

457
mlprec/impl/smoother/mld_z_as_smoother_apply.f90
Unescape Escape View File

@ -54,15 +54,15 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
complex(psb_dpk_),intent(inout), optional :: initu(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
complex(psb_dpk_), pointer :: ww(:), aux(:)
complex(psb_dpk_), allocatable :: tx(:),ty(:)
complex(psb_dpk_), pointer :: aux(:)
complex(psb_dpk_), allocatable :: tx(:),ty(:), ww(:)
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='z_as_smoother_apply', ch_err
call psb_erractionsave(err_act)
info = psb_success_
info = psb_success_
ictxt = desc_data%get_context()
call psb_info(ictxt,me,np)
@ -89,55 +89,26 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
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)
if ((6*isz) <= size(work)) then
ww => work(1:isz)
aux => work(3*isz+1:)
else if ((4*isz) <= size(work)) then
if ((4*isz) <= size(work)) then
aux => work(1:)
allocate(ww(isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/3*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
else if ((3*isz) <= size(work)) then
ww => work(1:isz)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
else
allocate(ww(isz), aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
endif
if (sweeps == 0) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -147,361 +118,86 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info)
call psb_geasb(ty,sm%desc_data,info)
call psb_geasb(ww,sm%desc_data,info)
call psb_geasb(tx,desc_data,info)
call psb_geasb(ty,desc_data,info)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
select case (init_)
case('Z')
tx(:) = zzero
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(zone,y,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,initu,zzero,tx,desc_data,info)
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,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
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=izero)
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(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
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.
do i=1, sweeps-1
!
! 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.
!
select case (init_)
case('Z')
ty = zzero
case('Y')
call psb_geaxpby(zone,y,zzero,ty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,initu,zzero,ty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
@ -509,17 +205,26 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
!
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
if (info ==0) deallocate(ww,tx,ty,stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

468
mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90
Unescape Escape View File

@ -54,12 +54,11 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
type(psb_z_vect_type),intent(inout), optional :: initu
integer(psb_ipk_) :: n_row,n_col, nrow_d, i
complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_dpk_), allocatable :: vx(:)
type(psb_z_vect_type) :: vtx, vty, vww
complex(psb_dpk_), pointer :: aux(:)
type(psb_z_vect_type) :: tx, ty, ww
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5)
character :: trans_, init_
character(len=20) :: name='z_as_smoother_apply', ch_err
character(len=20) :: name='z_as_smoother_apply_v', ch_err
call psb_erractionsave(err_act)
@ -95,55 +94,21 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
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,izero,izero,izero,izero/),&
& a_err='complex(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,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
if (4*isz <= size(work)) then
aux => work(:)
else
allocate(ww(isz),tx(isz),ty(isz),&
&aux(4*isz),stat=info)
allocate(aux(4*isz),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_alloc_request_,name,&
& i_err=(/4*isz,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
endif
if (sweeps == 0) then
!
! K^0 = I
! zero sweeps of any smoother is just the identity.
!
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then
if ((.not.sm%sv%is_iterative()).and.(sweeps == 1).and.(sm%novr==0)) then
!
! Shortcut: in this case it's just the same
! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
! Shortcut: in this case there is nothing else to be done.
!
call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info)
@ -153,382 +118,115 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
goto 9999
endif
else
else if (sweeps >= 0) then
!
!
! Apply multiple sweeps of an AS solver
! to compute an approximate solution of a linear system.
!
!
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vtx,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(vww,sm%desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call vtx%zero()
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(zone,y,zzero,vtx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,initu,zzero,vtx,desc_data,info)
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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=izero)
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(zone,vtx,zzero,vty,sm%desc_data,trans_,aux,info,init='Y')
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
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
do i=1, sweeps-1
!
! 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.
!
! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver
! to compute an approximate solution of a linear system.
!
!
select case (init_)
case('Z')
call vty%zero()
case('Y')
call psb_geaxpby(zone,y,zzero,vty,sm%desc_data,info)
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,initu,zzero,vty,sm%desc_data,info)
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
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=izero)
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(zone,vtx,zzero,vww,sm%desc_data,info)
call psb_spmm(-zone,sm%nd,vty,zone,vww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(zone,vww,zzero,vty,sm%desc_data,trans_,aux,info,init='Y')
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
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
else
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
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
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,vty,beta,y,desc_data,info)
end if
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
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)
info = psb_err_iarg_neg_
call psb_errpush(info,name,&
& i_err=(/itwo,sweeps,izero,izero,izero/))
goto 9999
endif
if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info)
endif
call vww%free(info)
call vtx%free(info)
call vty%free(info)
if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
call psb_errpush(info,name)
goto 9999
end if
call psb_erractionrestore(err_act)
return

150
mlprec/impl/smoother/mld_z_as_smoother_prol_a.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_z_as_smoother_prol_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_prol_a
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
complex(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='z_as_smther_prol_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_as_smoother_prol_a

150
mlprec/impl/smoother/mld_z_as_smoother_prol_v.f90
Unescape Escape View File

@ -0,0 +1,150 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_z_as_smoother_prol_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_prol_v
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
type(psb_z_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='z_as_smther_prol_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
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(x,sm%desc_data,info,&
! & update=sm%prol,work=work)
case(psb_sum_,psb_avg_)
!
! Update the overlap of x
!
call psb_ovrl(x,sm%desc_data,info,&
& update=sm%prol,work=work)
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(x,sm%desc_data,info,&
& update=psb_sum_,work=work)
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
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_as_smoother_prol_v

169
mlprec/impl/smoother/mld_z_as_smoother_restr_a.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_z_as_smoother_restr_a(sm,x,trans,work,info,data)
use psb_base_mod
use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_restr_a
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
complex(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='z_as_smther_restr_a', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_as_smoother_restr_a

169
mlprec/impl/smoother/mld_z_as_smoother_restr_v.f90
Unescape Escape View File

@ -0,0 +1,169 @@
!!$
!!$
!!$ MLD2P4 version 2.0
!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package
!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3)
!!$
!!$ (C) Copyright 2008, 2010, 2012, 2015
!!$
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$ Pasqua D'Ambra ICAR-CNR, Naples
!!$ Daniela di Serafino Second University of Naples
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the MLD2P4 group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine mld_z_as_smoother_restr_v(sm,x,trans,work,info,data)
use psb_base_mod
use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_restr_v
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
type(psb_z_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
!Local
integer(psb_ipk_) :: ictxt,np,me, err_act,isz,int_err(5), data_
character :: trans_
character(len=20) :: name='z_as_smther_restr_v', ch_err
call psb_erractionsave(err_act)
info = psb_success_
ictxt = sm%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
info = psb_err_iarg_invalid_i_
call psb_errpush(info,name)
goto 9999
end select
if (present(data)) then
data_ = data
else
data_ = psb_comm_ext_
end if
if (.not.allocated(sm%sv)) then
info = 1121
call psb_errpush(info,name)
goto 9999
end if
select case(trans_)
case('N')
!
! Get the overlap entries x
!
if (sm%restr == psb_halo_) then
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,work=work,data=data_)
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(x,sm%desc_data,info,&
& update=psb_avg_,work=work,mode=izero)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_ovrl'
goto 9999
end if
call psb_halo(x,sm%desc_data,info,work=work,data=data_)
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 psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_as_smoother_restr_v

67
mlprec/mld_c_as_smoother.f90
Unescape Escape View File

@ -62,6 +62,12 @@ module mld_c_as_smoother
procedure, pass(sm) :: clone => mld_c_as_smoother_clone
procedure, pass(sm) :: apply_v => mld_c_as_smoother_apply_vect
procedure, pass(sm) :: apply_a => mld_c_as_smoother_apply
procedure, pass(sm) :: restr_a => mld_c_as_smoother_restr_a
procedure, pass(sm) :: prol_a => mld_c_as_smoother_prol_a
procedure, pass(sm) :: restr_v => mld_c_as_smoother_restr_v
procedure, pass(sm) :: prol_v => mld_c_as_smoother_prol_v
generic, public :: apply_restr => restr_v, restr_a
generic, public :: apply_prol => prol_v, prol_a
procedure, pass(sm) :: free => mld_c_as_smoother_free
procedure, pass(sm) :: seti => mld_c_as_smoother_seti
procedure, pass(sm) :: setc => mld_c_as_smoother_setc
@ -95,6 +101,67 @@ module mld_c_as_smoother
end subroutine mld_c_as_smoother_check
end interface
interface
subroutine mld_c_as_smoother_restr_v(sm,x,trans,work,info,data)
import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
& psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
type(psb_c_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_c_as_smoother_restr_v
end interface
interface
subroutine mld_c_as_smoother_restr_a(sm,x,trans,work,info,data)
import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
& psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
complex(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_c_as_smoother_restr_a
end interface
interface
subroutine mld_c_as_smoother_prol_v(sm,x,trans,work,info,data)
import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
& psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
type(psb_c_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_c_as_smoother_prol_v
end interface
interface
subroutine mld_c_as_smoother_prol_a(sm,x,trans,work,info,data)
import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
& psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_c_as_smoother_type), intent(inout) :: sm
complex(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_c_as_smoother_prol_a
end interface
interface
subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu)

67
mlprec/mld_d_as_smoother.f90
Unescape Escape View File

@ -62,6 +62,12 @@ module mld_d_as_smoother
procedure, pass(sm) :: clone => mld_d_as_smoother_clone
procedure, pass(sm) :: apply_v => mld_d_as_smoother_apply_vect
procedure, pass(sm) :: apply_a => mld_d_as_smoother_apply
procedure, pass(sm) :: restr_a => mld_d_as_smoother_restr_a
procedure, pass(sm) :: prol_a => mld_d_as_smoother_prol_a
procedure, pass(sm) :: restr_v => mld_d_as_smoother_restr_v
procedure, pass(sm) :: prol_v => mld_d_as_smoother_prol_v
generic, public :: apply_restr => restr_v, restr_a
generic, public :: apply_prol => prol_v, prol_a
procedure, pass(sm) :: free => mld_d_as_smoother_free
procedure, pass(sm) :: seti => mld_d_as_smoother_seti
procedure, pass(sm) :: setc => mld_d_as_smoother_setc
@ -95,6 +101,67 @@ module mld_d_as_smoother
end subroutine mld_d_as_smoother_check
end interface
interface
subroutine mld_d_as_smoother_restr_v(sm,x,trans,work,info,data)
import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
& psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_d_as_smoother_restr_v
end interface
interface
subroutine mld_d_as_smoother_restr_a(sm,x,trans,work,info,data)
import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
& psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
real(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_d_as_smoother_restr_a
end interface
interface
subroutine mld_d_as_smoother_prol_v(sm,x,trans,work,info,data)
import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
& psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
type(psb_d_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_d_as_smoother_prol_v
end interface
interface
subroutine mld_d_as_smoother_prol_a(sm,x,trans,work,info,data)
import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
& psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_d_as_smoother_type), intent(inout) :: sm
real(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_d_as_smoother_prol_a
end interface
interface
subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu)

67
mlprec/mld_s_as_smoother.f90
Unescape Escape View File

@ -62,6 +62,12 @@ module mld_s_as_smoother
procedure, pass(sm) :: clone => mld_s_as_smoother_clone
procedure, pass(sm) :: apply_v => mld_s_as_smoother_apply_vect
procedure, pass(sm) :: apply_a => mld_s_as_smoother_apply
procedure, pass(sm) :: restr_a => mld_s_as_smoother_restr_a
procedure, pass(sm) :: prol_a => mld_s_as_smoother_prol_a
procedure, pass(sm) :: restr_v => mld_s_as_smoother_restr_v
procedure, pass(sm) :: prol_v => mld_s_as_smoother_prol_v
generic, public :: apply_restr => restr_v, restr_a
generic, public :: apply_prol => prol_v, prol_a
procedure, pass(sm) :: free => mld_s_as_smoother_free
procedure, pass(sm) :: seti => mld_s_as_smoother_seti
procedure, pass(sm) :: setc => mld_s_as_smoother_setc
@ -95,6 +101,67 @@ module mld_s_as_smoother
end subroutine mld_s_as_smoother_check
end interface
interface
subroutine mld_s_as_smoother_restr_v(sm,x,trans,work,info,data)
import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, &
& psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
type(psb_s_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_s_as_smoother_restr_v
end interface
interface
subroutine mld_s_as_smoother_restr_a(sm,x,trans,work,info,data)
import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, &
& psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
real(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_s_as_smoother_restr_a
end interface
interface
subroutine mld_s_as_smoother_prol_v(sm,x,trans,work,info,data)
import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, &
& psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
type(psb_s_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_s_as_smoother_prol_v
end interface
interface
subroutine mld_s_as_smoother_prol_a(sm,x,trans,work,info,data)
import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, &
& psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_s_as_smoother_type), intent(inout) :: sm
real(psb_spk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
real(psb_spk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_s_as_smoother_prol_a
end interface
interface
subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu)

67
mlprec/mld_z_as_smoother.f90
Unescape Escape View File

@ -62,6 +62,12 @@ module mld_z_as_smoother
procedure, pass(sm) :: clone => mld_z_as_smoother_clone
procedure, pass(sm) :: apply_v => mld_z_as_smoother_apply_vect
procedure, pass(sm) :: apply_a => mld_z_as_smoother_apply
procedure, pass(sm) :: restr_a => mld_z_as_smoother_restr_a
procedure, pass(sm) :: prol_a => mld_z_as_smoother_prol_a
procedure, pass(sm) :: restr_v => mld_z_as_smoother_restr_v
procedure, pass(sm) :: prol_v => mld_z_as_smoother_prol_v
generic, public :: apply_restr => restr_v, restr_a
generic, public :: apply_prol => prol_v, prol_a
procedure, pass(sm) :: free => mld_z_as_smoother_free
procedure, pass(sm) :: seti => mld_z_as_smoother_seti
procedure, pass(sm) :: setc => mld_z_as_smoother_setc
@ -95,6 +101,67 @@ module mld_z_as_smoother
end subroutine mld_z_as_smoother_check
end interface
interface
subroutine mld_z_as_smoother_restr_v(sm,x,trans,work,info,data)
import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &
& psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
type(psb_z_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_z_as_smoother_restr_v
end interface
interface
subroutine mld_z_as_smoother_restr_a(sm,x,trans,work,info,data)
import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &
& psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
complex(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_z_as_smoother_restr_a
end interface
interface
subroutine mld_z_as_smoother_prol_v(sm,x,trans,work,info,data)
import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &
& psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
type(psb_z_vect_type),intent(inout) :: x
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_z_as_smoother_prol_v
end interface
interface
subroutine mld_z_as_smoother_prol_a(sm,x,trans,work,info,data)
import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &
& psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_
implicit none
class(mld_z_as_smoother_type), intent(inout) :: sm
complex(psb_dpk_), intent(inout) :: x(:)
character(len=1),intent(in) :: trans
complex(psb_dpk_),target, intent(inout) :: work(:)
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: data
end subroutine mld_z_as_smoother_prol_a
end interface
interface
subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu)

Write Preview
Loading…
Cancel
Save