Fixes for sundial interfacing to VECT and LINSOLVE (rgmres)

kinsol-stop
sfilippone 4 days ago
parent a988b33cd5
commit 18a9ebbd8b

@ -1483,10 +1483,10 @@ contains
integer(psb_ipk_) :: i, n
info = 0
if (z%is_dev()) call z%sync()
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call z%div(x%v,y%v,info)
call x%set_host()
end subroutine c_base_div_v2
!
!> Function base_div_v_check
@ -1506,6 +1506,7 @@ contains
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call x%div(x%v,y%v,info,flag)
end subroutine c_base_div_v_check

@ -1490,10 +1490,10 @@ contains
integer(psb_ipk_) :: i, n
info = 0
if (z%is_dev()) call z%sync()
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call z%div(x%v,y%v,info)
call x%set_host()
end subroutine d_base_div_v2
!
!> Function base_div_v_check
@ -1513,6 +1513,7 @@ contains
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call x%div(x%v,y%v,info,flag)
end subroutine d_base_div_v_check

@ -418,8 +418,6 @@ contains
class(psb_d_vect_type), intent(in) :: x
integer(psb_ipk_) :: res
res = 0
write(0,*) allocated(x%v)
if (allocated(x%v)) write(0,*) allocated(x%v%v)
if (allocated(x%v)) res = x%v%get_nrows()
end function d_vect_get_nrows

@ -1490,10 +1490,10 @@ contains
integer(psb_ipk_) :: i, n
info = 0
if (z%is_dev()) call z%sync()
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call z%div(x%v,y%v,info)
call x%set_host()
end subroutine s_base_div_v2
!
!> Function base_div_v_check
@ -1513,6 +1513,7 @@ contains
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call x%div(x%v,y%v,info,flag)
end subroutine s_base_div_v_check

@ -1483,10 +1483,10 @@ contains
integer(psb_ipk_) :: i, n
info = 0
if (z%is_dev()) call z%sync()
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call z%div(x%v,y%v,info)
call x%set_host()
end subroutine z_base_div_v2
!
!> Function base_div_v_check
@ -1506,6 +1506,7 @@ contains
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
call x%div(x%v,y%v,info,flag)
end subroutine z_base_div_v_check

@ -21,32 +21,19 @@ contains
integer(psb_c_ipk_) :: info
res = -1
nullify(xp)
if (c_associated(cdh%item)) then
call c_f_pointer(cdh%item,descp)
else
return
end if
!!$ write(0,*) 'On entry to C_DGEALL xh->item'
call psb_c_print_pointer(xh%item)
if (c_associated(xh%item)) then
write(0,*) 'C associated on c_dgeall'
call psb_c_print_pointer(xh%item)
return
end if
allocate(xp,stat=info)
write(0,*) 'From DGEALL/ALLOCATE:',info
call psb_c_print_pointer(c_loc(xp))
allocate(xp)
call psb_geall(xp,descp,info)
xp%v%v(1) = 1.d0
write(0,*) 'c_dgeall out from geall xp ',xp%get_nrows()
xh%item = c_loc(xp)
res = min(0,info)
!!$ write(0,*) 'Check from C_DGEALL 1:',info
write(0,*) 'On end of C_DGEALL xh->item'
call psb_c_print_pointer(xh%item)
!!$ if (info==0) write(0,*) 'Check from C_DGEALL 2:',xp%get_nrows(),descp%get_local_cols()
return
end function psb_c_dgeall
@ -70,8 +57,6 @@ contains
return
end if
if (c_associated(xh%item)) then
write(0,*) 'C associated on c_dgeall_remote'
call psb_c_print_pointer(xh%item)
return
end if
allocate(xp)

@ -6,6 +6,7 @@ module psb_base_linsolve_cbind_mod
type, bind(c) :: solveroptions
integer(psb_c_ipk_) :: iter, itmax, itrace, irst, istop
real(c_double) :: eps, err
type(psb_c_object_type) :: s1, s2
end type solveroptions
contains
@ -20,7 +21,8 @@ contains
options%istop = 2
options%irst = 10
options%eps = 1.d-6
options%s1 = psb_c_get_new_object()
options%s2 = psb_c_get_new_object()
res = 0
end function psb_c_DefaultSolverOptions
@ -30,12 +32,15 @@ contains
type(solveroptions), value :: options
integer(psb_c_ipk_) :: res
write(*,*) 'PSBLAS C Interface Solver Options '
write(*,*) ' Maximum number of iterations :', options%itmax
write(*,*) ' Tracing :', options%itrace
write(*,*) ' Stopping Criterion :', options%istop
write(*,*) ' Restart :', options%irst
write(*,*) ' EPS (tolerance) :', options%eps
write(0,*) 'PSBLAS C Interface Solver Options '
write(0,*) ' Maximum number of iterations :', options%itmax
write(0,*) ' Tracing :', options%itrace
write(0,*) ' Stopping Criterion :', options%istop
write(0,*) ' Restart :', options%irst
write(0,*) ' EPS (tolerance) :', options%eps
write(0,*) ' S1 scaling :', c_associated(options%s1%item)
write(0,*) ' S2 scaling :', c_associated(options%s2%item)
res = 0
end function psb_c_PrintSolverOptions

@ -24,13 +24,14 @@ contains
res= psb_c_ckrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_ckrylov
function psb_c_ckrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istops1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
@ -49,10 +50,12 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_cspmat_type), pointer :: ap
type(psb_cprec_type), pointer :: precp
type(psb_c_vect_type), pointer :: xp, bp
type(psb_c_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
@ -84,6 +87,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call stringc2f(methd,fmethd)
@ -94,10 +107,27 @@ contains
fistop = istop
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info

@ -24,13 +24,14 @@ contains
res= psb_c_dkrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_dkrylov
function psb_c_dkrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istops1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
@ -49,10 +50,12 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_dspmat_type), pointer :: ap
type(psb_dprec_type), pointer :: precp
type(psb_d_vect_type), pointer :: xp, bp
type(psb_d_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
@ -84,6 +87,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call stringc2f(methd,fmethd)
@ -94,10 +107,27 @@ contains
fistop = istop
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info

@ -22,6 +22,8 @@ typedef struct psb_c_solveroptions {
int istop; /* Stopping criterion: 1:backward error 2: ||r||_2/||b||_2 */
double eps; /* Stopping tolerance */
double err; /* Convergence indicator on exit */
void *s1;
void *s2;
} psb_c_SolverOptions;
int psb_c_DefaultSolverOptions(psb_c_SolverOptions *opt);

@ -24,13 +24,14 @@ contains
res= psb_c_skrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_skrylov
function psb_c_skrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istops1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
@ -49,10 +50,12 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_sspmat_type), pointer :: ap
type(psb_sprec_type), pointer :: precp
type(psb_s_vect_type), pointer :: xp, bp
type(psb_s_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
@ -84,6 +87,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call stringc2f(methd,fmethd)
@ -94,10 +107,27 @@ contains
fistop = istop
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info

@ -24,13 +24,14 @@ contains
res= psb_c_zkrylov_opt(methd, ah, ph, bh, xh, options%eps,cdh, &
& itmax=options%itmax, iter=options%iter,&
& itrace=options%itrace, istop=options%istop,&
& irst=options%irst, err=options%err)
& irst=options%irst, err=options%err, s1=options%s1,s2=options%s2)
end function psb_c_zkrylov
function psb_c_zkrylov_opt(methd,&
& ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
& ah,ph,bh,xh,eps,cdh,itmax,iter,&
& err,itrace,irst,istops1,s2) bind(c) result(res)
use psb_base_mod
use psb_error_mod
use psb_prec_mod
@ -49,10 +50,12 @@ contains
integer(psb_c_ipk_) :: iter
real(c_double) :: err
character(c_char) :: methd(*)
type(psb_c_object_type) :: s1,s2
type(psb_desc_type), pointer :: descp
type(psb_zspmat_type), pointer :: ap
type(psb_zprec_type), pointer :: precp
type(psb_z_vect_type), pointer :: xp, bp
type(psb_z_vect_type), pointer :: xp, bp, s1p, s2p
integer(psb_c_ipk_) :: info,fitmax,fitrace,first,fistop,fiter,err_act
character(len=20) :: fmethd
@ -84,6 +87,16 @@ contains
else
return
end if
if (c_associated(s1%item)) then
call c_f_pointer(s1%item,s1p)
else
nullify(s1p)
end if
if (c_associated(s2%item)) then
call c_f_pointer(s2%item,s2p)
else
nullify(s2p)
end if
call stringc2f(methd,fmethd)
@ -94,10 +107,27 @@ contains
fistop = istop
err_act = psb_act_abort_
if (psb_errstatus_fatal()) call psb_error_handler(err_act)
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
if (associated(s1p).and.associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p,s2=s2p)
else if (associated(s1p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s1=s1p)
else if (associated(s2p)) then
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr,s2=s2p)
else
call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
& descp, info,&
& itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
& irst=first, err=ferr)
end if
iter = fiter
err = ferr
res = info

@ -382,8 +382,14 @@ subroutine psb_crgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(cone,a,w1,czero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
do k = 1, i

@ -382,8 +382,14 @@ subroutine psb_drgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(done,a,w1,dzero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
do k = 1, i

@ -382,8 +382,14 @@ subroutine psb_srgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(sone,a,w1,szero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
do k = 1, i

@ -382,8 +382,14 @@ subroutine psb_zrgmres_vect(a,prec,b,x,eps,desc_a,info,&
inner: Do i=1,nl
itx = itx + 1
call prec%apply(v(i),w1,desc_a,info)
if (present(s2)) then
call psb_gediv(v(i),s2,w,desc_a,info)
call prec%apply(w,w1,desc_a,info)
else
call prec%apply(v(i),w1,desc_a,info)
end if
call psb_spmm(zone,a,w1,zzero,w,desc_a,info,work=aux)
if (present(s1)) call psb_gemlt(s1,w,desc_a,info)
!
do k = 1, i

Loading…
Cancel
Save