fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-17 03:12:15 +00:00
Code Issues Packages Projects Releases Wiki Activity

psblas3:

Browse Source 
Makefile
 psb_d_bjacprec.f03
 psb_d_diagprec.f03
 psb_d_nullprec.f03
 psb_dprc_aply.f90
 psb_dprecbld.f90
 psb_dprecinit.f90
 psb_dprecset.f90
 psb_prec_type.f03

New preconditioner structure
This commit is contained in:
Salvatore Filippone
2009-10-20 14:25:43 +00:00
parent 89bec2256e
commit 390cbda059
9 changed files with 1811 additions and 222 deletions
Download Patch File Download Diff File Expand all files Collapse all files
prec/Makefile
+4 -5
View File
@@ -2,12 +2,10 @@ include ../Make.inc
LIBDIR=../lib
HERE=.
MODOBJS= psb_prec_type.o psb_prec_mod.o
F90OBJS= psb_dbjac_bld.o psb_dilu_fct.o\
MODOBJS= psb_prec_type.o psb_prec_mod.o psb_d_diagprec.o psb_d_nullprec.o \
psb_d_bjacprec.o
F90OBJS= psb_dilu_fct.o\
psb_dprecbld.o psb_dprecset.o psb_dprecinit.o \
psb_ddiagsc_bld.o \
psb_dprc_aply.o \
psb_dgprec_aply.o psb_dbjac_aply.o\
psb_sbjac_bld.o psb_silu_fct.o\
psb_sprecbld.o psb_sprecset.o psb_sprecinit.o \
psb_sdiagsc_bld.o \
@@ -39,6 +37,7 @@ $(OBJS): $(LIBDIR)/psb_base_mod$(.mod)
$(F90OBJS): $(MODOBJS)
psb_prec_mod.o: psb_prec_type.o
psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o: psb_prec_type.o psb_prec_mod.o
veryclean: clean
/bin/rm -f $(LIBNAME)
prec/psb_d_bjacprec.f03
+560
View File
@@ -0,0 +1,560 @@
module psb_d_bjacprec
use psb_prec_type
type, extends(psb_d_base_prec_type) :: psb_d_bjac_prec_type
integer, allocatable :: iprcparm(:)
type(psb_d_sparse_mat), allocatable :: av(:)
real(psb_dpk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => d_bjac_apply
procedure, pass(prec) :: precbld => d_bjac_precbld
procedure, pass(prec) :: precinit => d_bjac_precinit
procedure, pass(prec) :: d_base_precseti => d_bjac_precseti
procedure, pass(prec) :: d_base_precsetr => d_bjac_precsetr
procedure, pass(prec) :: d_base_precsetc => d_bjac_precsetc
procedure, pass(prec) :: precfree => d_bjac_precfree
procedure, pass(prec) :: precdescr => d_bjac_precdescr
procedure, pass(prec) :: sizeof => d_bjac_sizeof
end type psb_d_bjac_prec_type
character(len=15), parameter, private :: &
& fact_names(0:2)=(/'None ','ILU(n) ',&
& 'ILU(eps) '/)
contains
subroutine d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_d_bjac_prec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: alpha,beta
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
! Local variables
integer :: n_row,n_col
real(psb_dpk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name='d_bjac_prec_apply'
character(len=20) :: ch_err
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row = psb_cd_get_local_rows(desc_data)
n_col = psb_cd_get_local_cols(desc_data)
if (size(x) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/))
goto 9999
end if
if (size(y) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < n_row) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_apply
subroutine d_bjac_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_d_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_null_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_realloc(psb_ifpsz,prec%iprcparm,info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
end if
prec%iprcparm(:) = 0
prec%iprcparm(psb_p_type_) = psb_bjac_
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
prec%iprcparm(psb_ilu_fill_in_) = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precinit
subroutine d_bjac_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
use psb_prec_mod
Implicit None
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_d_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_d_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name='d_bjac_precbld'
character(len=20) :: ch_err
if(psb_get_errstatus() /= 0) return
info = 0
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(prec%av)) then
if (size(prec%av) < psb_bp_ilu_avsz) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
endif
end if
if (.not.allocated(prec%av)) then
allocate(prec%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = nrow_a
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(prec%d)) then
if (size(prec%d) < n_row) then
deallocate(prec%d)
endif
endif
if (.not.allocated(prec%d)) then
allocate(prec%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
t3 = psb_wtime()
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,prec%d,info)
if(info==0) then
call prec%av(psb_l_pr_)%mv_from(lf)
call prec%av(psb_u_pr_)%mv_from(uf)
call prec%av(psb_l_pr_)%set_asb()
call prec%av(psb_u_pr_)%set_asb()
call prec%av(psb_l_pr_)%trim()
call prec%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precbld
subroutine d_bjac_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_bjac_precset'
call psb_erractionsave(err_act)
info = 0
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
select case(what)
case (psb_f_type_)
if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.(prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precseti
subroutine d_bjac_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precsetr
subroutine d_bjac_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precsetc
subroutine d_bjac_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_d_bjac_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, i
character(len=20) :: name='d_bjac_precfree'
call psb_erractionsave(err_act)
info = 0
if (allocated(prec%av)) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
end if
if (allocated(prec%d)) then
deallocate(prec%d,stat=info)
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precfree
subroutine d_bjac_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_d_bjac_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='d_bjac_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
write(iout_,*) 'Block Jacobi with: ',&
& fact_names(prec%iprcparm(psb_f_type_))
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_bjac_precdescr
function d_bjac_sizeof(prec) result(val)
use psb_base_mod
class(psb_d_bjac_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
if (allocated(prec%d)) then
val = val + psb_sizeof_dp * size(prec%d)
endif
if (allocated(prec%av)) then
val = val + psb_sizeof(prec%av(psb_l_pr_))
val = val + psb_sizeof(prec%av(psb_u_pr_))
endif
return
end function d_bjac_sizeof
end module psb_d_bjacprec
prec/psb_d_diagprec.f03
+352
View File
@@ -0,0 +1,352 @@
module psb_d_diagprec
use psb_prec_type
type, extends(psb_d_base_prec_type) :: psb_d_diag_prec_type
real(psb_dpk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => d_diag_apply
procedure, pass(prec) :: precbld => d_diag_precbld
procedure, pass(prec) :: precinit => d_diag_precinit
procedure, pass(prec) :: d_base_precseti => d_diag_precseti
procedure, pass(prec) :: d_base_precsetr => d_diag_precsetr
procedure, pass(prec) :: d_base_precsetc => d_diag_precsetc
procedure, pass(prec) :: precfree => d_diag_precfree
procedure, pass(prec) :: precdescr => d_diag_precdescr
procedure, pass(prec) :: sizeof => d_diag_sizeof
end type psb_d_diag_prec_type
contains
subroutine d_diag_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_d_diag_prec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(in) :: alpha, beta
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='d_diag_prec_apply'
real(psb_dpk_), pointer :: ww(:)
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the DIAG preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < nrow) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(work) >= size(x)) then
ww => work
else
allocate(ww(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4025,name,i_err=(/size(x),0,0,0,0/),a_err='real(psb_dpk_)')
goto 9999
end if
end if
ww(1:nrow) = x(1:nrow)*prec%d(1:nrow)
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (size(work) < size(x)) then
deallocate(ww,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Deallocate')
goto 9999
end if
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_apply
subroutine d_diag_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_d_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_diag_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precinit
subroutine d_diag_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_d_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow,i
character(len=20) :: name='d_diag_precbld'
call psb_erractionsave(err_act)
info = 0
nrow = psb_cd_get_local_cols(desc_a)
if (allocated(prec%d)) then
if (size(prec%d) < nrow) then
deallocate(prec%d,stat=info)
end if
end if
if ((info == 0).and.(.not.allocated(prec%d))) then
allocate(prec%d(nrow), stat=info)
end if
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call a%get_diag(prec%d,info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name, a_err='get_diag')
goto 9999
end if
do i=1,nrow
if (prec%d(i) == dzero) then
prec%d(i) = done
else
prec%d(i) = done/prec%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precbld
subroutine d_diag_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precseti
subroutine d_diag_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precsetr
subroutine d_diag_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precsetc
subroutine d_diag_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_d_diag_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precfree
subroutine d_diag_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_d_diag_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='d_diag_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'Diagonal scaling'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_diag_precdescr
function d_diag_sizeof(prec) result(val)
use psb_base_mod
class(psb_d_diag_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
val = val + psb_sizeof_dp * size(prec%d)
return
end function d_diag_sizeof
end module psb_d_diagprec
prec/psb_d_nullprec.f03
+293
View File
@@ -0,0 +1,293 @@
module psb_d_nullprec
use psb_prec_type
type, extends(psb_d_base_prec_type) :: psb_d_null_prec_type
contains
procedure, pass(prec) :: apply => d_null_apply
procedure, pass(prec) :: precbld => d_null_precbld
procedure, pass(prec) :: precinit => d_null_precinit
procedure, pass(prec) :: d_base_precseti => d_null_precseti
procedure, pass(prec) :: d_base_precsetr => d_null_precsetr
procedure, pass(prec) :: d_base_precsetc => d_null_precsetc
procedure, pass(prec) :: precfree => d_null_precfree
procedure, pass(prec) :: precdescr => d_null_precdescr
procedure, pass(prec) :: sizeof => d_null_sizeof
end type psb_d_null_prec_type
contains
subroutine d_null_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_d_null_prec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(in) :: alpha, beta
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='d_null_prec_apply'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
if (info /= 0 ) then
info = 4010
call psb_errpush(infoi,name,a_err="psb_geaxpby")
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_apply
subroutine d_null_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_d_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_null_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precinit
subroutine d_null_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_d_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow
character(len=20) :: name='d_null_precbld'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precbld
subroutine d_null_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precseti
subroutine d_null_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precsetr
subroutine d_null_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precsetc
subroutine d_null_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_d_null_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precfree
subroutine d_null_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_d_null_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='d_null_precset'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'No preconditioning'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_null_precdescr
function d_null_sizeof(prec) result(val)
use psb_base_mod
class(psb_d_null_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
return
end function d_null_sizeof
end module psb_d_nullprec
prec/psb_dprc_aply.f90
+14 -4
View File
@@ -73,8 +73,13 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
end if
call psb_gprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info)
if (.not.allocated(prec%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
!!$ call psb_gprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info)
call prec%dprec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_)
if (present(work)) then
else
deallocate(work_)
@@ -139,7 +144,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
type(psb_desc_type),intent(in) :: desc_data
type(psb_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(inout) :: x(:)
real(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
@@ -161,12 +166,17 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
trans_='N'
end if
if (.not.allocated(prec%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
allocate(ww(size(x)),w1(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
call psb_dprc_aply(prec,x,ww,desc_data,info,trans_,work=w1)
call prec%dprec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999
x(:) = ww(:)
deallocate(ww,W1)
prec/psb_dprecbld.f90
+60 -49
View File
@@ -81,63 +81,74 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
! ALso should define symbolic names for the preconditioners.
!
call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',&
& psb_diag_,is_legal_prec)
call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(psb_p_type_))
case (psb_noprec_)
! Do nothing.
call psb_cdcpy(desc_a,p%desc_data,info)
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
call psb_errpush(info,name,a_err=ch_err)
if (.false.) then
!!$ call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',&
!!$ & psb_diag_,is_legal_prec)
!!$
!!$ call psb_nullify_desc(p%desc_data)
!!$
!!$ select case(p%iprcparm(psb_p_type_))
!!$ case (psb_noprec_)
!!$ ! Do nothing.
!!$ call psb_cdcpy(desc_a,p%desc_data,info)
!!$ if(info /= 0) then
!!$ info=4010
!!$ ch_err='psb_cdcpy'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ end if
!!$
!!$ case (psb_diag_)
!!$
!!$ call psb_diagsc_bld(a,desc_a,p,upd_,info)
!!$ if(info /= 0) then
!!$ info=4010
!!$ ch_err='psb_diagsc_bld'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ end if
!!$
!!$ case (psb_bjac_)
!!$
!!$ call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
!!$ & psb_f_ilu_n_,is_legal_ml_fact)
!!$
!!$ call psb_bjac_bld(a,desc_a,p,upd_,info)
!!$
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='psb_bjac_bld')
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ info=4010
!!$ ch_err='Unknown psb_p_type_'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$
!!$ end select
else
if (.not.allocated(p%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
call p%dprec%precbld(a,desc_a,info,upd)
if (info /= 0) goto 9999
case (psb_diag_)
call psb_diagsc_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
info=4010
ch_err='psb_diagsc_bld'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case (psb_bjac_)
call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
& psb_f_ilu_n_,is_legal_ml_fact)
call psb_bjac_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')
goto 9999
end if
case default
info=4010
ch_err='Unknown psb_p_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end if
return
end subroutine psb_dprecbld
end subroutine psb_dprecbld
prec/psb_dprecinit.f90
+51 -22
View File
@@ -33,6 +33,9 @@ subroutine psb_dprecinit(p,ptype,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dprecinit
use psb_d_nullprec
use psb_d_diagprec
use psb_d_bjacprec
implicit none
type(psb_dprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
@@ -40,31 +43,57 @@ subroutine psb_dprecinit(p,ptype,info)
info = 0
call psb_realloc(psb_ifpsz,p%iprcparm,info)
if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info)
if (info /= 0) return
if (.false.) then
!!$ call psb_realloc(psb_ifpsz,p%iprcparm,info)
!!$ if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info)
!!$ if (info /= 0) return
!!$
!!$ select case(psb_toupper(ptype(1:len_trim(ptype))))
!!$ case ('NONE','NOPREC')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_noprec_
!!$ p%iprcparm(psb_f_type_) = psb_f_none_
!!$
!!$ case ('DIAG')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_diag_
!!$ p%iprcparm(psb_f_type_) = psb_f_none_
!!$
!!$ case ('BJAC')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_bjac_
!!$ p%iprcparm(psb_f_type_) = psb_f_ilu_n_
!!$ p%iprcparm(psb_ilu_fill_in_) = 0
!!$
!!$ case default
!!$ write(0,*) 'Unknown preconditioner type request "',ptype,'"'
!!$ info = 2
!!$
!!$ end select
else
if (allocated(p%dprec) ) then
call p%dprec%precfree(info)
if (info == 0) deallocate(p%dprec,stat=info)
if (info /= 0) return
end if
select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_noprec_
p%iprcparm(psb_f_type_) = psb_f_none_
select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
allocate(psb_d_null_prec_type :: p%dprec, stat=info)
case ('DIAG')
allocate(psb_d_diag_prec_type :: p%dprec, stat=info)
case ('DIAG')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_diag_
p%iprcparm(psb_f_type_) = psb_f_none_
case ('BJAC')
allocate(psb_d_bjac_prec_type :: p%dprec, stat=info)
case ('BJAC')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_bjac_
p%iprcparm(psb_f_type_) = psb_f_ilu_n_
p%iprcparm(psb_ilu_fill_in_) = 0
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2
end select
end select
if (info == 0) call p%dprec%precinit(info)
end if
end subroutine psb_dprecinit
prec/psb_dprecset.f90
+63 -45
View File
@@ -37,50 +37,18 @@ subroutine psb_dprecseti(p,what,val,info)
type(psb_dprec_type), intent(inout) :: p
integer :: what, val
integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0
return
select case(what)
case (psb_f_type_)
if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_f_type_) = val
if (.not.allocated(p%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
case (psb_ilu_fill_in_)
if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return
end subroutine psb_dprecseti
subroutine psb_dprecsetd(p,what,val,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dprecsetd
implicit none
type(psb_dprec_type), intent(inout) :: p
integer :: what
real(psb_dpk_) :: val
integer, intent(out) :: info
!
! This will have to be changed if/when we put together an ILU(eps)
! factorization.
!
select case(what)
call p%dprec%precset(what,val,info)
!!$ select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
@@ -96,11 +64,61 @@ subroutine psb_dprecsetd(p,what,val,info)
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
!!$
!!$ case default
!!$ write(0,*) 'WHAT is invalid, ignoring user specification'
!!$
!!$ end select
return
end select
end subroutine psb_dprecseti
subroutine psb_dprecsetd(p,what,val,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dprecsetd
implicit none
type(psb_dprec_type), intent(inout) :: p
integer :: what
real(psb_dpk_) :: val
integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0
if (.not.allocated(p%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
call p%dprec%precset(what,val,info)
!!$!
!!$! This will have to be changed if/when we put together an ILU(eps)
!!$! factorization.
!!$!
!!$ select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_f_type_) = val
!!$
!!$ case (psb_ilu_fill_in_)
!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
!!$
!!$ case default
!!$ write(0,*) 'WHAT is invalid, ignoring user specification'
!!$
!!$ end select
return
end subroutine psb_dprecsetd
prec/psb_prec_type.f03
+414 -97
View File
@@ -78,14 +78,23 @@ module psb_prec_type
generic, public :: apply => s_apply2v, s_apply1v
end type psb_sprec_type
type psb_d_base_prec_type
contains
procedure, pass(prec) :: apply => d_base_apply
procedure, pass(prec) :: precbld => d_base_precbld
procedure, pass(prec) :: d_base_precseti
procedure, pass(prec) :: d_base_precsetr
procedure, pass(prec) :: d_base_precsetc
procedure, pass(prec) :: sizeof => d_base_sizeof
generic, public :: precset => d_base_precseti, d_base_precsetr, d_base_precsetc
procedure, pass(prec) :: precinit => d_base_precinit
procedure, pass(prec) :: precfree => d_base_precfree
procedure, pass(prec) :: precdescr => d_base_precdescr
end type psb_d_base_prec_type
type psb_dprec_type
type(psb_d_sparse_mat), allocatable :: av(:)
real(psb_dpk_), allocatable :: d(:)
type(psb_desc_type) :: desc_data
integer, allocatable :: iprcparm(:)
real(psb_dpk_), allocatable :: rprcparm(:)
integer, allocatable :: perm(:), invperm(:)
integer :: prec
class(psb_d_base_prec_type), allocatable :: dprec
contains
procedure, pass(prec) :: d_apply2v
procedure, pass(prec) :: d_apply1v
@@ -148,8 +157,10 @@ module psb_prec_type
end interface
interface psb_sizeof
module procedure psb_sprec_sizeof, psb_dprec_sizeof,&
& psb_cprec_sizeof, psb_zprec_sizeof
module procedure psb_sprec_sizeof, &
& psb_dprec_sizeof,&
& psb_cprec_sizeof, &
& psb_zprec_sizeof
end interface
@@ -161,8 +172,8 @@ module psb_prec_type
import psb_sprec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_sprec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
@@ -192,7 +203,7 @@ module psb_prec_type
import psb_dprec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(inout) :: x(:)
real(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine psb_dprc_aply1
@@ -201,8 +212,8 @@ module psb_prec_type
import psb_cprec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
@@ -212,7 +223,7 @@ module psb_prec_type
import psb_cprec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(inout) :: x(:)
complex(psb_spk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine psb_cprc_aply1
@@ -221,8 +232,8 @@ module psb_prec_type
import psb_zprec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
@@ -232,7 +243,7 @@ module psb_prec_type
import psb_zprec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(inout) :: x(:)
complex(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
end subroutine psb_zprc_aply1
@@ -245,9 +256,11 @@ contains
subroutine psb_file_prec_descr(p,iout)
use psb_base_mod
type(psb_dprec_type), intent(in) :: p
integer, intent(in), optional :: iout
integer :: iout_
character(len=20) :: name='prec_descr'
if (present(iout)) then
iout_ = iout
@@ -255,17 +268,12 @@ contains
iout_ = 6
end if
write(iout_,*) 'Preconditioner description'
select case(p%iprcparm(psb_p_type_))
case(psb_noprec_)
write(iout_,*) 'No preconditioning'
case(psb_diag_)
write(iout_,*) 'Diagonal scaling'
case(psb_bjac_)
write(iout_,*) 'Block Jacobi with: ',&
& fact_names(p%iprcparm(psb_f_type_))
end select
if (.not.allocated(p%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
end if
call p%dprec%precdescr(iout)
end subroutine psb_file_prec_descr
subroutine psb_sfile_prec_descr(p,iout)
@@ -496,40 +504,13 @@ contains
me=-1
! Actually we migh just deallocate the top level array, except
! for the inner UMFPACK or SLU stuff
if (allocated(p%d)) then
deallocate(p%d,stat=info)
if (allocated(p%dprec)) then
call p%dprec%precfree(info)
if (info /= 0) goto 9999
deallocate(p%dprec,stat=info)
if (info /= 0) goto 9999
end if
if (allocated(p%av)) then
do i=1,size(p%av)
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
end if
if (allocated(p%desc_data%matrix_data)) &
& call psb_cdfree(p%desc_data,info)
if (allocated(p%rprcparm)) then
deallocate(p%rprcparm,stat=info)
end if
if (allocated(p%perm)) then
deallocate(p%perm,stat=info)
endif
if (allocated(p%invperm)) then
deallocate(p%invperm,stat=info)
endif
if (allocated(p%iprcparm)) then
deallocate(p%iprcparm,stat=info)
end if
call psb_nullify_prec(p)
call psb_erractionrestore(err_act)
return
@@ -546,8 +527,6 @@ contains
subroutine psb_nullify_dprec(p)
type(psb_dprec_type), intent(inout) :: p
!!$ nullify(p%av,p%d,p%iprcparm,p%rprcparm,p%perm,p%invperm,p%mlia,&
!!$ & p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac)
end subroutine psb_nullify_dprec
@@ -695,18 +674,10 @@ contains
integer(psb_long_int_k_) :: val
integer :: i
val = 0
if (allocated(prec%iprcparm)) val = val + psb_sizeof_int * size(prec%iprcparm)
if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm)
if (allocated(prec%d)) val = val + psb_sizeof_dp * size(prec%d)
if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm)
if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm)
val = val + psb_sizeof(prec%desc_data)
if (allocated(prec%av)) then
do i=1,size(prec%av)
val = val + psb_sizeof(prec%av(i))
end do
end if
if (allocated(prec%dprec)) then
val = val + prec%dprec%sizeof()
end if
end function psb_dprec_sizeof
function psb_sprec_sizeof(prec) result(val)
@@ -851,23 +822,56 @@ contains
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act
character(len=20) :: name='d_prec_apply'
character :: trans_
real(psb_dpk_), pointer :: work_(:)
integer :: ictxt,np,me,err_act
character(len=20) :: name
name='psb_d_prc_apply'
info = 0
call psb_erractionsave(err_act)
select type(prec)
type is (psb_dprec_type)
call psb_precaply(prec,x,y,desc_data,info,trans,work)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
if (present(trans)) then
trans_=trans
else
trans_='N'
end if
if (present(work)) then
work_ => work
else
allocate(work_(4*psb_cd_get_local_cols(desc_data)),stat=info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
end if
if (.not.allocated(prec%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
call prec%dprec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_)
if (present(work)) then
else
deallocate(work_,stat=info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name,a_err='DeAllocate')
goto 9999
end if
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
@@ -885,24 +889,51 @@ contains
real(psb_dpk_),intent(inout) :: x(:)
integer, intent(out) :: info
character(len=1), optional :: trans
Integer :: err_act
character(len=20) :: name='d_prec_apply'
character :: trans_
integer :: ictxt,np,me, err_act
real(psb_dpk_), pointer :: WW(:), w1(:)
character(len=20) :: name
name='psb_d_apply1'
info = 0
call psb_erractionsave(err_act)
select type(prec)
type is (psb_dprec_type)
call psb_precaply(prec,x,desc_data,info,trans)
class default
info = 700
call psb_errpush(info,name)
goto 9999
end select
ictxt=psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
if (present(trans)) then
trans_=psb_toupper(trans)
else
trans_='N'
end if
if (.not.allocated(prec%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
allocate(ww(size(x)),w1(size(x)),stat=info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name,a_err='Allocate')
goto 9999
end if
call prec%dprec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999
x(:) = ww(:)
deallocate(ww,W1,stat=info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name,a_err='DeAllocate')
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
@@ -1052,4 +1083,290 @@ contains
end subroutine z_apply1v
subroutine d_base_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_d_base_prec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: alpha, beta
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='d_base_prec_apply'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_apply
subroutine d_base_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_d_base_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_base_precinit'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precinit
subroutine d_base_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_d_base_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow
character(len=20) :: name='d_base_precbld'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precbld
subroutine d_base_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_base_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_base_precseti'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precseti
subroutine d_base_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_base_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_base_precsetr'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precsetr
subroutine d_base_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_d_base_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_base_precsetc'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precsetc
subroutine d_base_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_d_base_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='d_base_precfree'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precfree
subroutine d_base_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_d_base_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='d_base_precdescr'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 700
call psb_errpush(info,name)
goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine d_base_precdescr
function d_base_sizeof(prec) result(val)
use psb_base_mod
class(psb_d_base_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
return
end function d_base_sizeof
end module psb_prec_type