Merge branch 'maint-1.2.0' of github.com:sfilippone/amg4psblas into maint-1.2.0

amgprec/amg_c_jac_smoother.f90

@@ -345,7 +345,6 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
-      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_c_prec_type.f90

@@ -665,6 +665,10 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
+    !
+    ! In the internals, do FREE on components,
+    ! but do not deallocate them
+    !
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free_smoothers(info)

amgprec/amg_d_jac_smoother.f90

@@ -345,7 +345,6 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
-      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_d_poly_smoother.f90

@@ -279,7 +279,6 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
-      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_d_prec_type.f90

@@ -665,6 +665,10 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
+    !
+    ! In the internals, do FREE on components,
+    ! but do not deallocate them
+    !
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free_smoothers(info)

amgprec/amg_s_jac_smoother.f90

@@ -345,7 +345,6 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
-      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_s_poly_smoother.f90

@@ -279,7 +279,6 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
-      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_s_prec_type.f90

@@ -665,6 +665,10 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
+    !
+    ! In the internals, do FREE on components,
+    ! but do not deallocate them
+    !
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free_smoothers(info)

amgprec/amg_z_jac_smoother.f90

@@ -345,7 +345,6 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
-      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_z_prec_type.f90

@@ -665,6 +665,10 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
+    !
+    ! In the internals, do FREE on components,
+    ! but do not deallocate them
+    !
     if (allocated(prec%precv)) then
       do i=1,size(prec%precv)
         call prec%precv(i)%free_smoothers(info)

amgprec/impl/smoother/amg_c_as_smoother_free.f90

@@ -53,7 +53,6 @@ subroutine amg_c_as_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info)
     if (info /= psb_success_) then 
       info = psb_err_alloc_dealloc_
       call psb_errpush(info,name)

amgprec/impl/smoother/amg_c_base_smoother_free.f90

@@ -52,7 +52,6 @@ subroutine amg_c_base_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info) 
   end if
   if (info /= psb_success_) then 
     info = psb_err_alloc_dealloc_

amgprec/impl/smoother/amg_d_as_smoother_free.f90

@@ -53,7 +53,6 @@ subroutine amg_d_as_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info)
     if (info /= psb_success_) then 
       info = psb_err_alloc_dealloc_
       call psb_errpush(info,name)

amgprec/impl/smoother/amg_d_base_smoother_free.f90

@@ -52,7 +52,6 @@ subroutine amg_d_base_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info) 
   end if
   if (info /= psb_success_) then 
     info = psb_err_alloc_dealloc_

amgprec/impl/smoother/amg_s_as_smoother_free.f90

@@ -53,7 +53,6 @@ subroutine amg_s_as_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info)
     if (info /= psb_success_) then 
       info = psb_err_alloc_dealloc_
       call psb_errpush(info,name)

amgprec/impl/smoother/amg_s_base_smoother_free.f90

@@ -52,7 +52,6 @@ subroutine amg_s_base_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info) 
   end if
   if (info /= psb_success_) then 
     info = psb_err_alloc_dealloc_

amgprec/impl/smoother/amg_z_as_smoother_free.f90

@@ -53,7 +53,6 @@ subroutine amg_z_as_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info)
     if (info /= psb_success_) then 
       info = psb_err_alloc_dealloc_
       call psb_errpush(info,name)

amgprec/impl/smoother/amg_z_base_smoother_free.f90

@@ -52,7 +52,6 @@ subroutine amg_z_base_smoother_free(sm,info)
 
   if (allocated(sm%sv)) then 
     call sm%sv%free(info)
-    if (info == psb_success_) deallocate(sm%sv,stat=info) 
   end if
   if (info /= psb_success_) then 
     info = psb_err_alloc_dealloc_

amgprec/impl/solver/amg_c_diag_solver_bld.f90

@@ -91,7 +91,7 @@ subroutine amg_c_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)
@@ -172,7 +172,7 @@ subroutine amg_c_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_c_jac_solver_bld.f90

@@ -100,7 +100,7 @@ subroutine amg_c_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_c_l1_jac_solver_bld.f90

@@ -103,7 +103,7 @@ subroutine amg_c_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = cone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_d_diag_solver_bld.f90

@@ -91,7 +91,7 @@ subroutine amg_d_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)
@@ -172,7 +172,7 @@ subroutine amg_d_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_d_jac_solver_bld.f90

@@ -100,7 +100,7 @@ subroutine amg_d_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_d_l1_jac_solver_bld.f90

@@ -103,7 +103,7 @@ subroutine amg_d_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = done/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_s_diag_solver_bld.f90

@@ -91,7 +91,7 @@ subroutine amg_s_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)
@@ -172,7 +172,7 @@ subroutine amg_s_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_s_jac_solver_bld.f90

@@ -100,7 +100,7 @@ subroutine amg_s_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_s_l1_jac_solver_bld.f90

@@ -103,7 +103,7 @@ subroutine amg_s_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = sone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_z_diag_solver_bld.f90

@@ -91,7 +91,7 @@ subroutine amg_z_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)
@@ -172,7 +172,7 @@ subroutine amg_z_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_z_jac_solver_bld.f90

@@ -100,7 +100,7 @@ subroutine amg_z_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

amgprec/impl/solver/amg_z_l1_jac_solver_bld.f90

@@ -103,7 +103,7 @@ subroutine amg_z_l1_jac_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
       sv%d(i) = zone/sv%d(i)
     end if
   end do
-  allocate(sv%dv,stat=info) 
+  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
   if (info == psb_success_) then 
     call sv%dv%bld(sv%d)
     if (present(vmold)) call sv%dv%cnv(vmold)

cbind/amgprec/amg_c_dprec.h

@@ -34,8 +34,9 @@ extern "C" {
   psb_i_t amg_c_ddescr(amg_c_dprec *ph);
 
   psb_i_t amg_c_dkrylov(const char *method, psb_c_dspmat *ah, amg_c_dprec *ph, 
-		  psb_c_dvector *bh, psb_c_dvector *xh,
-		  psb_c_descriptor *cdh, psb_c_SolverOptions *opt);
+			psb_c_dvector *bh, psb_c_dvector *xh,
+			psb_c_descriptor *cdh, psb_c_dvector *s1,
+			psb_c_dvector *s2, psb_c_SolverOptions *opt);
 
 
 #ifdef __cplusplus

cbind/amgprec/amg_c_zprec.h

@@ -35,8 +35,9 @@ extern "C" {
   psb_i_t amg_c_zdescr(amg_c_zprec *ph);
 
   psb_i_t amg_c_zkrylov(const char *method, psb_c_zspmat *ah, amg_c_zprec *ph, 
-		  psb_c_zvector *bh, psb_c_zvector *xh,
-		  psb_c_descriptor *cdh, psb_c_SolverOptions *opt);
+			psb_c_zvector *bh, psb_c_zvector *xh,
+			psb_c_descriptor *cdh, psb_c_zvector *s1,
+			psb_c_zvector *s2,  psb_c_SolverOptions *opt);
 
 
 #ifdef __cplusplus

cbind/amgprec/amg_dprec_cbind_mod.F90

@@ -3,6 +3,7 @@ module amg_dprec_cbind_mod
   use iso_c_binding
   use amg_prec_mod
   use psb_base_cbind_mod
+  use psb_dlinsolve_cbind_mod
 
   type, bind(c) :: amg_c_dprec
     type(c_ptr) :: item = c_null_ptr
@@ -18,7 +19,7 @@ contains
 #define AMGC_ERROR(MSG)
 #endif
 #define amg_success_ 0
-!#define AMGC_ERR_FILTER(INFO) min(0,INFO)
+  !#define AMGC_ERR_FILTER(INFO) min(0,INFO)
 #define AMGC_ERR_FILTER(INFO) (INFO)
 #define AMGC_ERR_HANDLE(INFO) if(INFO/=amg_success_)AMGC_ERROR("ERROR!")
   function  amg_c_dprecinit(cctxt,ph,ptype) bind(c) result(res)
@@ -172,15 +173,17 @@ contains
   end function amg_c_dprecbld
 
   function  amg_c_dhierarchy_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_dprec_type), pointer  :: precp
     type(psb_dspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret, act
+
 
     res = -1
 
@@ -204,11 +207,15 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
-
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
     return
   end function amg_c_dhierarchy_build
 
   function  amg_c_dsmoothers_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
@@ -217,7 +224,7 @@ contains
     type(psb_dspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
-    integer(psb_ipk_)     :: iret
+    integer(psb_ipk_)     :: iret, act
 
     res = -1
 
@@ -241,12 +248,16 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
 
     return
   end function amg_c_dsmoothers_build
 
   function  amg_c_dkrylov(methd,&
-       & ah,ph,bh,xh,cdh,options) bind(c) result(res)
+       & ah,ph,bh,xh,cdh,s1,s2,options) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
     use psb_linsolve_mod
@@ -254,20 +265,20 @@ contains
     use psb_dlinsolve_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1,s2
     character(c_char)       :: methd(*)
     type(solveroptions)     :: options
-
     res= amg_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=s1,s2=s2)
 
   end function amg_c_dkrylov
 
 
   function  amg_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,istop,s1,s2) bind(c) result(res)
     use psb_base_mod
     use psb_prec_mod
     use psb_linsolve_mod
@@ -276,16 +287,17 @@ contains
     use psb_base_string_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1, s2
     integer(psb_c_ipk_), value :: itmax,itrace,irst,istop
     real(c_double), value :: eps
     integer(psb_c_ipk_)        :: iter
     real(c_double)        :: err
     character(c_char)       :: methd(*)
+
     type(psb_desc_type), pointer   :: descp
     type(psb_dspmat_type), pointer :: ap
     type(amg_dprec_type), pointer  :: precp
-    type(psb_d_vect_type), pointer :: xp, bp
+    type(psb_d_vect_type), pointer :: xp, bp, s1p, s2p
 
     integer(psb_ipk_)     :: iret,fitmax,fitrace,first,fistop,fiter
     character(len=20)     :: fmethd
@@ -317,8 +329,17 @@ 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
+    !write(0,*) 'krylov: s1 s2 ',associated(s1p),associated(s2p)
     call stringc2f(methd,fmethd)
     feps    = eps
     fitmax  = itmax
@@ -326,10 +347,28 @@ contains
     first   = irst
     fistop  = istop
 
-    call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
-         & descp, iret,&
-         & 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, iret,&
+           & 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, iret,&
+           & 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, iret,&
+           & 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, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr)
+    end if
+    
     iter = fiter
     err  = ferr
     res = min(iret,0)
@@ -384,7 +423,7 @@ contains
     res = AMGC_ERR_FILTER(info)
     AMGC_ERR_HANDLE(res)
     return
-end function amg_c_dprecapply
+  end function amg_c_dprecapply
 
   function amg_c_dprecapply_opt(ph,bc,xc,cdh,ctrans) bind(c,name="amg_c_dprecapply_opt") result(res)
     use psb_base_mod
@@ -432,7 +471,7 @@ end function amg_c_dprecapply
 
     ! Convert transpose flag
     call stringc2f(ctrans,ftrans)
-    
+
     ! Apply preconditioner
     call precp%apply(bp,xp,descp,info,trans=ftrans)
 
@@ -440,7 +479,7 @@ end function amg_c_dprecapply
     res = AMGC_ERR_FILTER(info)
     AMGC_ERR_HANDLE(res)
     return
-end function amg_c_dprecapply_opt
+  end function amg_c_dprecapply_opt
 
   function  amg_c_dprecfree(ph) bind(c) result(res)
     implicit none

cbind/amgprec/amg_zprec_cbind_mod.F90

@@ -3,6 +3,7 @@ module amg_zprec_cbind_mod
   use iso_c_binding
   use amg_prec_mod
   use psb_base_cbind_mod
+  use psb_zlinsolve_cbind_mod
 
   type, bind(c) :: amg_c_zprec
     type(c_ptr) :: item = c_null_ptr
@@ -140,11 +141,11 @@ contains
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_zprec_type), pointer  :: precp
     type(psb_zspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret
 
     res = -1
 
@@ -173,15 +174,17 @@ contains
   end function amg_c_zprecbld
 
   function  amg_c_zhierarchy_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_zprec_type), pointer  :: precp
     type(psb_zspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret, act
+
 
     res = -1
 
@@ -205,20 +208,25 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
-
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
+      
     return
   end function amg_c_zhierarchy_build
 
   function  amg_c_zsmoothers_build(ah,cdh,ph) bind(c) result(res)
+    use psb_base_mod
     implicit none
 
     integer(psb_c_ipk_) :: res
     type(psb_c_object_type)  :: ph,ah,cdh
-    integer(psb_ipk_)     :: iret
     type(amg_zprec_type), pointer  :: precp
     type(psb_zspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
+    integer(psb_ipk_)     :: iret, act
 
     res = -1
 
@@ -242,37 +250,47 @@ contains
 
     res = AMGC_ERR_FILTER(iret)
     AMGC_ERR_HANDLE(res)
+    if (res /=0) then
+      act = psb_act_abort_
+      call psb_error_handler(act)
+    end if
 
     return
   end function amg_c_zsmoothers_build
 
   function  amg_c_zkrylov(methd,&
-       & ah,ph,bh,xh,cdh,options) bind(c) result(res)
+       & ah,ph,bh,xh,cdh,s1,s2,options) bind(c) result(res)
+    use psb_base_mod
+    use psb_prec_mod
     use psb_linsolve_mod
     use psb_prec_cbind_mod
     use psb_zlinsolve_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1,s2
     character(c_char)       :: methd(*)
     type(solveroptions)     :: options
 
     res= amg_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=s1,s2=s2)
 
   end function amg_c_zkrylov
 
 
   function  amg_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,istop,s1,s2) bind(c) result(res)
+    use psb_base_mod
+    use psb_prec_mod
     use psb_linsolve_mod
     use psb_objhandle_mod
+    use psb_prec_cbind_mod
     use psb_base_string_cbind_mod
     implicit none
     integer(psb_c_ipk_)          :: res
-    type(psb_c_object_type) :: ah,cdh,ph,bh,xh
+    type(psb_c_object_type) :: ah,cdh,ph,bh,xh,s1, s2
     integer(psb_c_ipk_), value :: itmax,itrace,irst,istop
     real(c_double), value :: eps
     integer(psb_c_ipk_)        :: iter
@@ -281,7 +299,7 @@ contains
     type(psb_desc_type), pointer   :: descp
     type(psb_zspmat_type), pointer :: ap
     type(amg_zprec_type), pointer  :: precp
-    type(psb_z_vect_type), pointer :: xp, bp
+    type(psb_z_vect_type), pointer :: xp, bp, s1p, s2p
 
     integer(psb_ipk_)     :: iret,fitmax,fitrace,first,fistop,fiter
     character(len=20)     :: fmethd
@@ -313,7 +331,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)
     feps    = eps
@@ -322,10 +349,28 @@ contains
     first   = irst
     fistop  = istop
 
-    call psb_krylov(fmethd, ap, precp, bp, xp, feps, &
-         & descp, iret,&
-         & 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, iret,&
+           & 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, iret,&
+           & 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, iret,&
+           & 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, iret,&
+           & itmax=fitmax,iter=fiter,itrace=fitrace,istop=fistop,&
+           & irst=first, err=ferr)
+    end if
+
     iter = fiter
     err  = ferr
     res = min(iret,0)

samples/advanced/pdegen/amg_d_genpde_mod.F90

@@ -176,7 +176,7 @@ contains
     else
       partition_ = 3
     end if
-    deltah   = done/(idim+2)
+    deltah   = done/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_dpk_* deltah
 
@@ -412,9 +412,9 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
           ! x, y, z coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
-          z = (iz-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
+          z = (iz)*deltah
           zt(k) = f_(x,y,z)
           ! internal point: build discretization
           !
@@ -643,7 +643,7 @@ contains
       f_ => d_null_func_2d
     end if
 
-    deltah   = done/(idim+2)
+    deltah   = done/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_dpk_* deltah
 
@@ -875,8 +875,8 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,glob_row,idim,idim)
           ! x, y coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
 
           zt(k) = f_(x,y)
           ! internal point: build discretization

samples/advanced/pdegen/amg_s_genpde_mod.F90

@@ -176,7 +176,7 @@ contains
     else
       partition_ = 3
     end if
-    deltah   = sone/(idim+2)
+    deltah   = sone/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_spk_* deltah
 
@@ -412,9 +412,9 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
           ! x, y, z coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
-          z = (iz-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
+          z = (iz)*deltah
           zt(k) = f_(x,y,z)
           ! internal point: build discretization
           !
@@ -643,7 +643,7 @@ contains
       f_ => s_null_func_2d
     end if
 
-    deltah   = sone/(idim+2)
+    deltah   = sone/(idim+1)
     sqdeltah = deltah*deltah
     deltah2  = 2.0_psb_spk_* deltah
 
@@ -875,8 +875,8 @@ contains
           ! compute gridpoint coordinates
           call idx2ijk(ix,iy,glob_row,idim,idim)
           ! x, y coordinates
-          x = (ix-1)*deltah
-          y = (iy-1)*deltah
+          x = (ix)*deltah
+          y = (iy)*deltah
 
           zt(k) = f_(x,y)
           ! internal point: build discretization

samples/simple/pdegen/amg_dpde_mod.f90

@@ -369,9 +369,9 @@ contains
         ! compute gridpoint coordinates
         call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
         ! x, y, z coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
-        z = (iz-1)*deltah
+        x = (ix)*deltah
+        y = (iy)*deltah
+        z = (iz)*deltah
         zt(k) = f_(x,y,z)
         ! internal point: build discretization
         !   

samples/simple/pdegen/amg_spde_mod.f90

@@ -369,9 +369,9 @@ contains
         ! compute gridpoint coordinates
         call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
         ! x, y, z coordinates
-        x = (ix-1)*deltah
-        y = (iy-1)*deltah
-        z = (iz-1)*deltah
+        x = (ix)*deltah
+        y = (iy)*deltah
+        z = (iz)*deltah
         zt(k) = f_(x,y,z)
         ! internal point: build discretization
         !