mirror of
https://github.com/sfilippone/amg4psblas.git
synced 2026-10-06 22:55:12 +00:00
Fixes for compilation with INTEL
This commit is contained in:
sfilippone
@@ -80,7 +80,7 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
& a,desc_a,ilaggr,nlaggr,op_prol,info)
|
||||
use psb_base_mod
|
||||
use amg_c_prec_type
|
||||
use amg_c_symdec_aggregator_mod, amg_protect_name => amg_c_symdec_aggregator_build_tprol
|
||||
use amg_c_symdec_aggregator_mod, only : amg_c_symdec_aggregator_type
|
||||
use amg_c_inner_mod
|
||||
implicit none
|
||||
class(amg_c_symdec_aggregator_type), target, intent(inout) :: ag
|
||||
@@ -103,6 +103,16 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
integer(psb_ipk_) :: debug_level, debug_unit
|
||||
logical :: clean_zeros
|
||||
|
||||
!!$ interface
|
||||
!!$ subroutine psb_caplusat(ain,aout,info)
|
||||
!!$ use psb_c_mat_mod, only : psb_cspmat_type
|
||||
!!$ import :: psb_ipk_
|
||||
!!$ implicit none
|
||||
!!$ type(psb_cspmat_type) :: ain, aout
|
||||
!!$ integer(psb_ipk_) :: info
|
||||
!!$ end subroutine psb_caplusat
|
||||
!!$ end interface
|
||||
!!$
|
||||
name='amg_c_symdec_aggregator_tprol'
|
||||
call psb_erractionsave(err_act)
|
||||
if (psb_errstatus_fatal()) then
|
||||
@@ -123,15 +133,11 @@ subroutine amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs)
|
||||
|
||||
nr = a%get_nrows()
|
||||
call a%csclip(atmp,info,imax=nr,jmax=nr,&
|
||||
call a%csclip(atrans,info,imax=nr,jmax=nr,&
|
||||
& rscale=.false.,cscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
if (info == psb_success_) call atmp%transp(atrans)
|
||||
if (info == psb_success_) call atrans%cscnv(info,type='COO')
|
||||
if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
call atrans%set_nrows(nr)
|
||||
call atrans%set_ncols(nr)
|
||||
call psb_aplusat(atrans,atmp,info)
|
||||
if (info == psb_success_) call atrans%free()
|
||||
if (info == psb_success_) call atmp%cscnv(info,type='CSR')
|
||||
|
||||
|
||||
@@ -80,7 +80,7 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
& a,desc_a,ilaggr,nlaggr,op_prol,info)
|
||||
use psb_base_mod
|
||||
use amg_d_prec_type
|
||||
use amg_d_symdec_aggregator_mod, amg_protect_name => amg_d_symdec_aggregator_build_tprol
|
||||
use amg_d_symdec_aggregator_mod, only : amg_d_symdec_aggregator_type
|
||||
use amg_d_inner_mod
|
||||
implicit none
|
||||
class(amg_d_symdec_aggregator_type), target, intent(inout) :: ag
|
||||
@@ -103,6 +103,16 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
integer(psb_ipk_) :: debug_level, debug_unit
|
||||
logical :: clean_zeros
|
||||
|
||||
!!$ interface
|
||||
!!$ subroutine psb_daplusat(ain,aout,info)
|
||||
!!$ use psb_d_mat_mod, only : psb_dspmat_type
|
||||
!!$ import :: psb_ipk_
|
||||
!!$ implicit none
|
||||
!!$ type(psb_dspmat_type) :: ain, aout
|
||||
!!$ integer(psb_ipk_) :: info
|
||||
!!$ end subroutine psb_daplusat
|
||||
!!$ end interface
|
||||
!!$
|
||||
name='amg_d_symdec_aggregator_tprol'
|
||||
call psb_erractionsave(err_act)
|
||||
if (psb_errstatus_fatal()) then
|
||||
@@ -123,15 +133,11 @@ subroutine amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs)
|
||||
|
||||
nr = a%get_nrows()
|
||||
call a%csclip(atmp,info,imax=nr,jmax=nr,&
|
||||
call a%csclip(atrans,info,imax=nr,jmax=nr,&
|
||||
& rscale=.false.,cscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
if (info == psb_success_) call atmp%transp(atrans)
|
||||
if (info == psb_success_) call atrans%cscnv(info,type='COO')
|
||||
if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
call atrans%set_nrows(nr)
|
||||
call atrans%set_ncols(nr)
|
||||
call psb_aplusat(atrans,atmp,info)
|
||||
if (info == psb_success_) call atrans%free()
|
||||
if (info == psb_success_) call atmp%cscnv(info,type='CSR')
|
||||
|
||||
|
||||
@@ -80,7 +80,7 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
& a,desc_a,ilaggr,nlaggr,op_prol,info)
|
||||
use psb_base_mod
|
||||
use amg_s_prec_type
|
||||
use amg_s_symdec_aggregator_mod, amg_protect_name => amg_s_symdec_aggregator_build_tprol
|
||||
use amg_s_symdec_aggregator_mod, only : amg_s_symdec_aggregator_type
|
||||
use amg_s_inner_mod
|
||||
implicit none
|
||||
class(amg_s_symdec_aggregator_type), target, intent(inout) :: ag
|
||||
@@ -103,6 +103,16 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
integer(psb_ipk_) :: debug_level, debug_unit
|
||||
logical :: clean_zeros
|
||||
|
||||
!!$ interface
|
||||
!!$ subroutine psb_saplusat(ain,aout,info)
|
||||
!!$ use psb_s_mat_mod, only : psb_sspmat_type
|
||||
!!$ import :: psb_ipk_
|
||||
!!$ implicit none
|
||||
!!$ type(psb_sspmat_type) :: ain, aout
|
||||
!!$ integer(psb_ipk_) :: info
|
||||
!!$ end subroutine psb_saplusat
|
||||
!!$ end interface
|
||||
!!$
|
||||
name='amg_s_symdec_aggregator_tprol'
|
||||
call psb_erractionsave(err_act)
|
||||
if (psb_errstatus_fatal()) then
|
||||
@@ -123,15 +133,11 @@ subroutine amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs)
|
||||
|
||||
nr = a%get_nrows()
|
||||
call a%csclip(atmp,info,imax=nr,jmax=nr,&
|
||||
call a%csclip(atrans,info,imax=nr,jmax=nr,&
|
||||
& rscale=.false.,cscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
if (info == psb_success_) call atmp%transp(atrans)
|
||||
if (info == psb_success_) call atrans%cscnv(info,type='COO')
|
||||
if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
call atrans%set_nrows(nr)
|
||||
call atrans%set_ncols(nr)
|
||||
call psb_aplusat(atrans,atmp,info)
|
||||
if (info == psb_success_) call atrans%free()
|
||||
if (info == psb_success_) call atmp%cscnv(info,type='CSR')
|
||||
|
||||
|
||||
@@ -80,7 +80,7 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
& a,desc_a,ilaggr,nlaggr,op_prol,info)
|
||||
use psb_base_mod
|
||||
use amg_z_prec_type
|
||||
use amg_z_symdec_aggregator_mod, amg_protect_name => amg_z_symdec_aggregator_build_tprol
|
||||
use amg_z_symdec_aggregator_mod, only : amg_z_symdec_aggregator_type
|
||||
use amg_z_inner_mod
|
||||
implicit none
|
||||
class(amg_z_symdec_aggregator_type), target, intent(inout) :: ag
|
||||
@@ -103,6 +103,16 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
integer(psb_ipk_) :: debug_level, debug_unit
|
||||
logical :: clean_zeros
|
||||
|
||||
!!$ interface
|
||||
!!$ subroutine psb_zaplusat(ain,aout,info)
|
||||
!!$ use psb_z_mat_mod, only : psb_zspmat_type
|
||||
!!$ import :: psb_ipk_
|
||||
!!$ implicit none
|
||||
!!$ type(psb_zspmat_type) :: ain, aout
|
||||
!!$ integer(psb_ipk_) :: info
|
||||
!!$ end subroutine psb_zaplusat
|
||||
!!$ end interface
|
||||
!!$
|
||||
name='amg_z_symdec_aggregator_tprol'
|
||||
call psb_erractionsave(err_act)
|
||||
if (psb_errstatus_fatal()) then
|
||||
@@ -123,15 +133,11 @@ subroutine amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,&
|
||||
call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs)
|
||||
|
||||
nr = a%get_nrows()
|
||||
call a%csclip(atmp,info,imax=nr,jmax=nr,&
|
||||
call a%csclip(atrans,info,imax=nr,jmax=nr,&
|
||||
& rscale=.false.,cscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
if (info == psb_success_) call atmp%transp(atrans)
|
||||
if (info == psb_success_) call atrans%cscnv(info,type='COO')
|
||||
if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.)
|
||||
call atmp%set_nrows(nr)
|
||||
call atmp%set_ncols(nr)
|
||||
call atrans%set_nrows(nr)
|
||||
call atrans%set_ncols(nr)
|
||||
call psb_aplusat(atrans,atmp,info)
|
||||
if (info == psb_success_) call atrans%free()
|
||||
if (info == psb_success_) call atmp%cscnv(info,type='CSR')
|
||||
|
||||
|
||||
@@ -81,8 +81,7 @@
|
||||
subroutine amg_c_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
!use amg_c_inner_mod
|
||||
use amg_c_prec_mod, amg_protect_name => amg_c_smoothers_bld
|
||||
use amg_c_prec_type
|
||||
|
||||
Implicit None
|
||||
|
||||
|
||||
@@ -81,8 +81,7 @@
|
||||
subroutine amg_d_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
!use amg_d_inner_mod
|
||||
use amg_d_prec_mod, amg_protect_name => amg_d_smoothers_bld
|
||||
use amg_d_prec_type
|
||||
|
||||
Implicit None
|
||||
|
||||
|
||||
@@ -81,8 +81,7 @@
|
||||
subroutine amg_s_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
!use amg_s_inner_mod
|
||||
use amg_s_prec_mod, amg_protect_name => amg_s_smoothers_bld
|
||||
use amg_s_prec_type
|
||||
|
||||
Implicit None
|
||||
|
||||
|
||||
@@ -81,8 +81,7 @@
|
||||
subroutine amg_z_smoothers_bld(a,desc_a,prec,info,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
!use amg_z_inner_mod
|
||||
use amg_z_prec_mod, amg_protect_name => amg_z_smoothers_bld
|
||||
use amg_z_prec_type
|
||||
|
||||
Implicit None
|
||||
|
||||
|
||||
@@ -54,6 +54,7 @@ amg_c_mumps_solver_apply.o \
|
||||
amg_c_mumps_solver_apply_vect.o \
|
||||
amg_c_mumps_solver_bld.o \
|
||||
amg_c_krm_solver_impl.o \
|
||||
amg_c_slu_solver_impl.o \
|
||||
amg_d_base_solver_apply.o \
|
||||
amg_d_base_solver_apply_vect.o \
|
||||
amg_d_base_solver_bld.o \
|
||||
@@ -101,6 +102,8 @@ amg_d_mumps_solver_apply.o \
|
||||
amg_d_mumps_solver_apply_vect.o \
|
||||
amg_d_mumps_solver_bld.o \
|
||||
amg_d_krm_solver_impl.o \
|
||||
amg_d_slu_solver_impl.o \
|
||||
amg_d_umf_solver_impl.o \
|
||||
amg_s_base_solver_apply.o \
|
||||
amg_s_base_solver_apply_vect.o \
|
||||
amg_s_base_solver_bld.o \
|
||||
@@ -148,6 +151,7 @@ amg_s_mumps_solver_apply.o \
|
||||
amg_s_mumps_solver_apply_vect.o \
|
||||
amg_s_mumps_solver_bld.o \
|
||||
amg_s_krm_solver_impl.o \
|
||||
amg_s_slu_solver_impl.o \
|
||||
amg_z_base_solver_apply.o \
|
||||
amg_z_base_solver_apply_vect.o \
|
||||
amg_z_base_solver_bld.o \
|
||||
@@ -303,6 +307,8 @@ amg_z_invk_solver_clone_settings.o \
|
||||
amg_z_invk_solver_cseti.o \
|
||||
amg_z_invk_solver_descr.o \
|
||||
amg_z_krm_solver_impl.o \
|
||||
amg_z_slu_solver_impl.o \
|
||||
amg_z_umf_solver_impl.o \
|
||||
amg_c_jac_solver_clone_settings.o \
|
||||
amg_c_jac_solver_clear_data.o \
|
||||
amg_c_jac_solver_cnv.o \
|
||||
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_c_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
use amg_c_slu_solver, amg_protect_name => amg_c_slu_solver_bld
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_cspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_c_slu_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_cspmat_type), intent(in), target, optional :: b
|
||||
class(psb_c_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_c_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_cspmat_type) :: atmp
|
||||
type(psb_c_csc_sparse_mat) :: acsc
|
||||
type(psb_c_coo_sparse_mat) :: acoo
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_slu_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_)
|
||||
nrow_a = atmp%get_nrows()
|
||||
call atmp%a%csclip(acoo,info,jmax=nrow_a)
|
||||
call acsc%mv_from_coo(acoo,info)
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entries to call C-base SuperLU
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_cslu_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%icp,acsc%ia,sv%lufactors)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_cslu_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_c_slu_solver_bld
|
||||
|
||||
subroutine amg_c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_c_slu_solver, amg_protect_name => amg_c_slu_solver_apply_vect
|
||||
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_c_slu_solver_type), intent(inout) :: sv
|
||||
type(psb_c_vect_type),intent(inout) :: x
|
||||
type(psb_c_vect_type),intent(inout) :: y
|
||||
complex(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_c_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_c_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='s_slu_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_c_slu_solver_apply_vect
|
||||
|
||||
subroutine amg_c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_c_slu_solver, amg_protect_name => amg_c_slu_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_c_slu_solver_type), intent(inout) :: sv
|
||||
complex(psb_spk_),intent(inout) :: x(:)
|
||||
complex(psb_spk_),intent(inout) :: y(:)
|
||||
complex(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_spk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
complex(psb_spk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
complex(psb_spk_), pointer :: ww(:)
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='s_slu_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='complex(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
ww(1:n_row) = x(1:n_row)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_cslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('T')
|
||||
info = amg_cslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('C')
|
||||
info = amg_cslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_, &
|
||||
& name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) &
|
||||
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_c_slu_solver_apply
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_c_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_c_umf_solver, amg_protect_name => amg_c_umf_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_c_umf_solver_type), intent(inout) :: sv
|
||||
complex(psb_spk_),intent(inout) :: x(:)
|
||||
complex(psb_spk_),intent(inout) :: y(:)
|
||||
complex(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_spk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
complex(psb_spk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
complex(psb_spk_), pointer :: ww(:)
|
||||
integer(psb_ipk_) :: i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='c_umf_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='complex(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_cumf_solve(0,n_row,ww,x,n_row,sv%numeric)
|
||||
case('T')
|
||||
!
|
||||
! Note: with UMF, 1 meand Ctranspose, 2 means transpose
|
||||
! even for complex data.
|
||||
!
|
||||
if (psb_c_is_complex_) then
|
||||
info = amg_cumf_solve(2,n_row,ww,x,n_row,sv%numeric)
|
||||
else
|
||||
info = amg_cumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
end if
|
||||
case('C')
|
||||
info = amg_cumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_c_umf_solver_apply
|
||||
|
||||
subroutine amg_c_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_c_umf_solver, amg_protect_name => amg_c_umf_solver_apply_vect
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_c_umf_solver_type), intent(inout) :: sv
|
||||
type(psb_c_vect_type),intent(inout) :: x
|
||||
type(psb_c_vect_type),intent(inout) :: y
|
||||
complex(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_c_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_c_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='c_umf_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_c_umf_solver_apply_vect
|
||||
|
||||
subroutine amg_c_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
use psb_base_mod
|
||||
use amg_c_umf_solver, amg_protect_name => amg_c_umf_solver_bld
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_cspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_c_umf_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_cspmat_type), intent(in), target, optional :: b
|
||||
class(psb_c_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_c_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_cspmat_type) :: atmp
|
||||
type(psb_c_csc_sparse_mat) :: acsc
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='c_umf_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
|
||||
call atmp%mv_to(acsc)
|
||||
nrow_a = acsc%get_nrows()
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entres to call C-base UMFPACK.
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_cumf_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%ia,acsc%icp,sv%symbolic,sv%numeric,&
|
||||
& sv%symbsize,sv%numsize)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_cumf_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_c_umf_solver_bld
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_d_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
use amg_d_slu_solver, amg_protect_name => amg_d_slu_solver_bld
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_dspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_d_slu_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_dspmat_type), intent(in), target, optional :: b
|
||||
class(psb_d_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_d_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_dspmat_type) :: atmp
|
||||
type(psb_d_csc_sparse_mat) :: acsc
|
||||
type(psb_d_coo_sparse_mat) :: acoo
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_slu_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_)
|
||||
nrow_a = atmp%get_nrows()
|
||||
call atmp%a%csclip(acoo,info,jmax=nrow_a)
|
||||
call acsc%mv_from_coo(acoo,info)
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entries to call C-base SuperLU
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_dslu_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%icp,acsc%ia,sv%lufactors)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_dslu_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_d_slu_solver_bld
|
||||
|
||||
subroutine amg_d_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_d_slu_solver, amg_protect_name => amg_d_slu_solver_apply_vect
|
||||
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_d_slu_solver_type), intent(inout) :: sv
|
||||
type(psb_d_vect_type),intent(inout) :: x
|
||||
type(psb_d_vect_type),intent(inout) :: y
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_d_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_d_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='s_slu_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_d_slu_solver_apply_vect
|
||||
|
||||
subroutine amg_d_slu_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_d_slu_solver, amg_protect_name => amg_d_slu_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_d_slu_solver_type), intent(inout) :: sv
|
||||
real(psb_dpk_),intent(inout) :: x(:)
|
||||
real(psb_dpk_),intent(inout) :: y(:)
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
real(psb_dpk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
real(psb_dpk_), pointer :: ww(:)
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='s_slu_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
ww(1:n_row) = x(1:n_row)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_dslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('T')
|
||||
info = amg_dslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('C')
|
||||
info = amg_dslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_, &
|
||||
& name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) &
|
||||
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_d_slu_solver_apply
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_d_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_d_umf_solver, amg_protect_name => amg_d_umf_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_d_umf_solver_type), intent(inout) :: sv
|
||||
real(psb_dpk_),intent(inout) :: x(:)
|
||||
real(psb_dpk_),intent(inout) :: y(:)
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
real(psb_dpk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
real(psb_dpk_), pointer :: ww(:)
|
||||
integer(psb_ipk_) :: i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='d_umf_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_dumf_solve(0,n_row,ww,x,n_row,sv%numeric)
|
||||
case('T')
|
||||
!
|
||||
! Note: with UMF, 1 meand Ctranspose, 2 means transpose
|
||||
! even for complex data.
|
||||
!
|
||||
if (psb_d_is_complex_) then
|
||||
info = amg_dumf_solve(2,n_row,ww,x,n_row,sv%numeric)
|
||||
else
|
||||
info = amg_dumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
end if
|
||||
case('C')
|
||||
info = amg_dumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_d_umf_solver_apply
|
||||
|
||||
subroutine amg_d_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_d_umf_solver, amg_protect_name => amg_d_umf_solver_apply_vect
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_d_umf_solver_type), intent(inout) :: sv
|
||||
type(psb_d_vect_type),intent(inout) :: x
|
||||
type(psb_d_vect_type),intent(inout) :: y
|
||||
real(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_d_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_d_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='d_umf_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_d_umf_solver_apply_vect
|
||||
|
||||
subroutine amg_d_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
use psb_base_mod
|
||||
use amg_d_umf_solver, amg_protect_name => amg_d_umf_solver_bld
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_dspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_d_umf_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_dspmat_type), intent(in), target, optional :: b
|
||||
class(psb_d_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_d_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_dspmat_type) :: atmp
|
||||
type(psb_d_csc_sparse_mat) :: acsc
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='d_umf_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
|
||||
call atmp%mv_to(acsc)
|
||||
nrow_a = acsc%get_nrows()
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entres to call C-base UMFPACK.
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_dumf_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%ia,acsc%icp,sv%symbolic,sv%numeric,&
|
||||
& sv%symbsize,sv%numsize)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_dumf_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_d_umf_solver_bld
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_s_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
use amg_s_slu_solver, amg_protect_name => amg_s_slu_solver_bld
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_sspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_s_slu_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_sspmat_type), intent(in), target, optional :: b
|
||||
class(psb_s_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_s_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_sspmat_type) :: atmp
|
||||
type(psb_s_csc_sparse_mat) :: acsc
|
||||
type(psb_s_coo_sparse_mat) :: acoo
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_slu_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_)
|
||||
nrow_a = atmp%get_nrows()
|
||||
call atmp%a%csclip(acoo,info,jmax=nrow_a)
|
||||
call acsc%mv_from_coo(acoo,info)
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entries to call C-base SuperLU
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_sslu_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%icp,acsc%ia,sv%lufactors)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_sslu_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_s_slu_solver_bld
|
||||
|
||||
subroutine amg_s_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_s_slu_solver, amg_protect_name => amg_s_slu_solver_apply_vect
|
||||
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_s_slu_solver_type), intent(inout) :: sv
|
||||
type(psb_s_vect_type),intent(inout) :: x
|
||||
type(psb_s_vect_type),intent(inout) :: y
|
||||
real(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_s_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_s_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='s_slu_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_s_slu_solver_apply_vect
|
||||
|
||||
subroutine amg_s_slu_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_s_slu_solver, amg_protect_name => amg_s_slu_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_s_slu_solver_type), intent(inout) :: sv
|
||||
real(psb_spk_),intent(inout) :: x(:)
|
||||
real(psb_spk_),intent(inout) :: y(:)
|
||||
real(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_spk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
real(psb_spk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
real(psb_spk_), pointer :: ww(:)
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='s_slu_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
ww(1:n_row) = x(1:n_row)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_sslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('T')
|
||||
info = amg_sslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('C')
|
||||
info = amg_sslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_, &
|
||||
& name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) &
|
||||
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_s_slu_solver_apply
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_s_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_s_umf_solver, amg_protect_name => amg_s_umf_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_s_umf_solver_type), intent(inout) :: sv
|
||||
real(psb_spk_),intent(inout) :: x(:)
|
||||
real(psb_spk_),intent(inout) :: y(:)
|
||||
real(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_spk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
real(psb_spk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
real(psb_spk_), pointer :: ww(:)
|
||||
integer(psb_ipk_) :: i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='s_umf_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='real(psb_spk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_sumf_solve(0,n_row,ww,x,n_row,sv%numeric)
|
||||
case('T')
|
||||
!
|
||||
! Note: with UMF, 1 meand Ctranspose, 2 means transpose
|
||||
! even for complex data.
|
||||
!
|
||||
if (psb_s_is_complex_) then
|
||||
info = amg_sumf_solve(2,n_row,ww,x,n_row,sv%numeric)
|
||||
else
|
||||
info = amg_sumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
end if
|
||||
case('C')
|
||||
info = amg_sumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_s_umf_solver_apply
|
||||
|
||||
subroutine amg_s_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_s_umf_solver, amg_protect_name => amg_s_umf_solver_apply_vect
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_s_umf_solver_type), intent(inout) :: sv
|
||||
type(psb_s_vect_type),intent(inout) :: x
|
||||
type(psb_s_vect_type),intent(inout) :: y
|
||||
real(psb_spk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
real(psb_spk_),target, intent(inout) :: work(:)
|
||||
type(psb_s_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_s_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='s_umf_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_s_umf_solver_apply_vect
|
||||
|
||||
subroutine amg_s_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
use psb_base_mod
|
||||
use amg_s_umf_solver, amg_protect_name => amg_s_umf_solver_bld
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_sspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_s_umf_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_sspmat_type), intent(in), target, optional :: b
|
||||
class(psb_s_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_s_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_sspmat_type) :: atmp
|
||||
type(psb_s_csc_sparse_mat) :: acsc
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_umf_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
|
||||
call atmp%mv_to(acsc)
|
||||
nrow_a = acsc%get_nrows()
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entres to call C-base UMFPACK.
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_sumf_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%ia,acsc%icp,sv%symbolic,sv%numeric,&
|
||||
& sv%symbsize,sv%numsize)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_sumf_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_s_umf_solver_bld
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_z_slu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
|
||||
use psb_base_mod
|
||||
use amg_z_slu_solver, amg_protect_name => amg_z_slu_solver_bld
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_zspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_z_slu_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_zspmat_type), intent(in), target, optional :: b
|
||||
class(psb_z_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_z_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_zspmat_type) :: atmp
|
||||
type(psb_z_csc_sparse_mat) :: acsc
|
||||
type(psb_z_coo_sparse_mat) :: acoo
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='s_slu_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='coo',dupl=psb_dupl_add_)
|
||||
nrow_a = atmp%get_nrows()
|
||||
call atmp%a%csclip(acoo,info,jmax=nrow_a)
|
||||
call acsc%mv_from_coo(acoo,info)
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entries to call C-base SuperLU
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_zslu_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%icp,acsc%ia,sv%lufactors)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_zslu_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_z_slu_solver_bld
|
||||
|
||||
subroutine amg_z_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_z_slu_solver, amg_protect_name => amg_z_slu_solver_apply_vect
|
||||
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_z_slu_solver_type), intent(inout) :: sv
|
||||
type(psb_z_vect_type),intent(inout) :: x
|
||||
type(psb_z_vect_type),intent(inout) :: y
|
||||
complex(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_z_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_z_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='s_slu_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_z_slu_solver_apply_vect
|
||||
|
||||
subroutine amg_z_slu_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_z_slu_solver, amg_protect_name => amg_z_slu_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_z_slu_solver_type), intent(inout) :: sv
|
||||
complex(psb_dpk_),intent(inout) :: x(:)
|
||||
complex(psb_dpk_),intent(inout) :: y(:)
|
||||
complex(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_dpk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
complex(psb_dpk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
complex(psb_dpk_), pointer :: ww(:)
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='s_slu_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='complex(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
ww(1:n_row) = x(1:n_row)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_zslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('T')
|
||||
info = amg_zslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
|
||||
case('C')
|
||||
info = amg_zslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_, &
|
||||
& name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) &
|
||||
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,&
|
||||
& name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_z_slu_solver_apply
|
||||
@@ -0,0 +1,202 @@
|
||||
subroutine amg_z_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_z_umf_solver, amg_protect_name => amg_z_umf_solver_apply
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_z_umf_solver_type), intent(inout) :: sv
|
||||
complex(psb_dpk_),intent(inout) :: x(:)
|
||||
complex(psb_dpk_),intent(inout) :: y(:)
|
||||
complex(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_dpk_),target, intent(inout) :: work(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
complex(psb_dpk_),intent(inout), optional :: initu(:)
|
||||
|
||||
integer :: n_row,n_col
|
||||
complex(psb_dpk_), pointer :: ww(:)
|
||||
integer(psb_ipk_) :: i, err_act
|
||||
character :: trans_
|
||||
character(len=20) :: name='z_umf_solver_apply'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
|
||||
trans_ = psb_toupper(trans)
|
||||
select case(trans_)
|
||||
case('N')
|
||||
case('T','C')
|
||||
case default
|
||||
call psb_errpush(psb_err_iarg_invalid_i_,name)
|
||||
goto 9999
|
||||
end select
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
n_row = desc_data%get_local_rows()
|
||||
n_col = desc_data%get_local_cols()
|
||||
|
||||
if (n_col <= size(work)) then
|
||||
ww => work(1:n_col)
|
||||
else
|
||||
allocate(ww(n_col),stat=info)
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_alloc_request_
|
||||
call psb_errpush(info,name,i_err=(/n_col,0,0,0,0/),&
|
||||
& a_err='complex(psb_dpk_)')
|
||||
goto 9999
|
||||
end if
|
||||
endif
|
||||
|
||||
select case(trans_)
|
||||
case('N')
|
||||
info = amg_zumf_solve(0,n_row,ww,x,n_row,sv%numeric)
|
||||
case('T')
|
||||
!
|
||||
! Note: with UMF, 1 meand Ctranspose, 2 means transpose
|
||||
! even for complex data.
|
||||
!
|
||||
if (psb_z_is_complex_) then
|
||||
info = amg_zumf_solve(2,n_row,ww,x,n_row,sv%numeric)
|
||||
else
|
||||
info = amg_zumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
end if
|
||||
case('C')
|
||||
info = amg_zumf_solve(1,n_row,ww,x,n_row,sv%numeric)
|
||||
case default
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
|
||||
goto 9999
|
||||
end select
|
||||
|
||||
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
|
||||
|
||||
|
||||
if (info /= psb_success_) then
|
||||
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
|
||||
goto 9999
|
||||
endif
|
||||
|
||||
if (n_col > size(work)) then
|
||||
deallocate(ww)
|
||||
endif
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_z_umf_solver_apply
|
||||
|
||||
subroutine amg_z_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
|
||||
& trans,work,wv,info,init,initu)
|
||||
use psb_base_mod
|
||||
use amg_z_umf_solver, amg_protect_name => amg_z_umf_solver_apply_vect
|
||||
implicit none
|
||||
type(psb_desc_type), intent(in) :: desc_data
|
||||
class(amg_z_umf_solver_type), intent(inout) :: sv
|
||||
type(psb_z_vect_type),intent(inout) :: x
|
||||
type(psb_z_vect_type),intent(inout) :: y
|
||||
complex(psb_dpk_),intent(in) :: alpha,beta
|
||||
character(len=1),intent(in) :: trans
|
||||
complex(psb_dpk_),target, intent(inout) :: work(:)
|
||||
type(psb_z_vect_type),intent(inout) :: wv(:)
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
character, intent(in), optional :: init
|
||||
type(psb_z_vect_type),intent(inout), optional :: initu
|
||||
|
||||
integer :: err_act
|
||||
character(len=20) :: name='z_umf_solver_apply_vect'
|
||||
|
||||
call psb_erractionsave(err_act)
|
||||
|
||||
info = psb_success_
|
||||
!
|
||||
! For non-iterative solvers, init and initu are ignored.
|
||||
!
|
||||
|
||||
call x%v%sync()
|
||||
call y%v%sync()
|
||||
call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info)
|
||||
call y%v%set_host()
|
||||
if (info /= 0) goto 9999
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
|
||||
end subroutine amg_z_umf_solver_apply_vect
|
||||
|
||||
subroutine amg_z_umf_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
|
||||
use psb_base_mod
|
||||
use amg_z_umf_solver, amg_protect_name => amg_z_umf_solver_bld
|
||||
|
||||
Implicit None
|
||||
|
||||
! Arguments
|
||||
type(psb_zspmat_type), intent(in), target :: a
|
||||
Type(psb_desc_type), Intent(inout) :: desc_a
|
||||
class(amg_z_umf_solver_type), intent(inout) :: sv
|
||||
integer(psb_ipk_), intent(out) :: info
|
||||
type(psb_zspmat_type), intent(in), target, optional :: b
|
||||
class(psb_z_base_sparse_mat), intent(in), optional :: amold
|
||||
class(psb_z_base_vect_type), intent(in), optional :: vmold
|
||||
class(psb_i_base_vect_type), intent(in), optional :: imold
|
||||
! Local variables
|
||||
type(psb_zspmat_type) :: atmp
|
||||
type(psb_z_csc_sparse_mat) :: acsc
|
||||
integer :: n_row,n_col, nrow_a, nztota
|
||||
type(psb_ctxt_type) :: ctxt
|
||||
integer :: np,me,i, err_act, debug_unit, debug_level
|
||||
character(len=20) :: name='z_umf_solver_bld', ch_err
|
||||
|
||||
info=psb_success_
|
||||
call psb_erractionsave(err_act)
|
||||
debug_unit = psb_get_debug_unit()
|
||||
debug_level = psb_get_debug_level()
|
||||
ctxt = desc_a%get_context()
|
||||
call psb_info(ctxt, me, np)
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' start'
|
||||
|
||||
|
||||
n_row = desc_a%get_local_rows()
|
||||
n_col = desc_a%get_local_cols()
|
||||
|
||||
call a%cscnv(atmp,info,type='coo')
|
||||
call psb_rwextd(n_row,atmp,info,b=b)
|
||||
call atmp%cscnv(info,type='csc',dupl=psb_dupl_add_)
|
||||
call atmp%mv_to(acsc)
|
||||
nrow_a = acsc%get_nrows()
|
||||
nztota = acsc%get_nzeros()
|
||||
! Fix the entres to call C-base UMFPACK.
|
||||
acsc%ia(:) = acsc%ia(:) - 1
|
||||
acsc%icp(:) = acsc%icp(:) - 1
|
||||
info = amg_zumf_fact(nrow_a,nztota,acsc%val,&
|
||||
& acsc%ia,acsc%icp,sv%symbolic,sv%numeric,&
|
||||
& sv%symbsize,sv%numsize)
|
||||
|
||||
if (info /= psb_success_) then
|
||||
info=psb_err_from_subroutine_
|
||||
ch_err='amg_zumf_fact'
|
||||
call psb_errpush(info,name,a_err=ch_err)
|
||||
goto 9999
|
||||
end if
|
||||
|
||||
call acsc%free()
|
||||
call atmp%free()
|
||||
|
||||
if (debug_level >= psb_debug_outer_) &
|
||||
& write(debug_unit,*) me,' ',trim(name),' end'
|
||||
|
||||
call psb_erractionrestore(err_act)
|
||||
return
|
||||
|
||||
9999 call psb_error_handler(err_act)
|
||||
return
|
||||
end subroutine amg_z_umf_solver_bld
|
||||
Reference in New Issue
Block a user