Added string trimming to improve passage of arguments from C interfaces

.VERSION

@@ -1,4 +1,5 @@
 $Format:%d%n%n$
+# Fall back version, probably last release:
 1.2.1
 
 # AMG4PSBLAS version file.

LICENSE

@@ -7,7 +7,7 @@
     (C) Copyright 2025  Salvatore Filippone
     (C) Copyright 2025  Pasqua D'Ambra
     (C) Copyright 2025  Fabio Durastante
-
+    
     Redistribution and use in source and binary forms, with or without
     modification, are permitted provided that the following conditions
     are met:

amgprec/CMakeLists.txt

@@ -509,6 +509,7 @@ set(AMG_amgprec_source_files
   impl/smoother/amg_s_base_smoother_free.f90
   impl/smoother/amg_d_jac_smoother_clone.f90
   impl/smoother/amg_d_jac_smoother_apply_vect.f90
+  impl/smoother/amg_d_richards_smoother_impl.f90
   impl/smoother/amg_c_as_smoother_clear_data.f90
   impl/smoother/amg_s_poly_smoother_descr.f90
   impl/smoother/amg_z_as_smoother_cseti.f90
@@ -773,6 +774,7 @@ set(AMG_amgprec_source_files
   amg_prec_mod.f90
   amg_d_gs_solver.f90
   amg_d_jac_smoother.f90
+  amg_d_richards_smoother.f90
   amg_z_symdec_aggregator_mod.f90
   amg_s_gs_solver.f90
   amg_ainv_mod.f90

amgprec/Makefile

@@ -8,7 +8,7 @@ FINCLUDES=$(FMFLAG)$(HERE) $(FMFLAG)$(INCDIR) $(PSBLAS_INCLUDES)
 
 
 DMODOBJS=amg_d_prec_type.o \
-	amg_d_inner_mod.o amg_d_ilu_solver.o amg_d_diag_solver.o amg_d_jac_smoother.o amg_d_as_smoother.o \
+	amg_d_inner_mod.o amg_d_ilu_solver.o amg_d_diag_solver.o amg_d_jac_smoother.o amg_d_as_smoother.o amg_d_richards_smoother.o \
         amg_d_poly_smoother.o amg_d_poly_coeff_mod.o\
 	amg_d_umf_solver.o amg_d_slu_solver.o amg_d_sludist_solver.o amg_d_id_solver.o\
 	amg_d_base_solver_mod.o amg_d_base_smoother_mod.o amg_d_onelev_mod.o \
@@ -159,7 +159,7 @@ amg_d_umf_solver.o amg_d_diag_solver.o amg_d_ilu_solver.o amg_d_jac_solver.o: am
 
 #amg_d_ilu_fact_mod.o: amg_base_prec_type.o  amg_d_base_solver_mod.o
 #amg_d_ilu_solver.o amg_d_iluk_fact.o: amg_d_ilu_fact_mod.o
-amg_d_as_smoother.o amg_d_jac_smoother.o: amg_d_base_smoother_mod.o
+amg_d_as_smoother.o amg_d_jac_smoother.o amg_d_richards_smoother.o: amg_d_base_smoother_mod.o
 amg_d_jac_smoother.o: amg_d_diag_solver.o
 amg_dprecinit.o amg_dprecset.o: amg_d_diag_solver.o amg_d_ilu_solver.o \
 		amg_d_umf_solver.o amg_d_as_smoother.o amg_d_jac_smoother.o \

amgprec/amg_base_prec_type.F90

@@ -216,7 +216,8 @@ module amg_base_prec_type
   integer(psb_ipk_), parameter :: amg_l1_gs_    = 7
   integer(psb_ipk_), parameter :: amg_l1_fbgs_  = 8
   integer(psb_ipk_), parameter :: amg_poly_     = 9
-  integer(psb_ipk_), parameter :: amg_max_prec_ = 9
+  integer(psb_ipk_), parameter :: amg_richardson_ = 10
+  integer(psb_ipk_), parameter :: amg_max_prec_ = 10
   !
   ! Constants for pre/post signaling. Now only used internally
   !
@@ -229,7 +230,8 @@ module amg_base_prec_type
   !
   ! Legal values for entry: amg_sub_solve_
   !
-  integer(psb_ipk_), parameter :: amg_slv_delta_     = amg_max_prec_+1
+  ! Keep this fixed so sub-solver numeric IDs remain stable.
+  integer(psb_ipk_), parameter :: amg_slv_delta_     = 10
   integer(psb_ipk_), parameter :: amg_f_none_        = amg_slv_delta_+0
   integer(psb_ipk_), parameter :: amg_diag_scale_    = amg_slv_delta_+1
   integer(psb_ipk_), parameter :: amg_l1_diag_scale_ = amg_slv_delta_+2
@@ -409,7 +411,7 @@ module amg_base_prec_type
        & 'none          ','Jacobi        ',&
        & 'L1-Jacobi     ','none          ','none          ',&
        & 'none          ','none          ','L1-GS         ',&
-       & 'L1-FBGS       ','Polynomial    ','none          ','Point Jacobi  ',&
+      & 'L1-FBGS       ','Polynomial    ', 'Richards      ','Point Jacobi  ',&
        & 'L1-Jacobi     ','Gauss-Seidel  ','ILU(n)        ',&
        & 'MILU(n)       ','ILU(t,n)      ',&
        & 'SuperLU       ','UMFPACK LU    ',&
@@ -575,6 +577,8 @@ contains
       val = amg_as_
     case('POLY')
       val = amg_poly_
+    case('RICHARDSON','RICHARDS')
+      val = amg_richardson_
     case('CHEB_4')
       val = amg_cheb_4_
     case('CHEB_4_OPT')
@@ -1202,6 +1206,8 @@ contains
       pr_to_str='BJAC'
     case(amg_as_)
       pr_to_str='AS'
+    case(amg_richardson_)
+      pr_to_str='RICHARDS'
     end select
 
   end function pr_to_str

amgprec/amg_c_jac_smoother.f90

@@ -345,6 +345,7 @@ 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

@@ -671,10 +671,6 @@ 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_config.h.in

@@ -17,6 +17,8 @@
 @CHAVEMUMPS@
 @CHAVEMUMPSMODULES@
 @CHAVEMUMPSINCLUDES@
+@CHAVEMUMPSVERSION@
+@CHAVEMUMPSVERSIONSTRING@
 @CXXMATCHBOXBIT@
 
 

amgprec/amg_d_jac_smoother.f90

@@ -345,6 +345,7 @@ 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_mumps_solver.F90

@@ -420,7 +420,7 @@ subroutine d_mumps_solver_cseti(sv,what,val,info,idx)
   info = psb_success_
   call psb_erractionsave(err_act)
 
-  select case(psb_toupper(what))
+  select case(psb_toupper(trim(what)))
 #if defined(AMG_HAVE_MUMPS)
   case('MUMPS_LOC_GLOB')
     sv%ipar(1) = val
@@ -466,7 +466,7 @@ subroutine d_mumps_solver_csetr(sv,what,val,info,idx)
   info = psb_success_
   call psb_erractionsave(err_act)
 
-  select case(psb_toupper(what))
+  select case(psb_toupper(trim(what)))
 #if defined(AMG_HAVE_MUMPS)
   case('MUMPS_RPAR_ENTRY')
     if(present(idx)) then 

amgprec/amg_d_poly_smoother.f90

@@ -279,6 +279,7 @@ 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

@@ -671,10 +671,6 @@ 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_richards_smoother.f90

@@ -0,0 +1,365 @@
+!
+!
+!                             AMG4PSBLAS version 1.0
+!    Algebraic Multigrid Package
+!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Daniela di Serafino
+!
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the AMG4PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific prior written permission.
+!
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!
+!
+!
+! File: amg_d_richards_smoother.f90
+!
+! Module: amg_d_richards_smoother
+!
+!  This module defines:
+!    the amg_d_richards_smoother_type data structure containing the
+!    smoother for a preconditioned Richards iteration.
+!  The smoother applies the iterative method:
+!    x_{k+1} = x_k + omega * P^{-1} * (b - A * x_k)
+!  where P is a preconditioner (the solver component) acting on the global matrix A.
+!  This allows using distributed solvers like MUMPS on the full system.
+!
+module amg_d_richards_smoother
+
+  use amg_d_base_smoother_mod
+
+  type, extends(amg_d_base_smoother_type) :: amg_d_richards_smoother_type
+    ! The local solver component is inherited from the
+    ! parent type, but acts on the global matrix.
+    !    class(amg_d_base_solver_type), allocatable :: sv
+    !
+    type(psb_dspmat_type), pointer  :: pa => null()
+    integer(psb_lpk_)               :: global_nnz_tot
+    logical                         :: checkres
+    logical                         :: printres
+    integer(psb_ipk_)               :: checkiter
+    integer(psb_ipk_)               :: printiter
+    real(psb_dpk_)                  :: tol
+    real(psb_dpk_)                  :: omega
+  contains
+    procedure, pass(sm) :: apply_v => amg_d_richards_smoother_apply_vect
+    procedure, pass(sm) :: apply_a => amg_d_richards_smoother_apply
+    procedure, pass(sm) :: dump    => amg_d_richards_smoother_dmp
+    procedure, pass(sm) :: build   => amg_d_richards_smoother_bld
+    procedure, pass(sm) :: cnv     => amg_d_richards_smoother_cnv
+    procedure, pass(sm) :: clone   => amg_d_richards_smoother_clone
+    procedure, pass(sm) :: clone_settings => amg_d_richards_smoother_clone_settings
+    procedure, pass(sm) :: clear_data     => amg_d_richards_smoother_clear_data
+    procedure, pass(sm) :: free    => d_richards_smoother_free
+    procedure, pass(sm) :: cseti   => amg_d_richards_smoother_cseti
+    procedure, pass(sm) :: csetc   => amg_d_richards_smoother_csetc
+    procedure, pass(sm) :: csetr   => amg_d_richards_smoother_csetr
+    procedure, pass(sm) :: descr   => amg_d_richards_smoother_descr
+    procedure, pass(sm) :: sizeof  => d_richards_smoother_sizeof
+    procedure, pass(sm) :: default => d_richards_smoother_default
+    procedure, pass(sm) :: get_nzeros => d_richards_smoother_get_nzeros
+    procedure, pass(sm) :: get_wrksz => d_richards_smoother_get_wrksize
+    procedure, nopass   :: get_fmt    => d_richards_smoother_get_fmt
+    procedure, nopass   :: get_id     => d_richards_smoother_get_id
+  end type amg_d_richards_smoother_type
+
+  private :: d_richards_smoother_free, &
+       & d_richards_smoother_sizeof, d_richards_smoother_get_nzeros, &
+       & d_richards_smoother_get_fmt, d_richards_smoother_get_id, &
+       & d_richards_smoother_get_wrksize
+
+  interface
+    subroutine amg_d_richards_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
+         & sweeps,work,wv,info,init,initu)
+      import :: psb_desc_type, amg_d_richards_smoother_type, psb_d_vect_type, psb_dpk_, &
+           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
+           & psb_ipk_
+
+      type(psb_desc_type), intent(in)                 :: desc_data
+      class(amg_d_richards_smoother_type), intent(inout) :: sm
+      type(psb_d_vect_type),intent(inout)           :: x
+      type(psb_d_vect_type),intent(inout)           :: y
+      real(psb_dpk_),intent(in)                      :: alpha,beta
+      character(len=1),intent(in)                     :: trans
+      integer(psb_ipk_), intent(in)                   :: sweeps
+      real(psb_dpk_),target, intent(inout)           :: work(:)
+      type(psb_d_vect_type),intent(inout)           :: wv(:)
+      integer(psb_ipk_), intent(out)                  :: info
+      character, intent(in), optional                :: init
+      type(psb_d_vect_type),intent(inout), optional   :: initu
+    end subroutine amg_d_richards_smoother_apply_vect
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
+         & sweeps,work,info,init,initu)
+      import :: psb_desc_type, amg_d_richards_smoother_type, psb_d_vect_type, psb_dpk_, &
+           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, &
+           & psb_ipk_
+      type(psb_desc_type), intent(in)      :: desc_data
+      class(amg_d_richards_smoother_type), intent(inout) :: sm
+      real(psb_dpk_),intent(inout)         :: x(:)
+      real(psb_dpk_),intent(inout)         :: y(:)
+      real(psb_dpk_),intent(in)            :: alpha,beta
+      character(len=1),intent(in)           :: trans
+      integer(psb_ipk_), intent(in)         :: sweeps
+      real(psb_dpk_),target, intent(inout) :: work(:)
+      integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      real(psb_dpk_),intent(inout), optional :: initu(:)
+    end subroutine amg_d_richards_smoother_apply
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
+      import :: psb_desc_type, amg_d_richards_smoother_type, psb_d_vect_type, psb_dpk_, &
+           & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
+           & psb_ipk_, psb_i_base_vect_type
+      type(psb_dspmat_type), intent(inout), target        :: a
+      Type(psb_desc_type), Intent(inout)                  :: desc_a
+      class(amg_d_richards_smoother_type), intent(inout)       :: sm
+      integer(psb_ipk_), intent(out)                      :: info
+      class(psb_d_base_sparse_mat), intent(in), optional :: amold
+      class(psb_d_base_vect_type), intent(in), optional  :: vmold
+      class(psb_i_base_vect_type), intent(in), optional  :: imold
+    end subroutine amg_d_richards_smoother_bld
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_cnv(sm,info,amold,vmold,imold)
+      import :: amg_d_richards_smoother_type, psb_dpk_, &
+           & psb_d_base_sparse_mat, psb_d_base_vect_type,&
+           & psb_ipk_, psb_i_base_vect_type
+      class(amg_d_richards_smoother_type), intent(inout)       :: sm
+      integer(psb_ipk_), intent(out)                      :: info
+      class(psb_d_base_sparse_mat), intent(in), optional :: amold
+      class(psb_d_base_vect_type), intent(in), optional  :: vmold
+      class(psb_i_base_vect_type), intent(in), optional  :: imold
+    end subroutine amg_d_richards_smoother_cnv
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num)
+      import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
+           & psb_dpk_, amg_d_richards_smoother_type, psb_epk_, psb_desc_type, &
+           & psb_ipk_
+      implicit none
+      class(amg_d_richards_smoother_type), intent(in) :: sm
+      type(psb_desc_type), intent(in)               :: desc
+      integer(psb_ipk_), intent(in)               :: level
+      integer(psb_ipk_), intent(out)              :: info
+      character(len=*), intent(in), optional :: prefix, head
+      logical, optional, intent(in)    :: smoother, solver, global_num
+    end subroutine amg_d_richards_smoother_dmp
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_clone(sm,smout,info)
+      import :: amg_d_richards_smoother_type, psb_dpk_, &
+           & amg_d_base_smoother_type, psb_ipk_
+      class(amg_d_richards_smoother_type), intent(inout)               :: sm
+      class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout
+      integer(psb_ipk_), intent(out)                :: info
+    end subroutine amg_d_richards_smoother_clone
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_clone_settings(sm,smout,info)
+      import :: amg_d_richards_smoother_type, psb_dpk_, &
+           & amg_d_base_smoother_type, psb_ipk_
+      class(amg_d_richards_smoother_type), intent(inout)  :: sm
+      class(amg_d_base_smoother_type), intent(inout) :: smout
+      integer(psb_ipk_), intent(out)                :: info
+    end subroutine amg_d_richards_smoother_clone_settings
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_clear_data(sm,info)
+      import :: amg_d_richards_smoother_type, psb_dpk_, &
+           & amg_d_base_smoother_type, psb_ipk_
+      class(amg_d_richards_smoother_type), intent(inout)               :: sm
+      integer(psb_ipk_), intent(out)                :: info
+    end subroutine amg_d_richards_smoother_clear_data
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_descr(sm,info,iout,coarse,prefix)
+      import :: amg_d_richards_smoother_type, psb_ipk_
+      class(amg_d_richards_smoother_type), intent(in) :: sm
+      integer(psb_ipk_), intent(out)               :: info
+      integer(psb_ipk_), intent(in), optional      :: iout
+      logical, intent(in), optional                :: coarse
+      character(len=*), intent(in), optional       :: prefix
+    end subroutine amg_d_richards_smoother_descr
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_cseti(sm,what,val,info,idx)
+      import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
+           & psb_dpk_, amg_d_richards_smoother_type, psb_epk_, psb_desc_type, psb_ipk_
+      implicit none
+      class(amg_d_richards_smoother_type), intent(inout) :: sm
+      character(len=*), intent(in)                   :: what
+      integer(psb_ipk_), intent(in)                  :: val
+      integer(psb_ipk_), intent(out)                 :: info
+      integer(psb_ipk_), intent(in), optional        :: idx
+    end subroutine amg_d_richards_smoother_cseti
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_csetc(sm,what,val,info,idx)
+      import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
+           & psb_dpk_, amg_d_richards_smoother_type, psb_epk_, psb_desc_type, psb_ipk_
+      implicit none
+      class(amg_d_richards_smoother_type), intent(inout) :: sm
+      character(len=*), intent(in)                   :: what
+      character(len=*), intent(in)                   :: val
+      integer(psb_ipk_), intent(out)                 :: info
+      integer(psb_ipk_), intent(in), optional        :: idx
+    end subroutine amg_d_richards_smoother_csetc
+  end interface
+
+  interface
+    subroutine amg_d_richards_smoother_csetr(sm,what,val,info,idx)
+      import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
+           & psb_dpk_, amg_d_richards_smoother_type, psb_epk_, psb_desc_type, psb_ipk_
+      implicit none
+      class(amg_d_richards_smoother_type), intent(inout) :: sm
+      character(len=*), intent(in)                   :: what
+      real(psb_dpk_), intent(in)                   :: val
+      integer(psb_ipk_), intent(out)                 :: info
+      integer(psb_ipk_), intent(in), optional        :: idx
+    end subroutine amg_d_richards_smoother_csetr
+  end interface
+
+contains
+
+  subroutine d_richards_smoother_free(sm,info)
+    Implicit None
+
+    ! Arguments
+    class(amg_d_richards_smoother_type), intent(inout) :: sm
+    integer(psb_ipk_), intent(out)                  :: info
+    integer(psb_ipk_) :: err_act
+    character(len=20)  :: name='d_richards_smoother_free'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    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)
+        goto 9999
+      end if
+    end if
+    sm%pa => null()
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 call psb_error_handler(err_act)
+    return
+  end subroutine d_richards_smoother_free
+
+  function d_richards_smoother_sizeof(sm) result(val)
+    implicit none
+    ! Arguments
+    class(amg_d_richards_smoother_type), intent(in) :: sm
+    integer(psb_epk_) :: val
+    integer(psb_ipk_)        :: i
+
+    val = psb_sizeof_lp
+    if (allocated(sm%sv)) val = val + sm%sv%sizeof()
+
+    return
+  end function d_richards_smoother_sizeof
+
+  subroutine d_richards_smoother_default(sm)
+    Implicit None
+
+    ! Arguments
+    class(amg_d_richards_smoother_type), intent(inout) :: sm
+
+    !
+    ! Default: Richards iteration with omega=1.0 and no residual check
+    !
+    sm%checkres = .false.
+    sm%printres = .false.
+    sm%checkiter = -1
+    sm%printiter = -1
+    sm%tol = 0
+    sm%omega = 1.0d0
+
+    if (allocated(sm%sv)) then
+      call sm%sv%default()
+    end if
+
+    return
+  end subroutine d_richards_smoother_default
+
+  function d_richards_smoother_get_nzeros(sm) result(val)
+    implicit none
+    ! Arguments
+    class(amg_d_richards_smoother_type), intent(in) :: sm
+    integer(psb_epk_) :: val
+    integer(psb_ipk_)        :: i
+
+    val = 0
+    if (allocated(sm%sv)) val = val + sm%sv%get_nzeros()
+
+    return
+  end function d_richards_smoother_get_nzeros
+
+  function d_richards_smoother_get_wrksize(sm) result(val)
+    implicit none
+    class(amg_d_richards_smoother_type), intent(inout) :: sm
+    integer(psb_ipk_)  :: val
+
+    ! apply_vect uses tx, ty, tz mapped to wv(1:3)
+    val = 3
+    if (allocated(sm%sv)) val = val + sm%sv%get_wrksz()
+
+  end function d_richards_smoother_get_wrksize
+
+  function d_richards_smoother_get_fmt() result(val)
+    implicit none
+    character(len=32)  :: val
+
+    val = "Richards smoother"
+  end function d_richards_smoother_get_fmt
+
+  function d_richards_smoother_get_id() result(val)
+    implicit none
+    integer(psb_ipk_)  :: val
+
+    val = amg_richardson_
+  end function d_richards_smoother_get_id
+
+end module amg_d_richards_smoother

amgprec/amg_s_jac_smoother.f90

@@ -345,6 +345,7 @@ 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,6 +279,7 @@ 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

@@ -671,10 +671,6 @@ 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,6 +345,7 @@ 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

@@ -671,10 +671,6 @@ 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/level/amg_d_base_onelev_csetc.F90

@@ -46,6 +46,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
   use amg_d_poly_smoother
   use amg_d_jac_smoother
   use amg_d_as_smoother
+  use amg_d_richards_smoother
   use amg_d_diag_solver
   use amg_d_l1_diag_solver
   use amg_d_jac_solver
@@ -97,6 +98,7 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
   type(amg_d_invk_solver_type)     ::  amg_d_invk_solver_mold
   type(amg_d_invt_solver_type)     ::  amg_d_invt_solver_mold
   type(amg_d_poly_smoother_type)   ::  amg_d_poly_smoother_mold
+  type(amg_d_richards_smoother_type) :: amg_d_richards_smoother_mold
 #if defined(AMG_HAVE_UMF)
   type(amg_d_umf_solver_type)     ::  amg_d_umf_solver_mold
 #endif
@@ -161,6 +163,11 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
     case ('POLY')
       call lv%set(amg_d_poly_smoother_mold,info,pos=pos)
       if (info == 0) call lv%set(amg_d_l1_diag_solver_mold,info,pos=pos)
+
+    case ('RICHARDS','RICHARDSON')
+      call lv%set(amg_d_richards_smoother_mold,info,pos=pos)
+      if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos)
+
     case ('GS','FWGS')
       call lv%set(amg_d_jac_smoother_mold,info,pos='pre')
       if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre')

amgprec/impl/level/amg_d_base_onelev_cseti.F90

@@ -45,6 +45,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx)
   use amg_d_parmatch_aggregator_mod
   use amg_d_jac_smoother
   use amg_d_as_smoother
+  use amg_d_richards_smoother
   use amg_d_diag_solver
   use amg_d_l1_diag_solver
   use amg_d_ilu_solver
@@ -79,6 +80,7 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx)
   type(amg_d_jac_smoother_type)    ::  amg_d_jac_smoother_mold
   type(amg_d_l1_jac_smoother_type) ::  amg_d_l1_jac_smoother_mold
   type(amg_d_as_smoother_type)     ::  amg_d_as_smoother_mold
+  type(amg_d_richards_smoother_type) :: amg_d_richards_smoother_mold
   type(amg_d_diag_solver_type)     ::  amg_d_diag_solver_mold
   type(amg_d_l1_diag_solver_type)  ::  amg_d_l1_diag_solver_mold
   type(amg_d_ilu_solver_type)      ::  amg_d_ilu_solver_mold
@@ -141,6 +143,10 @@ subroutine amg_d_base_onelev_cseti(lv,what,val,info,pos,idx)
       call lv%set(amg_d_as_smoother_mold,info,pos=pos)
       if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos)
 
+    case (amg_richardson_)
+      call lv%set(amg_d_richards_smoother_mold,info,pos=pos)
+      if (info == 0) call lv%set(amg_d_ilu_solver_mold,info,pos=pos)
+
     case (amg_fbgs_)
       call lv%set(amg_d_jac_smoother_mold,info,pos='pre')
       if (info == 0) call lv%set(amg_d_gs_solver_mold,info,pos='pre')

amgprec/impl/smoother/Makefile

@@ -94,6 +94,7 @@ amg_d_jac_smoother_cnv.o \
 amg_d_jac_smoother_csetc.o \
 amg_d_jac_smoother_cseti.o \
 amg_d_jac_smoother_csetr.o \
+amg_d_richards_smoother_impl.o \
 amg_d_l1_jac_smoother_bld.o \
 amg_d_l1_jac_smoother_descr.o \
 amg_d_l1_jac_smoother_clone.o \

amgprec/impl/smoother/amg_c_as_smoother_free.f90

@@ -53,6 +53,7 @@ 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,6 +52,7 @@ 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,6 +53,7 @@ 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,6 +52,7 @@ 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_d_richards_smoother_impl.f90

@@ -0,0 +1,644 @@
+!
+!
+!                             AMG4PSBLAS version 1.0
+!    Algebraic Multigrid Package
+!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
+!
+!
+subroutine amg_d_richards_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
+     & sweeps,work,wv,info,init,initu)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_apply_vect
+  implicit none
+
+  type(psb_desc_type), intent(in)                    :: desc_data
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  type(psb_d_vect_type), intent(inout)               :: x
+  type(psb_d_vect_type), intent(inout)               :: y
+  real(psb_dpk_), intent(in)                         :: alpha, beta
+  character(len=1), intent(in)                       :: trans
+  integer(psb_ipk_), intent(in)                      :: sweeps
+  real(psb_dpk_), target, intent(inout)              :: work(:)
+  type(psb_d_vect_type), intent(inout)               :: wv(:)
+  integer(psb_ipk_), intent(out)                     :: info
+  character, intent(in), optional                    :: init
+  type(psb_d_vect_type), intent(inout), optional     :: initu
+
+  integer(psb_ipk_)   :: n_col, err_act, i
+  character           :: trans_, init_
+  real(psb_dpk_), pointer :: aux(:)
+  character(len=32)   :: name='d_richards_smoother_apply_v'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  trans_ = psb_toupper(trans)
+  if ((trans_ /= 'N').and.(trans_ /= 'T').and.(trans_ /= 'C')) then
+    call psb_errpush(psb_err_iarg_invalid_i_,name)
+    goto 9999
+  end if
+
+  if (sweeps < 0) then
+    info = psb_err_iarg_neg_
+    call psb_errpush(info,name,i_err=(/itwo,sweeps,izero,izero,izero/))
+    goto 9999
+  end if
+
+  if (.not.allocated(sm%sv)) then
+    info = psb_err_invalid_cd_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  if (.not.associated(sm%pa)) then
+    info = psb_err_invalid_cd_state_
+    call psb_errpush(info,name,a_err='matrix pointer not associated')
+    goto 9999
+  end if
+
+  if (size(wv) < 3) then
+    info = psb_err_internal_error_
+    call psb_errpush(info,name,a_err='workspace vectors too small')
+    goto 9999
+  end if
+
+  init_ = 'Z'
+  if (present(init)) init_ = psb_toupper(init)
+
+  n_col = desc_data%get_local_cols()
+  if (4*n_col <= size(work)) then
+    aux => work(:)
+  else
+    allocate(aux(4*n_col),stat=info)
+    if (info /= psb_success_) then
+      info = psb_err_alloc_request_
+      call psb_errpush(info,name,i_err=(/4*n_col,izero,izero,izero,izero/),&
+           & a_err='real(psb_dpk_)')
+      goto 9999
+    end if
+  end if
+
+  associate(tx => wv(1), ty => wv(2), tz => wv(3))
+
+    select case(init_)
+    case('Z')
+      call psb_geaxpby(dzero,x,dzero,ty,desc_data,info)
+    case('Y')
+      call psb_geaxpby(done,y,dzero,ty,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='missing initu to smoother_apply_vect')
+        goto 9999
+      end if
+      call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='wrong init to smoother_apply_vect')
+      goto 9999
+    end select
+
+    do i = 1, sweeps
+      call psb_geaxpby(done,x,dzero,tx,desc_data,info)
+      call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
+      if (info /= psb_success_) exit
+
+      call sm%sv%apply(done,tx,dzero,tz,desc_data,trans_,aux,wv(4:),info,init='Y')
+      if (info /= psb_success_) exit
+
+      call psb_geaxpby(sm%omega,tz,done,ty,desc_data,info)
+      if (info /= psb_success_) exit
+    end do
+
+    if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
+
+  end associate
+
+  if (.not.(4*n_col <= size(work))) deallocate(aux)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_apply_vect
+
+
+subroutine amg_d_richards_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
+     & sweeps,work,info,init,initu)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_apply
+  implicit none
+
+  type(psb_desc_type), intent(in)                    :: desc_data
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  real(psb_dpk_), intent(inout)                      :: x(:)
+  real(psb_dpk_), intent(inout)                      :: y(:)
+  real(psb_dpk_), intent(in)                         :: alpha, beta
+  character(len=1), intent(in)                       :: trans
+  integer(psb_ipk_), intent(in)                      :: sweeps
+  real(psb_dpk_), target, intent(inout)              :: work(:)
+  integer(psb_ipk_), intent(out)                     :: info
+  character, intent(in), optional                    :: init
+  real(psb_dpk_), intent(inout), optional            :: initu(:)
+
+  integer(psb_ipk_)      :: n_col, err_act, i
+  character              :: trans_, init_
+  real(psb_dpk_), pointer :: aux(:)
+  real(psb_dpk_), allocatable :: tx(:), ty(:), tz(:)
+  character(len=30)      :: name='d_richards_smoother_apply'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  trans_ = psb_toupper(trans)
+  if ((trans_ /= 'N').and.(trans_ /= 'T').and.(trans_ /= 'C')) then
+    call psb_errpush(psb_err_iarg_invalid_i_,name)
+    goto 9999
+  end if
+
+  if (sweeps < 0) then
+    info = psb_err_iarg_neg_
+    call psb_errpush(info,name,i_err=(/itwo,sweeps,izero,izero,izero/))
+    goto 9999
+  end if
+
+  if (.not.allocated(sm%sv)) then
+    info = psb_err_invalid_cd_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  if (.not.associated(sm%pa)) then
+    info = psb_err_invalid_cd_state_
+    call psb_errpush(info,name,a_err='matrix pointer not associated')
+    goto 9999
+  end if
+
+  n_col = desc_data%get_local_cols()
+  if (4*n_col <= size(work)) then
+    aux => work(:)
+  else
+    allocate(aux(4*n_col),stat=info)
+    if (info /= psb_success_) then
+      info = psb_err_alloc_request_
+      call psb_errpush(info,name,i_err=(/4*n_col,izero,izero,izero,izero/),&
+           & a_err='real(psb_dpk_)')
+      goto 9999
+    end if
+  end if
+
+  allocate(tx(size(x)),ty(size(y)),tz(size(y)),stat=info)
+  if (info /= 0) then
+    info = psb_err_alloc_request_
+    call psb_errpush(info,name,a_err='alloc tx/ty/tz')
+    goto 9999
+  end if
+
+  init_ = 'Z'
+  if (present(init)) init_ = psb_toupper(init)
+
+  select case(init_)
+  case('Z')
+    ty(:) = dzero
+  case('Y')
+    ty(:) = y(:)
+  case('U')
+    if (.not.present(initu)) then
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    ty(:) = initu(:)
+  case default
+    info = psb_err_internal_error_
+    call psb_errpush(info,name,a_err='wrong init to smoother_apply')
+    goto 9999
+  end select
+
+  do i = 1, sweeps
+    tx(:) = x(:)
+    call psb_spmm(-done,sm%pa,ty,done,tx,desc_data,info,work=aux,trans=trans_)
+    if (info /= psb_success_) exit
+
+    call sm%sv%apply(done,tx,dzero,tz,desc_data,trans_,aux,info,init='Y')
+    if (info /= psb_success_) exit
+
+    ty(:) = ty(:) + sm%omega*tz(:)
+  end do
+
+  if (info == psb_success_) y(:) = beta*y(:) + alpha*ty(:)
+
+  deallocate(tx,ty,tz,stat=info)
+  if (.not.(4*n_col <= size(work))) deallocate(aux)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_apply
+
+
+subroutine amg_d_richards_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_bld
+  implicit none
+
+  type(psb_dspmat_type), intent(inout), target         :: a
+  type(psb_desc_type), intent(inout)                   :: desc_a
+  class(amg_d_richards_smoother_type), intent(inout)   :: sm
+  integer(psb_ipk_), intent(out)                       :: info
+  class(psb_d_base_sparse_mat), intent(in), optional   :: amold
+  class(psb_d_base_vect_type), intent(in), optional    :: vmold
+  class(psb_i_base_vect_type), intent(in), optional    :: imold
+
+  integer(psb_ipk_) :: err_act
+  character(len=28) :: name='d_richards_smoother_bld'
+
+  info = psb_success_
+  call psb_erractionsave(err_act)
+
+  sm%pa => a
+  if (sm%omega == dzero) sm%omega = done
+
+  if (.not.allocated(sm%sv)) then
+    info = psb_err_invalid_cd_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call sm%sv%build(a,desc_a,info,amold=amold,vmold=vmold,imold=imold)
+  if (info /= psb_success_) then
+    call psb_errpush(psb_err_from_subroutine_,name,a_err='solver build')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_bld
+
+
+subroutine amg_d_richards_smoother_cnv(sm,info,amold,vmold,imold)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_cnv
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout)   :: sm
+  integer(psb_ipk_), intent(out)                       :: info
+  class(psb_d_base_sparse_mat), intent(in), optional   :: amold
+  class(psb_d_base_vect_type), intent(in), optional    :: vmold
+  class(psb_i_base_vect_type), intent(in), optional    :: imold
+
+  integer(psb_ipk_) :: err_act
+  character(len=28) :: name='d_richards_smoother_cnv'
+
+  info = psb_success_
+  call psb_erractionsave(err_act)
+
+  if (allocated(sm%sv)) call sm%sv%cnv(info,amold=amold,vmold=vmold,imold=imold)
+  if (info /= psb_success_) then
+    call psb_errpush(psb_err_from_subroutine_,name,a_err='solver cnv')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_cnv
+
+
+subroutine amg_d_richards_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_dmp
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(in)      :: sm
+  type(psb_desc_type), intent(in)                      :: desc
+  integer(psb_ipk_), intent(in)                        :: level
+  integer(psb_ipk_), intent(out)                       :: info
+  character(len=*), intent(in), optional               :: prefix, head
+  logical, intent(in), optional                        :: smoother, solver, global_num
+
+  info = psb_success_
+  if (allocated(sm%sv)) call sm%sv%dump(desc,level,info,solver=solver,prefix=prefix,global_num=global_num)
+
+end subroutine amg_d_richards_smoother_dmp
+
+
+subroutine amg_d_richards_smoother_clone(sm,smout,info)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_clone
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout)               :: sm
+  class(amg_d_base_smoother_type), allocatable, intent(inout)      :: smout
+  integer(psb_ipk_), intent(out)                                   :: info
+
+  integer(psb_ipk_) :: err_act
+
+  info = psb_success_
+  call psb_erractionsave(err_act)
+
+  if (allocated(smout)) then
+    call smout%free(info)
+    if (info == psb_success_) deallocate(smout, stat=info)
+  end if
+  if (info == psb_success_) allocate(amg_d_richards_smoother_type :: smout, stat=info)
+  if (info /= 0) then
+    info = psb_err_alloc_dealloc_
+    goto 9999
+  end if
+
+  select type(smo => smout)
+  type is (amg_d_richards_smoother_type)
+    smo%global_nnz_tot = sm%global_nnz_tot
+    smo%checkres = sm%checkres
+    smo%printres = sm%printres
+    smo%checkiter = sm%checkiter
+    smo%printiter = sm%printiter
+    smo%tol = sm%tol
+    smo%omega = sm%omega
+    smo%pa => sm%pa
+    if (allocated(sm%sv)) then
+      allocate(smo%sv,mold=sm%sv,stat=info)
+      if (info == psb_success_) call sm%sv%clone(smo%sv,info)
+    end if
+  class default
+    info = psb_err_internal_error_
+  end select
+
+  if (info /= 0) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_clone
+
+
+subroutine amg_d_richards_smoother_clone_settings(sm,smout,info)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_clone_settings
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  class(amg_d_base_smoother_type), intent(inout)     :: smout
+  integer(psb_ipk_), intent(out)                     :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=36) :: name='d_richards_smoother_clone_settings'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  select type(smout)
+  class is(amg_d_richards_smoother_type)
+    smout%pa => null()
+    smout%global_nnz_tot = 0
+    smout%checkres = sm%checkres
+    smout%printres = sm%printres
+    smout%checkiter = sm%checkiter
+    smout%printiter = sm%printiter
+    smout%tol = sm%tol
+    smout%omega = sm%omega
+
+    if (allocated(smout%sv)) then
+      if (.not.same_type_as(sm%sv,smout%sv)) then
+        call smout%sv%free(info)
+        if (info == 0) deallocate(smout%sv,stat=info)
+      end if
+    end if
+
+    if (info == 0) then
+      if (allocated(smout%sv)) then
+        if (same_type_as(sm%sv,smout%sv)) then
+          call sm%sv%clone_settings(smout%sv,info)
+        else
+          info = psb_err_internal_error_
+        end if
+      else
+        allocate(smout%sv,mold=sm%sv,stat=info)
+        if (info == 0) call sm%sv%clone_settings(smout%sv,info)
+      end if
+    end if
+  class default
+    info = psb_err_internal_error_
+  end select
+
+  if (info /= 0) then
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_clone_settings
+
+
+subroutine amg_d_richards_smoother_clear_data(sm,info)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_clear_data
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  integer(psb_ipk_), intent(out)                     :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=31) :: name='d_richards_smoother_clear_data'
+
+  call psb_erractionsave(err_act)
+
+  info = psb_success_
+  sm%pa => null()
+  sm%global_nnz_tot = 0
+  if ((info == 0).and.allocated(sm%sv)) call sm%sv%clear_data(info)
+  if (info /= 0) then
+    info = psb_err_internal_error_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_clear_data
+
+
+subroutine amg_d_richards_smoother_descr(sm,info,iout,coarse,prefix)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_descr
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(in) :: sm
+  integer(psb_ipk_), intent(out)                  :: info
+  integer(psb_ipk_), intent(in), optional         :: iout
+  logical, intent(in), optional                   :: coarse
+  character(len=*), intent(in), optional          :: prefix
+
+  integer(psb_ipk_)   :: iout_
+  character(1024)     :: prefix_
+
+  info = psb_success_
+  iout_ = psb_out_unit
+  if (present(iout)) iout_ = iout
+  prefix_ = ''
+  if (present(prefix)) prefix_ = prefix
+
+  write(iout_,*) trim(prefix_), '  Richardson smoother'
+  write(iout_,*) trim(prefix_), '       Relaxation omega: ', sm%omega
+  if (allocated(sm%sv)) then
+    write(iout_,*) trim(prefix_), '       Preconditioner details:'
+    call sm%sv%descr(info,iout_,coarse=coarse,prefix=prefix)
+  end if
+
+end subroutine amg_d_richards_smoother_descr
+
+
+subroutine amg_d_richards_smoother_cseti(sm,what,val,info,idx)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_cseti
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  character(len=*), intent(in)                       :: what
+  integer(psb_ipk_), intent(in)                      :: val
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_ipk_), intent(in), optional            :: idx
+
+  integer(psb_ipk_) :: err_act
+
+  info = psb_success_
+  call psb_erractionsave(err_act)
+
+  select case(psb_toupper(what))
+  case('SMOOTHER_RESIDUAL')
+    sm%checkiter = val
+  case('SMOOTHER_ITRACE')
+    sm%printiter = val
+  case default
+    call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx)
+  end select
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_cseti
+
+
+subroutine amg_d_richards_smoother_csetc(sm,what,val,info,idx)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_csetc
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  character(len=*), intent(in)                       :: what
+  character(len=*), intent(in)                       :: val
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_ipk_), intent(in), optional            :: idx
+
+  integer(psb_ipk_) :: err_act
+  character(len=30) :: name='d_richards_smoother_csetc'
+
+  info = psb_success_
+  call psb_erractionsave(err_act)
+
+  select case(psb_toupper(trim(what)))
+  case('SMOOTHER_STOP')
+    select case(psb_toupper(trim(val)))
+    case('T','TRUE')
+      sm%checkres = .true.
+    case('F','FALSE')
+      sm%checkres = .false.
+    end select
+  case('SMOOTHER_TRACE')
+    select case(psb_toupper(trim(val)))
+    case('T','TRUE')
+      sm%printres = .true.
+    case('F','FALSE')
+      sm%printres = .false.
+    end select
+  case default
+    call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx)
+  end select
+
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_csetc
+
+
+subroutine amg_d_richards_smoother_csetr(sm,what,val,info,idx)
+
+  use psb_base_mod
+  use amg_d_richards_smoother, amg_protect_name => amg_d_richards_smoother_csetr
+  implicit none
+
+  class(amg_d_richards_smoother_type), intent(inout) :: sm
+  character(len=*), intent(in)                       :: what
+  real(psb_dpk_), intent(in)                         :: val
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_ipk_), intent(in), optional            :: idx
+
+  integer(psb_ipk_) :: err_act
+
+  info = psb_success_
+  call psb_erractionsave(err_act)
+
+  select case(psb_toupper(what))
+  case('SMOOTHER_STOPTOL')
+    sm%tol = val
+  case('RICHARDSON_OMEGA','SMOOTHER_OMEGA','OMEGA')
+    sm%omega = val
+  case default
+    call sm%amg_d_base_smoother_type%set(what,val,info,idx=idx)
+  end select
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
+
+end subroutine amg_d_richards_smoother_csetr

amgprec/impl/smoother/amg_s_as_smoother_free.f90

@@ -53,6 +53,7 @@ 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,6 +52,7 @@ 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,6 +53,7 @@ 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,6 +52,7 @@ 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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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 +176,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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 +176,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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_mumps_solver_apply.F90

@@ -59,14 +59,18 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
   integer(psb_lpk_)  :: nglob
   integer(psb_epk_)  :: eng
   real(psb_dpk_), allocatable     :: ww(:)
+  real(psb_dpk_), allocatable, target :: rhs_loc(:), sol_loc(:)
   real(psb_dpk_), allocatable, target :: gx(:)
+  integer(psb_lpk_), allocatable  :: gidx(:)
+  integer, allocatable, target    :: irhs_loc(:), isol_loc(:)
   integer(psb_ipk_)   :: i, err_act
   character           :: trans_
   character(len=20)   :: name='d_mumps_solver_apply'
 
   call psb_erractionsave(err_act)
 
-#if defined(AMG_HAVE_MUMPS) 
+#if defined(AMG_HAVE_MUMPS)
+#if defined(AMG_MUMPS_VERSION) && (AMG_MUMPS_VERSION <= 590)
   info = psb_success_
   trans_ = psb_toupper(trans)
   select case(trans_)
@@ -154,6 +158,147 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
   call psb_erractionrestore(err_act)
   return
 
+#else
+
+  ! Snapshot MUMPS branch (e.g. snapshot06052026 and newer snapshots).
+  info = psb_success_
+  trans_ = psb_toupper(trans)
+  select case(trans_)
+  case('N')
+  case('T')
+  case default
+    call psb_errpush(psb_err_iarg_invalid_i_,name)
+    goto 9999
+  end select
+
+  nglob = desc_data%get_global_rows()
+  n_row = desc_data%get_local_rows()
+  n_col = desc_data%get_local_cols()
+
+  if (sv%ipar(1) == amg_local_solver_ ) then
+    gx = x
+    sv%id%icntl(20) = 0
+    sv%id%icntl(21) = 0
+  else if (sv%ipar(1) == amg_global_solver_ ) then
+
+    if (n_col <= size(work)) then
+      ww = work(1:n_col)
+    else
+      allocate(ww(n_col),stat=info)
+      if (info /= psb_success_) then
+        info=psb_err_alloc_request_
+        call psb_errpush(info,name,i_err=(/n_col/),&
+             & a_err='real(psb_dpk_)')
+        goto 9999
+      end if
+    end if
+    allocate(rhs_loc(max(n_row,1)),sol_loc(max(n_row,1)),&
+         & irhs_loc(max(n_row,1)),isol_loc(max(n_row,1)),stat=info)
+    if (info /= psb_success_) then
+      info=psb_err_alloc_request_; eng = max(n_row,1)
+      call psb_errpush(info,name,e_err=(/eng/),&
+           & a_err='real(psb_dpk_)')
+      goto 9999
+    end if
+    ww = 0.0_psb_dpk_
+    sol_loc = 0.0_psb_dpk_
+    rhs_loc = 0.0_psb_dpk_
+    gidx = desc_data%get_global_indices(owned=.true.)
+    if (size(gidx) /= n_row) then
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,&
+           & a_err='Invalid local row distribution in MUMPS')
+      goto 9999
+    end if
+    if ((n_row > 0) .and. ((minval(gidx) < 1_psb_lpk_) .or. &
+         & (maxval(gidx) > int(huge(0),psb_lpk_)))) then
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,&
+           & a_err='Overflow in distributed MUMPS RHS indices')
+      goto 9999
+    end if
+    if (n_row > 0) then
+      rhs_loc(1:n_row) = x(1:n_row)
+      irhs_loc(1:n_row) = int(gidx(1:n_row))
+      isol_loc(1:n_row) = int(gidx(1:n_row))
+    end if
+  else
+      info=psb_err_internal_error_
+      call psb_errpush(info,name,&
+           & a_err='Invalid local/global solver in MUMPS')
+      goto 9999
+  end if
+
+  select case(trans_)
+  case('N')
+    sv%id%icntl(9) = 1
+  case('T')
+    sv%id%icntl(9) = 2
+  case default
+    call psb_errpush(psb_err_internal_error_,&
+         & name,a_err='Invalid TRANS in subsolve')
+    goto 9999
+  end select
+
+  sv%id%nrhs =  1
+  if (sv%ipar(1) == amg_local_solver_ ) then
+    sv%id%nloc_rhs  = 0
+    sv%id%lrhs_loc  = 0
+    sv%id%nsol_loc  = 0
+    sv%id%lsol_loc  = 0
+    nullify(sv%id%rhs_loc)
+    nullify(sv%id%irhs_loc)
+    nullify(sv%id%sol_loc)
+    nullify(sv%id%isol_loc)
+    sv%id%rhs  => gx
+  else
+    nullify(sv%id%rhs)
+    sv%id%icntl(20) = 10
+    sv%id%nloc_rhs  = n_row
+    sv%id%lrhs_loc  = n_row
+    sv%id%rhs_loc   => rhs_loc
+    sv%id%irhs_loc  => irhs_loc
+    sv%id%icntl(21) = 2
+    sv%id%nsol_loc  = n_row
+    sv%id%lsol_loc  = n_row
+    sv%id%sol_loc   => sol_loc
+    sv%id%isol_loc  => isol_loc
+  end if
+  ! Snapshot path: preserve stream settings coming from build,
+  ! only enforce silent mode if the user did not set them explicitly.
+  if (sv%id%icntl(1) == 0) sv%id%icntl(1) = -1
+  if (sv%id%icntl(2) == 0) sv%id%icntl(2) = -1
+  if (sv%id%icntl(3) == 0) sv%id%icntl(3) = -1
+  if (sv%id%icntl(4) == 0) sv%id%icntl(4) = -1
+  sv%id%job = 3
+  call dmumps(sv%id)
+
+  if (sv%ipar(1) == amg_local_solver_ ) then
+    call psb_geaxpby(alpha,gx,beta,y,desc_data,info)
+  else
+    if (n_row > 0) ww(1:n_row) = sol_loc(1:n_row)
+    call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
+  end if
+
+  nullify(sv%id%rhs)
+  nullify(sv%id%rhs_loc)
+  nullify(sv%id%irhs_loc)
+  nullify(sv%id%sol_loc)
+  nullify(sv%id%isol_loc)
+
+  if (info /= psb_success_) then
+    call psb_errpush(psb_err_internal_error_,&
+         & name,a_err='Error in subsolve')
+    goto 9999
+  endif
+
+  if (allocated(ww)) deallocate(ww)
+
+  call psb_erractionrestore(err_act)
+  return
+
+#endif
+
 9999 continue
   call psb_erractionrestore(err_act)
   if (err_act == psb_act_abort_) then

amgprec/impl/solver/amg_d_mumps_solver_bld.F90

@@ -69,7 +69,8 @@ subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   integer(psb_ipk_)  :: np, iam, me, i, err_act, debug_unit, debug_level
   character(len=20)  :: name='d_mumps_solver_bld', ch_err
 
-#if defined(AMG_HAVE_MUMPS) 
+#if defined(AMG_HAVE_MUMPS)
+#if defined(AMG_MUMPS_VERSION) && (AMG_MUMPS_VERSION <= 590)
 
   info=psb_success_
 
@@ -245,6 +246,171 @@ subroutine d_mumps_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
   call psb_erractionrestore(err_act)
   return
 
+#else
+
+  ! Snapshot MUMPS branch (e.g. snapshot06052026 and newer snapshots).
+
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+  ctxt       = desc_a%get_context()
+  call psb_info(ctxt, iam, np)
+  if (sv%ipar(1) == amg_local_solver_ ) then
+    call psb_init(ctxt1,np=1,basectxt=ctxt,ids=(/iam/))
+  else if (sv%ipar(1) == amg_global_solver_ ) then
+    call psb_init(ctxt1,basectxt=ctxt)
+  else
+    info = psb_err_internal_error_
+    call psb_errpush(info,name,&
+         & a_err='Invalid local/global solver in MUMPS')
+    goto 9999
+  end if
+  icomm = ctxt1%get_mpic()
+  sv%local_ctxt = ctxt1
+  call psb_info(ctxt1, me, npr)
+  npc  = 1
+  if (debug_level >= psb_debug_outer_) &
+       & write(debug_unit,*) me,' ',trim(name),' start'
+
+  if(.not.allocated(sv%id)) then
+    allocate(sv%id,stat=info)
+    if (info /= psb_success_) then
+      info=psb_err_alloc_dealloc_
+      call psb_errpush(info,name,a_err='amg_dmumps_default')
+      goto 9999
+    end if
+  end if
+
+  ! Snapshot path: bind directly to the local communicator.
+  sv%id%comm = icomm
+
+  sv%id%job     = -1
+  sv%id%par     =  1
+  if (sv%ipar(3) == 2) then
+    sv%id%sym = 2
+  else
+    sv%id%sym = 0
+  end if
+
+  call dmumps(sv%id)
+
+  if (allocated(sv%icntl)) then
+    do i=1,amg_mumps_icntl_size
+      if (allocated(sv%icntl(i)%item)) then
+        sv%id%icntl(i) = sv%icntl(i)%item
+      end if
+    end do
+  end if
+  if (allocated(sv%rcntl)) then
+    do i=1,amg_mumps_rcntl_size
+      if (allocated(sv%rcntl(i)%item)) sv%id%cntl(i) = sv%rcntl(i)%item
+    end do
+  end if
+  sv%id%icntl(5)=0
+  sv%id%icntl(3)=sv%ipar(2)
+
+  nglob  = desc_a%get_global_rows()
+  nrow_a = a%get_nrows()
+  if (sv%ipar(1) == amg_local_solver_ ) then
+    nglobrec=desc_a%get_local_rows()
+    if (sv%ipar(3) == 2) then
+      call a%triu(c,info,jmax=nrow_a)
+      call c%set_symmetric()
+    else
+      call a%csclip(c,info,jmax=nrow_a)
+    end if
+    call c%cp_to(acoo)
+    nglob = c%get_nrows()
+    if (nglobrec /= nglob) then
+      write(*,*)'WARNING: MUMPS solver does not allow overlap in AS yet. '
+      write(*,*)'A zero-overlap is used instead'
+    end if
+  else
+    call a%cp_to(acoo)
+  end if
+  nza = acoo%get_nzeros()
+
+  if ((sv%ipar(1) == amg_global_solver_) .and. (nza > 0)) then
+#if defined(PSB_IPK4) && defined(PSB_LPK8)
+    if (nglob > huge(1_psb_ipk_)) then
+      write(0,*) iam,' ',trim(name),': Error: overflow of local indices '
+      info=psb_err_internal_error_
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    gia = acoo%ia(1:nza)
+    gja = acoo%ja(1:nza)
+    call psb_loc_to_glob(gia(1:nza), desc_a, info, iact='I')
+    call psb_loc_to_glob(gja(1:nza), desc_a, info, iact='I')
+    acoo%ia(1:nza) = gia(1:nza)
+    acoo%ja(1:nza) = gja(1:nza)
+#else
+    call psb_loc_to_glob(acoo%ja(1:nza), desc_a, info, iact='I')
+    call psb_loc_to_glob(acoo%ia(1:nza), desc_a, info, iact='I')
+#endif
+
+    if (sv%ipar(3) == 2 ) then
+      block
+        integer(psb_ipk_) :: j,nz
+        nz = 0
+        do j=1,nza
+          if (acoo%ja(j) >= acoo%ia(j)) then
+            nz = nz + 1
+            acoo%ia(nz)  = acoo%ia(j)
+            acoo%ja(nz)  = acoo%ja(j)
+            acoo%val(nz) = acoo%val(j)
+          end if
+        end do
+        call acoo%set_nzeros(nz)
+        call acoo%set_triangle()
+        call acoo%set_upper()
+        call acoo%set_symmetric()
+        nza = nz
+      end block
+    end if
+  end if
+
+  ! Snapshot input mode: distributed assembled matrix.
+  ! Each process contributes its own non-overlapping local triplets.
+  sv%id%irn_loc   => acoo%ia
+  sv%id%jcn_loc   => acoo%ja
+  sv%id%a_loc     => acoo%val
+  sv%id%icntl(18) = 3
+  sv%id%n      = nglob
+  sv%id%nnz_loc = acoo%get_nzeros()
+  sv%id%nnz     = acoo%get_nzeros()
+  sv%id%job    = 4
+  if (sv%ipar(1) == amg_global_solver_ ) then
+    call psb_sum(ctxt,sv%id%nnz)
+  end if
+
+  call dmumps(sv%id)
+  info = sv%id%infog(1)
+  if (info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='amg_dmumps_fact '
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  nullify(sv%id%irn)
+  nullify(sv%id%jcn)
+  nullify(sv%id%a)
+
+  call acoo%free()
+  sv%built=.true.
+
+  if (debug_level >= psb_debug_outer_) &
+       & write(debug_unit,*) iam,' ',trim(name),' end'
+
+  call psb_erractionrestore(err_act)
+  return
+
+#endif
+
 9999 continue
   call psb_erractionrestore(err_act)
   if (err_act == psb_act_abort_) then

amgprec/impl/solver/amg_s_diag_solver_bld.f90

@@ -91,7 +91,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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 +176,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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 +176,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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,11 @@ 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
-  if (.not.allocated(sv%dv)) allocate(sv%dv,stat=info) 
+  if (allocated(sv%dv)) then
+    call sv%dv%free(info)
+    deallocate(sv%dv)
+  end if
+  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.c

@@ -11,7 +11,7 @@ amg_c_dprec* amg_c_dprec_new()
 }
 
 
-psb_c_i_t amg_c_dprec_delete(amg_c_dprec* p)
+psb_i_t amg_c_dprec_delete(amg_c_dprec* p)
 {
   int iret;
   iret=amg_c_dprecfree(p);

cbind/amgprec/amg_c_dprec.h

@@ -17,25 +17,37 @@ extern "C" {
   } amg_c_dprec; 
   
   amg_c_dprec* amg_c_dprec_new();
-  psb_c_i_t amg_c_dprec_delete(amg_c_dprec* p);
+  psb_i_t amg_c_dprec_delete(amg_c_dprec* p);
  
-  psb_c_i_t amg_c_dprecinit(psb_c_ctxt cctxt, amg_c_dprec *ph, const char *ptype);
-  psb_c_i_t amg_c_dprecseti(amg_c_dprec *ph, const char *what, psb_c_i_t val);
-  psb_c_i_t amg_c_dprecsetc(amg_c_dprec *ph, const char *what, const char *val);
-  psb_c_i_t amg_c_dprecsetr(amg_c_dprec *ph, const char *what, double val);
-  psb_c_i_t amg_c_dprecbld(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph);
-  psb_c_i_t amg_c_dhierarchy_build(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph);
-  psb_c_i_t amg_c_dsmoothers_build(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph);
-  psb_c_i_t amg_c_dsmoothers_build_opt(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph, const char *afmt, const char *chfmt);
-  psb_c_i_t amg_c_dprecapply(amg_c_dprec *ph, psb_c_dvector *bh, psb_c_dvector *xh, psb_c_descriptor *cdh);
-  psb_c_i_t amg_c_dprecapply_opt(amg_c_dprec *ph, psb_c_dvector *bh, psb_c_dvector *xh, psb_c_descriptor *cdh, const char *ctrans);
-  psb_c_i_t amg_c_dprecfree(amg_c_dprec *ph);
-  psb_c_i_t amg_c_dprecbld_opt(psb_c_dspmat *ah, psb_c_descriptor *cdh, 
+  psb_i_t amg_c_dprecinit(psb_c_ctxt cctxt, amg_c_dprec *ph, const char *ptype);
+  psb_i_t amg_c_dprecseti(amg_c_dprec *ph, const char *what, psb_i_t val);
+  psb_i_t amg_c_dprecsetc(amg_c_dprec *ph, const char *what, const char *val);
+  psb_i_t amg_c_dprecsetr(amg_c_dprec *ph, const char *what, double val);
+  psb_i_t amg_c_dpreccseti_idx(amg_c_dprec *ph, const char *what, psb_i_t val, psb_i_t idx);
+  psb_i_t amg_c_dpreccsetr_idx(amg_c_dprec *ph, const char *what, double val, psb_i_t idx);
+  psb_i_t amg_c_dpreccsetc_idx(amg_c_dprec *ph, const char *what, const char *val, psb_i_t idx);
+  psb_i_t amg_c_dpreccseti_pos(amg_c_dprec *ph, const char *what, psb_i_t val, const char *pos);
+  psb_i_t amg_c_dpreccsetr_pos(amg_c_dprec *ph, const char *what, double val, const char *pos);
+  psb_i_t amg_c_dpreccsetc_pos(amg_c_dprec *ph, const char *what, const char *val, const char *pos);
+  psb_i_t amg_c_dpreccseti_lev(amg_c_dprec *ph, const char *what, psb_i_t val, psb_i_t ilev, psb_i_t ilmax);
+  psb_i_t amg_c_dpreccsetr_lev(amg_c_dprec *ph, const char *what, double val, psb_i_t ilev, psb_i_t ilmax);
+  psb_i_t amg_c_dpreccsetc_lev(amg_c_dprec *ph, const char *what, const char *val, psb_i_t ilev, psb_i_t ilmax);
+  psb_i_t amg_c_dpreccseti_opt(amg_c_dprec *ph, const char *what, psb_i_t val, psb_i_t ilev, psb_i_t ilmax, const char *pos, psb_i_t idx);
+  psb_i_t amg_c_dpreccsetr_opt(amg_c_dprec *ph, const char *what, double val, psb_i_t ilev, psb_i_t ilmax, const char *pos, psb_i_t idx);
+  psb_i_t amg_c_dpreccsetc_opt(amg_c_dprec *ph, const char *what, const char *val, psb_i_t ilev, psb_i_t ilmax, const char *pos, psb_i_t idx);
+  psb_i_t amg_c_dprecbld(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph);
+  psb_i_t amg_c_dhierarchy_build(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph);
+  psb_i_t amg_c_dsmoothers_build(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph);
+  psb_i_t amg_c_dsmoothers_build_opt(psb_c_dspmat *ah, psb_c_descriptor *cdh, amg_c_dprec *ph, const char *afmt, const char *chfmt);
+  psb_i_t amg_c_dprecapply(amg_c_dprec *ph, psb_c_dvector *bh, psb_c_dvector *xh, psb_c_descriptor *cdh);
+  psb_i_t amg_c_dprecapply_opt(amg_c_dprec *ph, psb_c_dvector *bh, psb_c_dvector *xh, psb_c_descriptor *cdh, const char *ctrans);
+  psb_i_t amg_c_dprecfree(amg_c_dprec *ph);
+  psb_i_t amg_c_dprecbld_opt(psb_c_dspmat *ah, psb_c_descriptor *cdh, 
 			  amg_c_dprec *ph, const char *afmt);
-  psb_c_i_t amg_c_ddescr(amg_c_dprec *ph);
-  psb_c_i_t amg_c_dallocate_wrk(amg_c_dprec *ph, const char *chfmt);
+  psb_i_t amg_c_ddescr(amg_c_dprec *ph);
+  psb_i_t amg_c_dallocate_wrk(amg_c_dprec *ph, const char *chfmt);
 
-  psb_c_i_t amg_c_dkrylov(const char *method, psb_c_dspmat *ah, 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);
 

cbind/amgprec/amg_c_zprec.c

@@ -11,7 +11,7 @@ amg_c_zprec* amg_c_new_zprec()
 }
 
 
-psb_c_i_t amg_c_delete_zprec(amg_c_zprec* p)
+psb_i_t amg_c_delete_zprec(amg_c_zprec* p)
 {
   int iret;
   iret=amg_c_zprecfree(p);

cbind/amgprec/amg_c_zprec.h

@@ -19,26 +19,38 @@ extern "C"
   } amg_c_zprec;
 
   amg_c_zprec *amg_c_zprec_new();
-  psb_c_i_t amg_c_zprec_delete(amg_c_zprec *p);
+  psb_i_t amg_c_zprec_delete(amg_c_zprec *p);
 
-  psb_c_i_t amg_c_zprecinit(psb_c_ctxt cctxt, amg_c_zprec *ph, const char *ptype);
-  psb_c_i_t amg_c_zprecseti(amg_c_zprec *ph, const char *what, psb_c_i_t val);
-  psb_c_i_t amg_c_zprecsetc(amg_c_zprec *ph, const char *what, const char *val);
-  psb_c_i_t amg_c_zprecsetr(amg_c_zprec *ph, const char *what, double val);
-  psb_c_i_t amg_c_zprecbld(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph);
-  psb_c_i_t amg_c_zhierarchy_build(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph);
-  psb_c_i_t amg_c_zsmoothers_build(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph);
-  psb_c_i_t amg_c_zsmoothers_build_opt(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph, const char *afmt, const char *chfmt);
-  psb_c_i_t amg_c_zprecapply(amg_c_zprec *ph, psb_c_zvector *bh, psb_c_zvector *xh, psb_c_descriptor *cdh);
-  psb_c_i_t amg_c_zprecapply_opt(amg_c_zprec *ph, psb_c_zvector *bh, psb_c_zvector *xh, psb_c_descriptor *cdh, const char *ctrans);
-  psb_c_i_t amg_c_zprecfree(amg_c_zprec *ph);
-  psb_c_i_t amg_c_zprecbld_opt(psb_c_zspmat *ah, psb_c_descriptor *cdh,
+  psb_i_t amg_c_zprecinit(psb_c_ctxt cctxt, amg_c_zprec *ph, const char *ptype);
+  psb_i_t amg_c_zprecseti(amg_c_zprec *ph, const char *what, psb_i_t val);
+  psb_i_t amg_c_zprecsetc(amg_c_zprec *ph, const char *what, const char *val);
+  psb_i_t amg_c_zprecsetr(amg_c_zprec *ph, const char *what, double val);
+  psb_i_t amg_c_zpreccseti_idx(amg_c_zprec *ph, const char *what, psb_i_t val, psb_i_t idx);
+  psb_i_t amg_c_zpreccsetr_idx(amg_c_zprec *ph, const char *what, double val, psb_i_t idx);
+  psb_i_t amg_c_zpreccsetc_idx(amg_c_zprec *ph, const char *what, const char *val, psb_i_t idx);
+  psb_i_t amg_c_zpreccseti_pos(amg_c_zprec *ph, const char *what, psb_i_t val, const char *pos);
+  psb_i_t amg_c_zpreccsetr_pos(amg_c_zprec *ph, const char *what, double val, const char *pos);
+  psb_i_t amg_c_zpreccsetc_pos(amg_c_zprec *ph, const char *what, const char *val, const char *pos);  
+  psb_i_t amg_c_zpreccseti_lev(amg_c_zprec *ph, const char *what, psb_i_t val, psb_i_t ilev, psb_i_t ilmax);
+  psb_i_t amg_c_zpreccsetr_lev(amg_c_zprec *ph, const char *what, double val, psb_i_t ilev, psb_i_t ilmax);
+  psb_i_t amg_c_zpreccsetc_lev(amg_c_zprec *ph, const char *what, const char *val, psb_i_t ilev, psb_i_t ilmax);
+  psb_i_t amg_c_zpreccseti_opt(amg_c_zprec *ph, const char *what, psb_i_t val, psb_i_t ilev, psb_i_t ilmax, const char *pos, psb_i_t idx);
+  psb_i_t amg_c_zpreccsetr_opt(amg_c_zprec *ph, const char *what, double val, psb_i_t ilev, psb_i_t ilmax, const char *pos, psb_i_t idx);
+  psb_i_t amg_c_zpreccsetc_opt(amg_c_zprec *ph, const char *what, const char *val, psb_i_t ilev, psb_i_t ilmax, const char *pos, psb_i_t idx);
+  psb_i_t amg_c_zprecbld(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph);
+  psb_i_t amg_c_zhierarchy_build(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph);
+  psb_i_t amg_c_zsmoothers_build(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph);
+  psb_i_t amg_c_zsmoothers_build_opt(psb_c_zspmat *ah, psb_c_descriptor *cdh, amg_c_zprec *ph, const char *afmt, const char *chfmt);
+  psb_i_t amg_c_zprecapply(amg_c_zprec *ph, psb_c_zvector *bh, psb_c_zvector *xh, psb_c_descriptor *cdh);
+  psb_i_t amg_c_zprecapply_opt(amg_c_zprec *ph, psb_c_zvector *bh, psb_c_zvector *xh, psb_c_descriptor *cdh, const char *ctrans);
+  psb_i_t amg_c_zprecfree(amg_c_zprec *ph);
+  psb_i_t amg_c_zprecbld_opt(psb_c_zspmat *ah, psb_c_descriptor *cdh,
                              amg_c_zprec *ph, const char *afmt);
 
-  psb_c_i_t amg_c_zdescr(amg_c_zprec *ph);
-  psb_c_i_t amg_c_zallocate_wrk(amg_c_zprec *ph, const char *chfmt);
+  psb_i_t amg_c_zdescr(amg_c_zprec *ph);
+  psb_i_t amg_c_zallocate_wrk(amg_c_zprec *ph, const char *chfmt);
 
-  psb_c_i_t amg_c_zkrylov(const char *method, psb_c_zspmat *ah, 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);
 

cbind/amgprec/amg_dprec_cbind_mod.F90

@@ -134,6 +134,351 @@ contains
     return
   end function amg_c_dprecsetc
 
+function amg_c_dpreccseti_idx(ph,what,val,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    integer(psb_c_ipk_), value :: val
+    integer(psb_c_ipk_), value :: idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%cseti(fwhat,val,iret,idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccseti_idx
+
+
+  function amg_c_dpreccsetr_idx(ph,what,val,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    real(c_double), value :: val
+    integer(psb_c_ipk_), value :: idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%csetr(fwhat,val,iret,idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetr_idx
+
+
+  function amg_c_dpreccsetc_idx(ph,what,val,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*)
+    integer(psb_c_ipk_), value :: idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat,fval
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+
+    call precp%csetc(fwhat,fval,iret,idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetc_idx
+
+  function amg_c_dpreccseti_pos(ph,what,val,pos) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    integer(psb_c_ipk_), value :: val
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%cseti(fwhat,val,iret,pos=trim(fpos))
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccseti_pos
+
+
+  function amg_c_dpreccsetr_pos(ph,what,val,pos) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    real(c_double), value :: val
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetr(fwhat,val,iret,pos=trim(fpos))
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetr_pos
+
+
+  function amg_c_dpreccsetc_pos(ph,what,val,pos) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*), pos(*)
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fval, fpos
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetc(fwhat,fval,iret,pos=trim(fpos))
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetc_pos
+
+  function amg_c_dpreccseti_lev(ph,what,val,ilev,ilmax) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    integer(psb_c_ipk_), value :: val, ilev, ilmax
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%cseti(fwhat,val,iret,ilev=ilev,ilmax=ilmax)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccseti_lev
+
+
+  function amg_c_dpreccsetr_lev(ph,what,val,ilev,ilmax) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    real(c_double), value :: val
+    integer(psb_c_ipk_), value :: ilev, ilmax
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%csetr(fwhat,val,iret,ilev=ilev,ilmax=ilmax)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetr_lev
+
+
+  function amg_c_dpreccsetc_lev(ph,what,val,ilev,ilmax) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*)
+    integer(psb_c_ipk_), value :: ilev, ilmax
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fval
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+
+    call precp%csetc(fwhat,fval,iret,ilev=ilev,ilmax=ilmax)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetc_lev
+
+  function amg_c_dpreccseti_opt(ph,what,val,ilev,ilmax,pos,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    integer(psb_c_ipk_), value :: val, ilev, ilmax, idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%cseti(fwhat,val,iret,ilev=ilev,ilmax=ilmax,pos=trim(fpos),idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccseti_opt
+
+
+  function amg_c_dpreccsetr_opt(ph,what,val,ilev,ilmax,pos,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    real(c_double), value :: val
+    integer(psb_c_ipk_), value :: ilev, ilmax, idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetr(fwhat,val,iret,ilev=ilev,ilmax=ilmax,pos=trim(fpos),idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetr_opt
+
+
+  function amg_c_dpreccsetc_opt(ph,what,val,ilev,ilmax,pos,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*), pos(*)
+    integer(psb_c_ipk_), value :: ilev, ilmax, idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fval, fpos
+    type(amg_dprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetc(fwhat,fval,iret,ilev=ilev,ilmax=ilmax,pos=trim(fpos),idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_dpreccsetc_opt
+
   function  amg_c_dprecbld(ah,cdh,ph) bind(c) result(res)
     implicit none
 
@@ -172,17 +517,15 @@ 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
 
@@ -206,16 +549,11 @@ 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
@@ -224,7 +562,7 @@ contains
     type(psb_dspmat_type), pointer :: ap
     type(psb_desc_type), pointer   :: descp
     character(len=80)     :: fptype
-    integer(psb_ipk_)     :: iret, act
+    integer(psb_ipk_)     :: iret
 
     res = -1
 
@@ -248,10 +586,7 @@ 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
 
@@ -381,7 +716,7 @@ contains
     type(psb_c_object_type) :: ah,cdh,ph,bh,xh
     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,&

cbind/amgprec/amg_zprec_cbind_mod.F90

@@ -135,16 +135,361 @@ contains
     return
   end function amg_c_zprecsetc
 
+  function amg_c_zpreccseti_idx(ph,what,val,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    integer(psb_c_ipk_), value :: val
+    integer(psb_c_ipk_), value :: idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%cseti(fwhat,val,iret,idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccseti_idx
+
+
+  function amg_c_zpreccsetr_idx(ph,what,val,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    real(c_double), value :: val
+    integer(psb_c_ipk_), value :: idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%csetr(fwhat,val,iret,idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetr_idx
+
+
+  function amg_c_zpreccsetc_idx(ph,what,val,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*)
+    integer(psb_c_ipk_), value :: idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat,fval
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+
+    call precp%csetc(fwhat,fval,iret,idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetc_idx
+
+  function amg_c_zpreccseti_pos(ph,what,val,pos) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    integer(psb_c_ipk_), value :: val
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%cseti(fwhat,val,iret,pos=trim(fpos))
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccseti_pos
+
+
+  function amg_c_zpreccsetr_pos(ph,what,val,pos) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    real(c_double), value :: val
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetr(fwhat,val,iret,pos=trim(fpos))
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetr_pos
+
+
+  function amg_c_zpreccsetc_pos(ph,what,val,pos) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*), pos(*)
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fval, fpos
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetc(fwhat,fval,iret,pos=trim(fpos))
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetc_pos
+
+  function amg_c_zpreccseti_lev(ph,what,val,ilev,ilmax) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    integer(psb_c_ipk_), value :: val, ilev, ilmax
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%cseti(fwhat,val,iret,ilev=ilev,ilmax=ilmax)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccseti_lev
+
+
+  function amg_c_zpreccsetr_lev(ph,what,val,ilev,ilmax) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*)
+    real(c_double), value :: val
+    integer(psb_c_ipk_), value :: ilev, ilmax
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+
+    call precp%csetr(fwhat,val,iret,ilev=ilev,ilmax=ilmax)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetr_lev
+
+
+  function amg_c_zpreccsetc_lev(ph,what,val,ilev,ilmax) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*)
+    integer(psb_c_ipk_), value :: ilev, ilmax
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fval
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+
+    call precp%csetc(fwhat,fval,iret,ilev=ilev,ilmax=ilmax)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetc_lev
+
+  function amg_c_zpreccseti_opt(ph,what,val,ilev,ilmax,pos,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    integer(psb_c_ipk_), value :: val, ilev, ilmax, idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%cseti(fwhat,val,iret,ilev=ilev,ilmax=ilmax,pos=trim(fpos),idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccseti_opt
+
+
+  function amg_c_zpreccsetr_opt(ph,what,val,ilev,ilmax,pos,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), pos(*)
+    real(c_double), value :: val
+    integer(psb_c_ipk_), value :: ilev, ilmax, idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fpos
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetr(fwhat,val,iret,ilev=ilev,ilmax=ilmax,pos=trim(fpos),idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetr_opt
+
+
+  function amg_c_zpreccsetc_opt(ph,what,val,ilev,ilmax,pos,idx) bind(c) result(res)
+    implicit none
+
+    integer(psb_c_ipk_) :: res
+    type(psb_c_object_type) :: ph
+    character(c_char)     :: what(*), val(*), pos(*)
+    integer(psb_c_ipk_), value :: ilev, ilmax, idx
+    integer(psb_ipk_)     :: iret
+    character(len=80)     :: fwhat, fval, fpos
+    type(amg_zprec_type), pointer  :: precp
+
+    res = -1
+    if (c_associated(ph%item)) then
+      call c_f_pointer(ph%item,precp)
+    else
+      return
+    end if
+
+    call psb_stringc2f(what,fwhat)
+    call psb_stringc2f(val,fval)
+    call psb_stringc2f(pos,fpos)
+
+    call precp%csetc(fwhat,fval,iret,ilev=ilev,ilmax=ilmax,pos=trim(fpos),idx=idx)
+
+    res = AMGC_ERR_FILTER(iret)
+    AMGC_ERR_HANDLE(res)
+    return
+  end function amg_c_zpreccsetc_opt
+
   function  amg_c_zprecbld(ah,cdh,ph) bind(c) result(res)
     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
 
     res = -1
 
@@ -173,17 +518,15 @@ 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
 
@@ -208,24 +551,19 @@ 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
 
@@ -250,11 +588,6 @@ 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
 

cbind/test/pargen/amgec.c

@@ -123,16 +123,16 @@ double g(double x, double y, double z)
 
 #define NBMAX       20
 
-psb_c_i_t matgen(psb_c_ctxt cctxt, psb_c_i_t nl, psb_c_i_t idim, psb_c_l_t vl[],
+psb_i_t matgen(psb_c_ctxt cctxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
 	       psb_c_dspmat *ah,psb_c_descriptor *cdh,
 	       psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
 {
-  psb_c_i_t iam, np;
-  psb_c_l_t ix, iy, iz, el,glob_row;
-  psb_c_i_t i, k, info,ret;
+  psb_i_t iam, np;
+  psb_l_t ix, iy, iz, el,glob_row;
+  psb_i_t i, k, info,ret;
   double x, y, z, deltah, sqdeltah, deltah2;
   double val[10*NBMAX], zt[NBMAX];
-  psb_c_l_t irow[10*NBMAX], icol[10*NBMAX];
+  psb_l_t irow[10*NBMAX], icol[10*NBMAX];
 
   info = 0;
   psb_c_info(cctxt,&iam,&np);
@@ -254,15 +254,15 @@ void get_hparm(FILE *fp, char *val)
 int main(int argc, char *argv[])
 {
   psb_c_ctxt *cctxt;
-  psb_c_i_t iam, np;
+  psb_i_t iam, np;
   char methd[40], ptype[40], afmt[8];
-  psb_c_i_t nparms;
-  psb_c_i_t idim,info,istop,itmax,itrace,irst,iter,ret;
+  psb_i_t nparms;
+  psb_i_t idim,info,istop,itmax,itrace,irst,iter,ret;
   amg_c_dprec *ph;
   psb_c_dspmat *ah;
   psb_c_dvector *bh, *xh, *rh; 
-  psb_c_i_t nb,nlr, nl;
-  psb_c_l_t i,ng, *vl, k;
+  psb_i_t nb,nlr, nl;
+  psb_l_t i,ng, *vl, k;
   double t1,t2,eps,err;
   double *xv, *bv, *rv;
   double one=1.0, zero=0.0, res2;
@@ -305,16 +305,16 @@ int main(int argc, char *argv[])
   psb_c_set_index_base(0);
 
   /* Simple minded BLOCK data distribution */ 
-  ng = ((psb_c_l_t) idim)*idim*idim;  
+  ng = ((psb_l_t) idim)*idim*idim;  
   nb = (ng+np-1)/np;
   nl = nb;
   if ( (ng -iam*nb) < nl) nl = ng -iam*nb; 
   fprintf(stderr,"%d: Input data %d %ld %d %d\n",iam,idim,ng,nb, nl);
-  if ((vl=malloc(nb*sizeof(psb_c_l_t)))==NULL) {
+  if ((vl=malloc(nb*sizeof(psb_l_t)))==NULL) {
     fprintf(stderr,"On %d: malloc failure\n",iam);
     psb_c_abort(*cctxt);
   }
-  i = ((psb_c_l_t)iam) * nb;
+  i = ((psb_l_t)iam) * nb;
   for (k=0; k<nl; k++)
     vl[k] = i+k; 
 

cbind/test/pargen/amgecgpu.c

@@ -126,17 +126,17 @@ double g(double x, double y, double z)
 
 #define NBMAX 20
 
-psb_c_i_t matgen(psb_c_ctxt cctxt, psb_c_i_t nl, psb_c_i_t idim, psb_c_l_t vl[],
+psb_i_t matgen(psb_c_ctxt cctxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
                psb_c_dspmat *ah, const char *afmt,
                psb_c_descriptor *cdh, const char *cdfmt,
                psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
 {
-  psb_c_i_t iam, np;
-  psb_c_l_t ix, iy, iz, el, glob_row;
-  psb_c_i_t i, k, info, ret;
+  psb_i_t iam, np;
+  psb_l_t ix, iy, iz, el, glob_row;
+  psb_i_t i, k, info, ret;
   double x, y, z, deltah, sqdeltah, deltah2;
   double val[10 * NBMAX], zt[NBMAX];
-  psb_c_l_t irow[10 * NBMAX], icol[10 * NBMAX];
+  psb_l_t irow[10 * NBMAX], icol[10 * NBMAX];
 
   info = 0;
   psb_c_info(cctxt, &iam, &np);
@@ -360,15 +360,15 @@ void get_hparm(FILE *fp, char *val)
 int main(int argc, char *argv[])
 {
   psb_c_ctxt *cctxt;
-  psb_c_i_t iam, np;
+  psb_i_t iam, np;
   char methd[40], ptype[40], afmt[8], cdfmt[8];
-  psb_c_i_t nparms;
-  psb_c_i_t idim, info, istop, itmax, itrace, irst, iter, ret;
+  psb_i_t nparms;
+  psb_i_t idim, info, istop, itmax, itrace, irst, iter, ret;
   amg_c_dprec *ph;
   psb_c_dspmat *ah;
   psb_c_dvector *bh, *xh, *rh;
-  psb_c_i_t nb, nlr, nl;
-  psb_c_l_t i, ng, *vl, k;
+  psb_i_t nb, nlr, nl;
+  psb_l_t i, ng, *vl, k;
   double t1, t2, eps, err;
   double *xv, *bv, *rv;
   double one = 1.0, zero = 0.0, res2;
@@ -439,18 +439,18 @@ int main(int argc, char *argv[])
   psb_c_set_index_base(0);
 
   /* Simple minded BLOCK data distribution */
-  ng = ((psb_c_l_t)idim) * idim * idim;
+  ng = ((psb_l_t)idim) * idim * idim;
   nb = (ng + np - 1) / np;
   nl = nb;
   if ((ng - iam * nb) < nl)
     nl = ng - iam * nb;
   fprintf(stderr, "%d: Input data %d %ld %d %d\n", iam, idim, ng, nb, nl);
-  if ((vl = malloc(nb * sizeof(psb_c_l_t))) == NULL)
+  if ((vl = malloc(nb * sizeof(psb_l_t))) == NULL)
   {
     fprintf(stderr, "On %d: malloc failure\n", iam);
     psb_c_abort(*cctxt);
   }
-  i = ((psb_c_l_t)iam) * nb;
+  i = ((psb_l_t)iam) * nb;
   for (k = 0; k < nl; k++)
     vl[k] = i + k;
 

configure.ac

@@ -601,7 +601,7 @@ fi
 ###############################################################################
 #	Parachute rules for ar and ranlib ... (could cause problems)
 ###############################################################################
-AC_PROG_AR
+AC_CHECK_TOOL([AR],[ar],[ar])
 AR="${AR} -cr"
 AC_PROG_RANLIB
 
@@ -749,9 +749,13 @@ AC_LANG([C])
 PAC_CHECK_MUMPS
 #
 #    1. Enable even with LPK=8, internally it will check if
-#    	the problem size fits into 4 bytes, very likely since we
-#    	are mostly using MUMPS at coarse level.
+#       the problem size fits into 4 bytes, very likely since we
+#       are mostly using MUMPS at coarse level.
 #
+amg4psblas_cv_mumps_version="unknown"
+amg4psblas_cv_mumps_version_num=0
+CHAVEMUMPSVERSION=""
+CHAVEMUMPSVERSIONSTRING=""
 dnl if test "x$amg4psblas_cv_have_mumps" == "xyes" ; then
 dnl    if test "x$pac_cv_psblas_ipk" == "x8" ; then
 dnl       AC_MSG_NOTICE([PSBLAS defines PSB_IPK_ as $pac_cv_psblas_ipk. MUMPS interfacing disabled. ])
@@ -761,24 +765,64 @@ dnl       amg4psblas_cv_have_mumps=no;
 dnl    fi
 dnl fi   
 if test "x$amg4psblas_cv_have_mumps" == "xyes" ; then
+   amg_mumps_incdir=`echo "$MUMPS_INCLUDES" | sed -e 's/^-I//'`
+   amg_mumps_header=""
+   for amg_mumps_candidate in \
+      "$amg_mumps_incdir/dmumps_c.h" \
+      "$amg4psblas_cv_mumpsincdir/dmumps_c.h" \
+      "$amg4psblas_cv_mumpsdir/dmumps_c.h" \
+      "$amg4psblas_cv_mumpsdir/include/dmumps_c.h" \
+      "$amg4psblas_cv_mumpsdir/Include/dmumps_c.h"
+   do
+      if test -f "$amg_mumps_candidate" ; then
+         amg_mumps_header="$amg_mumps_candidate"
+         break
+      fi
+   done
+   if test "x$amg_mumps_header" != "x" ; then
+      amg4psblas_cv_mumps_version=`awk '/^#[ \t]*define[ \t]+MUMPS_VERSION[ \t]+/ {v=$0; sub(/^#[ \t]*define[ \t]+MUMPS_VERSION[ \t]+/,"",v); gsub(/"/,"",v); print v; exit}' "$amg_mumps_header" | tr -d '\r'`
+      if test "x$amg4psblas_cv_mumps_version" = "x" ; then
+         amg4psblas_cv_mumps_version=`awk '/part of MUMPS/ {for(i=1;i<=NF;i++) if($i=="MUMPS") {v=$(i+1); sub(/,.*/,"",v); print v; exit}}' "$amg_mumps_header" | tr -d '\r'`
+      fi
+      if test "x$amg4psblas_cv_mumps_version" = "x" ; then
+         amg4psblas_cv_mumps_version="unknown"
+      fi
+
+      # Strip leading "snapshot" (case-insensitive) and any delimiter
+      case "$amg4psblas_cv_mumps_version" in
+         [[sS]][[nN]][[aA]][[pP]][[sS]][[hH]][[oO]][[tT]]*)
+            amg4psblas_cv_mumps_version=`echo "$amg4psblas_cv_mumps_version" | sed -e 's/^[[sS]][[nN]][[aA]][[pP]][[sS]][[hH]][[oO]][[tT]][[-_ ]]*//'`
+            ;;
+      esac
+
+      amg4psblas_cv_mumps_version_num=`echo "$amg4psblas_cv_mumps_version" | sed -e 's/[^0-9]//g' -e 's/^0*//'`
+      if test "x$amg4psblas_cv_mumps_version_num" = "x" ; then
+         amg4psblas_cv_mumps_version_num=0
+      fi
+      AC_MSG_NOTICE([Configuring with MUMPS version $amg4psblas_cv_mumps_version (numeric flag$amg4psblas_cv_mumps_version_num)])
+   else
+      AC_MSG_NOTICE([Could not locate dmumps_c.h to extract the MUMPS version.])
+   fi
+   CHAVEMUMPSVERSION="#define  AMG_MUMPS_VERSION $amg4psblas_cv_mumps_version_num"
+   CHAVEMUMPSVERSIONSTRING="#define  AMG_MUMPS_VERSION_STRING \"$amg4psblas_cv_mumps_version\""
    if test "x$pac_cv_psblas_lpk" == "x8" ; then
       AC_MSG_NOTICE([PSBLAS defines PSB_LPK_ as $pac_cv_psblas_lpk. MUMPS interfacing will fail when called in global mode on very large matrices. ])
    fi
    MUMPS_LIBS="-lsmumps -ldmumps -lcmumps -lzmumps -lmumps_common -lpord"
    if test "x$amg4psblas_cv_mumpslibdir" != "x" ; then
       MUMPS_LIBS="${MUMPS_LIBS} -L$amg4psblas_cv_mumpslibdir"
-   fi 		   
+   fi            
    if test "x$pac_mumps_fmods_ok" == "xyes" ; then 
-      FDEFINES="$amg_cv_define_prepend-DAMG_HAVE_MUMPS $amg_cv_define_prepend-DAMG_HAVE_MUMPS_MODULES $MUMPS_MODULES $FDEFINES"
-       MUMPS_FLAGS="-DAMG_HAVE_MUMPS $MUMPS_MODULES"
+        FDEFINES="$amg_cv_define_prepend-DAMG_HAVE_MUMPS$amg_cv_define_prepend -DAMG_HAVE_MUMPS_MODULES $amg_cv_define_prepend-DAMG_MUMPS_VERSION=$amg4psblas_cv_mumps_version_num $MUMPS_MODULES $FDEFINES"
+        MUMPS_FLAGS="-DAMG_HAVE_MUMPS -DAMG_MUMPS_VERSION=$amg4psblas_cv_mumps_version_num $MUMPS_MODULES"
        CHAVEMUMPS="#define  AMG_HAVE_MUMPS"
        CHAVEMUMPSMODULES="#define  AMG_HAVE_MUMPS_MODULES"
        
    elif test "x$pac_mumps_fincs_ok" == "xyes" ; then
-      FDEFINES="$amg_cv_define_prepend-DAMG_HAVE_MUMPS $amg_cv_define_prepend-DAMG_HAVE_MUMPS_INCLUDES $MUMPS_FINCLUDES $FDEFINES"
-       MUMPS_FLAGS="-DAMG_HAVE_MUMPS $MUMPS_FINCLUDES"
+        FDEFINES="$amg_cv_define_prepend-DAMG_HAVE_MUMPS$amg_cv_define_prepend -DAMG_HAVE_MUMPS_INCLUDES $amg_cv_define_prepend-DAMG_MUMPS_VERSION=$amg4psblas_cv_mumps_version_num $MUMPS_FINCLUDES $FDEFINES"
+        MUMPS_FLAGS="-DAMG_HAVE_MUMPS -DAMG_MUMPS_VERSION=$amg4psblas_cv_mumps_version_num $MUMPS_FINCLUDES "
        CHAVEMUMPS="#define  AMG_HAVE_MUMPS"
-       CHAVEMUMPSINCLUDES="#define  AMG_HAVE_MUMPS_INCLUDES"	
+       CHAVEMUMPSINCLUDES="#define  AMG_HAVE_MUMPS_INCLUDES"    
    else
       # This should not happen
       MUMPS_FLAGS=""
@@ -869,6 +913,8 @@ AC_SUBST(CSLUDISTVERSION)
 AC_SUBST(CHAVEMUMPS)
 AC_SUBST(CHAVEMUMPSMODULES)
 AC_SUBST(CHAVEMUMPSINCLUDES)
+AC_SUBST(CHAVEMUMPSVERSION)
+AC_SUBST(CHAVEMUMPSVERSIONSTRING)
 AC_SUBST(CXXMATCHBOXBIT)
 
 

samples/advanced/pdegen/Makefile

@@ -3,7 +3,8 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec 
+#AMG_LIBS = -L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec $(MUMPSLIBS) $(EXTRALIBS)
+AMG_LIBS=-L$(AMGLIBDIR) -lamg_prec
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
 
 LINKOPT=
@@ -11,7 +12,8 @@ EXEDIR=./runs
 DGEN2D=amg_d_pde2d_poisson_mod.o amg_d_pde2d_exp_mod.o \
   amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
 DGEN3D=amg_d_pde3d_poisson_mod.o amg_d_pde3d_exp_mod.o \
-  amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o 	
+  amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o \
+  amg_d_pde3d_aniso_mod.o	
 SGEN2D=amg_s_pde2d_poisson_mod.o amg_s_pde2d_exp_mod.o \
   amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
 SGEN3D=amg_s_pde3d_poisson_mod.o amg_s_pde3d_exp_mod.o \

samples/advanced/pdegen/amg_d_genpde_mod.F90

@@ -53,6 +53,10 @@ module amg_d_genpde_mod
     module procedure  amg_d_gen_pde3d
   end interface amg_gen_pde3d
 
+  interface amg_gen_pde3d_aniso
+    module procedure  amg_d_gen_aniso_poisson3d
+  end interface
+
   interface
     function d_func_2d(x,y) result(val)
       import :: psb_dpk_
@@ -1006,4 +1010,645 @@ contains
     end if
     return
   end subroutine amg_d_gen_pde2d
+
+    subroutine amg_d_gen_aniso_poisson3d(ctxt,idim,a,bv,xv,desc_a,afmt,&
+       & k11,k22,k33,k12,k13,k23,g,info,f,amold,vmold,partition,nrl,iv)
+    use psb_base_mod
+    use psb_util_mod
+#if defined(PSB_OPENMP)
+    use omp_lib
+#endif
+    !
+    !   Discretizes the anisotropic Poisson equation
+    !
+    !       - div(K grad(u)) = f
+    !
+    !   with Dirichlet boundary conditions u = g on the unit cube.
+    !
+    !   The symmetric diffusion tensor is
+    !
+    !       K = [ k11 k12 k13 ]
+    !           [ k12 k22 k23 ]
+    !           [ k13 k23 k33 ].
+    !
+    !   The coefficient functions are evaluated at the grid point.  The
+    !   19-point stencil below is the centered discretization of the
+    !   constant-tensor operator
+    !
+    !       -k11 u_xx - k22 u_yy - k33 u_zz
+    !       -2 k12 u_xy - 2 k13 u_xz - 2 k23 u_yz.
+    !
+    !   For variable coefficient functions this is therefore a
+    !   non-divergence-form discretization; for the intended anisotropic
+    !   Poisson benchmarks the tensor is normally constant.
+    !
+    implicit none
+    procedure(d_func_3d)  :: k11,k22,k33,k12,k13,k23,g
+    integer(psb_ipk_)     :: idim
+    type(psb_dspmat_type) :: a
+    type(psb_d_vect_type) :: xv,bv
+    type(psb_desc_type)   :: desc_a
+    integer(psb_ipk_)     :: info
+    type(psb_ctxt_type)   :: ctxt
+    character             :: afmt*5
+    procedure(d_func_3d), optional :: f
+    class(psb_d_base_sparse_mat), optional :: amold
+    class(psb_d_base_vect_type), optional :: vmold
+    integer(psb_ipk_), optional :: partition, nrl,iv(:)
+
+    ! Local variables.
+
+    integer(psb_ipk_), parameter :: nb=20
+    type(psb_d_csc_sparse_mat)  :: acsc
+    type(psb_d_coo_sparse_mat)  :: acoo
+    type(psb_d_csr_sparse_mat)  :: acsr
+    integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
+    integer(psb_lpk_) :: m,n,glob_row,nt
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    ! For 3D partition
+    ! Note: integer control variables going directly into an MPI call
+    ! must be 4 bytes, i.e. psb_mpk_
+    integer(psb_mpk_) :: npdims(3), npp, minfo
+    integer(psb_ipk_) :: npx,npy,npz, iamx,iamy,iamz,mynx,myny,mynz
+    integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:)
+    ! Process grid
+    integer(psb_ipk_) :: np, iam
+    integer(psb_ipk_) :: icoeff
+    integer(psb_lpk_), allocatable     :: myidx(:)
+    ! deltah dimension of each grid cell
+    ! deltat discretization time
+    real(psb_dpk_)            :: deltah, sqdeltah, deltah2
+    real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero
+    real(psb_dpk_)    :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
+    integer(psb_ipk_) :: err_act
+    procedure(d_func_3d), pointer :: f_
+    character(len=20)  :: name, ch_err,tmpfmt
+
+    info = psb_success_
+    name = 'd_create_aniso_poisson3d'
+    call psb_erractionsave(err_act)
+
+    call psb_info(ctxt, iam, np)
+
+
+    if (present(f)) then
+      f_ => f
+    else
+      f_ => d_null_func_3d
+    end if
+
+    if (present(partition)) then
+      if ((1<= partition).and.(partition <= 3)) then
+        partition_ = partition
+      else
+        write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
+        partition_ = 3
+      end if
+    else
+      partition_ = 3
+    end if
+    deltah   = done/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.0_psb_dpk_* deltah
+
+    if (present(partition)) then
+      if ((1<= partition).and.(partition <= 3)) then
+        partition_ = partition
+      else
+        write(*,*) 'Invalid partition choice ',partition,' defaulting to 3'
+        partition_ = 3
+      end if
+    else
+      partition_ = 3
+    end if
+
+    ! initialize array descriptor and sparse matrix storage. provide an
+    ! estimate of the number of non zeroes
+
+    m   = (1_psb_lpk_*idim)*idim*idim
+    n   = m
+    ! A full symmetric tensor diffusion operator uses a 19-point stencil.
+    nnz = 19*((n+np-1)/np)
+    if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
+    t0 = psb_wtime()
+    select case(partition_)
+    case(1)
+      ! A BLOCK partition
+      if (present(nrl)) then
+        nr = nrl
+      else
+        !
+        ! Using a simple BLOCK distribution.
+        !
+        nt = (m+np-1)/np
+        nr = max(0,min(nt,m-(iam*nt)))
+      end if
+
+      nt = nr
+      call psb_sum(ctxt,nt)
+      if (nt /= m) then
+        write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m
+        info = -1
+        call psb_barrier(ctxt)
+        call psb_abort(ctxt)
+        return
+      end if
+
+      !
+      ! First example  of use of CDALL: specify for each process a number of
+      ! contiguous rows
+      !
+      call psb_cdall(ctxt,desc_a,info,nl=nr)
+      if (info /=0) goto 9999
+      myidx = desc_a%get_global_indices()
+      nlr = size(myidx)
+
+    case(2)
+      ! A  partition  defined by the user through IV
+
+      if (present(iv)) then
+        if (size(iv) /= m) then
+          write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m
+          info = -1
+          call psb_barrier(ctxt)
+          call psb_abort(ctxt)
+          return
+        end if
+      else
+        write(psb_err_unit,*) iam, 'Initialization error: IV not present'
+        info = -1
+        call psb_barrier(ctxt)
+        call psb_abort(ctxt)
+        return
+      end if
+
+      !
+      ! Second example  of use of CDALL: specify for each row the
+      ! process that owns it
+      !
+      call psb_cdall(ctxt,desc_a,info,vg=iv)
+      if (info /=0) goto 9999
+      myidx = desc_a%get_global_indices()
+      nlr = size(myidx)
+
+    case(3)
+      ! A 3-dimensional partition
+
+      ! A nifty MPI function will split the process list
+      npdims = 0
+#if defined(PSB_SERIAL_MPI)
+      npdims = 1
+#else
+      npp = np
+      call mpi_dims_create(npp,3,npdims,minfo)
+#endif
+      npx = npdims(1)
+      npy = npdims(2)
+      npz = npdims(3)
+
+      allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz))
+      ! We can reuse idx2ijk for process indices as well.
+      call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=mzero)
+      ! Now let's split the 3D cube in hexahedra
+      call dist1Didx(bndx,idim,npx)
+      mynx = bndx(iamx+1)-bndx(iamx)
+      call dist1Didx(bndy,idim,npy)
+      myny = bndy(iamy+1)-bndy(iamy)
+      call dist1Didx(bndz,idim,npz)
+      mynz = bndz(iamz+1)-bndz(iamz)
+
+      ! How many indices do I own?
+      nlr = mynx*myny*mynz
+      allocate(myidx(nlr))
+      ! Now, let's generate the list of indices I own
+      nr = 0
+      do i=bndx(iamx),bndx(iamx+1)-1
+        do j=bndy(iamy),bndy(iamy+1)-1
+          do k=bndz(iamz),bndz(iamz+1)-1
+            nr = nr + 1
+            call ijk2idx(myidx(nr),i,j,k,idim,idim,idim)
+          end do
+        end do
+      end do
+      if (nr /= nlr) then
+        write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',&
+             & nr,nlr,mynx,myny,mynz
+        info = -1
+        call psb_barrier(ctxt)
+        call psb_abort(ctxt)
+      end if
+
+      !
+      ! Third example  of use of CDALL: specify for each process
+      ! the set of global indices it owns.
+      !
+      call psb_cdall(ctxt,desc_a,info,vl=myidx)
+      if (info /=0) goto 9999
+
+      !
+      ! Specify process topology
+      !
+      block
+        !
+        ! Use adjcncy methods 
+        ! 
+        integer(psb_ipk_), allocatable :: neighbours(:)
+        integer(psb_mpk_) :: cnt
+        logical, parameter :: debug_adj=.true.
+        if (debug_adj.and.(np > 1)) then 
+          cnt = 0
+          allocate(neighbours(np))
+          if (iamx < npx-1) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx+1,iamy,iamz,npx,npy,npz,base=mzero)
+          end if
+          if (iamy < npy-1) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy+1,iamz,npx,npy,npz,base=mzero)
+          end if
+          if (iamz < npz-1) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy,iamz+1,npx,npy,npz,base=mzero)
+          end if
+          if (iamx >0) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx-1,iamy,iamz,npx,npy,npz,base=mzero)
+          end if
+          if (iamy >0) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy-1,iamz,npx,npy,npz,base=mzero)
+          end if
+          if (iamz >0) then
+            cnt = cnt + 1 
+            call ijk2idx(neighbours(cnt),iamx,iamy,iamz-1,npx,npy,npz,base=mzero)
+          end if
+          call psb_realloc(cnt, neighbours,info)
+          call desc_a%set_p_adjcncy(neighbours)
+          !write(0,*) iam,' Check on neighbours: ',desc_a%get_p_adjcncy()
+        end if
+      end block
+      
+    case default
+      write(psb_err_unit,*) iam, 'Initialization error: should not get here'
+      info = -1
+      call psb_barrier(ctxt)
+      call psb_abort(ctxt)
+      return
+    end select
+
+
+    if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz)
+    ! define  rhs from boundary conditions; also build initial guess
+    if (info == psb_success_) call psb_geall(xv,desc_a,info)
+    if (info == psb_success_) call psb_geall(bv,desc_a,info)
+
+    call psb_barrier(ctxt)
+    talc = psb_wtime()-t0
+
+    call psb_barrier(ctxt)
+    t1 = psb_wtime()
+    ! Disable OMP here for the time being
+    !
+    ! For a symmetric diffusion tensor the centered discretization has a
+    ! 19-point stencil: 1 center, 6 face neighbours and 12 edge neighbours.
+    block
+      integer(psb_ipk_) :: i,j,k,ii,ib,icoeff, ix,iy,iz, ith,nth
+      integer(psb_lpk_) :: glob_row
+      integer(psb_lpk_), allocatable :: irow(:),icol(:)
+      real(psb_dpk_), allocatable :: val(:)
+      real(psb_dpk_) :: x,y,z,zt(nb)
+      real(psb_dpk_) :: v
+      real(psb_dpk_) :: c11,c22,c33,c12,c13,c23
+#if defined(PSB_OPENMP)
+      nth = omp_get_num_threads()
+      ith = omp_get_thread_num()
+#else
+      nth = 1
+      ith = 0
+#endif
+      allocate(val(20*nb),irow(20*nb),icol(20*nb),stat=info)
+      if (info /= psb_success_) then
+        info=psb_err_alloc_dealloc_
+        call psb_errpush(info,name)
+      endif
+
+      do ii=1,nlr,nb
+        if (info /= psb_success_) cycle
+        ib = min(nb,nlr-ii+1)
+        icoeff = 1
+
+        do k=1,ib
+          i=ii+k-1
+          glob_row=myidx(i)
+
+          call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim)
+          x = (ix-1)*deltah
+          y = (iy-1)*deltah
+          z = (iz-1)*deltah
+
+          zt(k) = f_(x,y,z)
+
+          c11 = k11(x,y,z)
+          c22 = k22(x,y,z)
+          c33 = k33(x,y,z)
+          c12 = k12(x,y,z)
+          c13 = k13(x,y,z)
+          c23 = k23(x,y,z)
+
+          !
+          ! x-direction neighbours
+          !
+          v = -c11/sqdeltah
+          if (ix == 1) then
+            zt(k) = g(dzero,y,z)*(-v) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          endif
+
+          v = -c11/sqdeltah
+          if (ix == idim) then
+            zt(k) = g(done,y,z)*(-v) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          endif
+
+          !
+          ! y-direction neighbours
+          !
+          v = -c22/sqdeltah
+          if (iy == 1) then
+            zt(k) = g(x,dzero,z)*(-v) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          endif
+
+          v = -c22/sqdeltah
+          if (iy == idim) then
+            zt(k) = g(x,done,z)*(-v) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          endif
+
+          !
+          ! z-direction neighbours
+          !
+          v = -c33/sqdeltah
+          if (iz == 1) then
+            zt(k) = g(x,y,dzero)*(-v) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          endif
+
+          v = -c33/sqdeltah
+          if (iz == idim) then
+            zt(k) = g(x,y,done)*(-v) + zt(k)
+          else
+            call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          endif
+
+          !
+          ! xy mixed derivative:
+          !
+          !   -2 c12 u_xy
+          !
+          ! gives
+          !
+          !   -c12/(2h^2) [u(i+1,j+1)-u(i+1,j-1)
+          !                 -u(i-1,j+1)+u(i-1,j-1)].
+          !
+          v = -c12/(2.0_psb_dpk_*sqdeltah)
+
+          ! (ix+1,iy+1,iz)
+          if ((ix < idim).and.(iy < idim)) then
+            call ijk2idx(icol(icoeff),ix+1,iy+1,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(done,x,ix==idim),merge(done,y,iy==idim),z)*(-v)+zt(k)
+          endif
+
+          ! (ix+1,iy-1,iz)
+          if ((ix < idim).and.(iy > 1)) then
+            call ijk2idx(icol(icoeff),ix+1,iy-1,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=-v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(done,x,ix==idim),merge(dzero,y,iy==1),z)*v+zt(k)
+          endif
+
+          ! (ix-1,iy+1,iz)
+          if ((ix > 1).and.(iy < idim)) then
+            call ijk2idx(icol(icoeff),ix-1,iy+1,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=-v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(dzero,x,ix==1),merge(done,y,iy==idim),z)*v+zt(k)
+          endif
+
+          ! (ix-1,iy-1,iz)
+          if ((ix > 1).and.(iy > 1)) then
+            call ijk2idx(icol(icoeff),ix-1,iy-1,iz,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(dzero,x,ix==1),merge(dzero,y,iy==1),z)*(-v)+zt(k)
+          endif
+
+          !
+          ! xz mixed derivative
+          !
+          v = -c13/(2.0_psb_dpk_*sqdeltah)
+
+          ! (ix+1,iy,iz+1)
+          if ((ix < idim).and.(iz < idim)) then
+            call ijk2idx(icol(icoeff),ix+1,iy,iz+1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(done,x,ix==idim),y,merge(done,z,iz==idim))*(-v)+zt(k)
+          endif
+
+          ! (ix+1,iy,iz-1)
+          if ((ix < idim).and.(iz > 1)) then
+            call ijk2idx(icol(icoeff),ix+1,iy,iz-1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=-v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(done,x,ix==idim),y,merge(dzero,z,iz==1))*v+zt(k)
+          endif
+
+          ! (ix-1,iy,iz+1)
+          if ((ix > 1).and.(iz < idim)) then
+            call ijk2idx(icol(icoeff),ix-1,iy,iz+1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=-v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(dzero,x,ix==1),y,merge(done,z,iz==idim))*v+zt(k)
+          endif
+
+          ! (ix-1,iy,iz-1)
+          if ((ix > 1).and.(iz > 1)) then
+            call ijk2idx(icol(icoeff),ix-1,iy,iz-1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(merge(dzero,x,ix==1),y,merge(dzero,z,iz==1))*(-v)+zt(k)
+          endif
+
+          !
+          ! yz mixed derivative
+          !
+          v = -c23/(2.0_psb_dpk_*sqdeltah)
+
+          ! (ix,iy+1,iz+1)
+          if ((iy < idim).and.(iz < idim)) then
+            call ijk2idx(icol(icoeff),ix,iy+1,iz+1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(x,merge(done,y,iy==idim),merge(done,z,iz==idim))*(-v)+zt(k)
+          endif
+
+          ! (ix,iy+1,iz-1)
+          if ((iy < idim).and.(iz > 1)) then
+            call ijk2idx(icol(icoeff),ix,iy+1,iz-1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=-v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(x,merge(done,y,iy==idim),merge(dzero,z,iz==1))*v+zt(k)
+          endif
+
+          ! (ix,iy-1,iz+1)
+          if ((iy > 1).and.(iz < idim)) then
+            call ijk2idx(icol(icoeff),ix,iy-1,iz+1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=-v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(x,merge(dzero,y,iy==1),merge(done,z,iz==idim))*v+zt(k)
+          endif
+
+          ! (ix,iy-1,iz-1)
+          if ((iy > 1).and.(iz > 1)) then
+            call ijk2idx(icol(icoeff),ix,iy-1,iz-1,idim,idim,idim)
+            irow(icoeff)=glob_row
+            val(icoeff)=v
+            icoeff=icoeff+1
+          else
+            zt(k)=g(x,merge(dzero,y,iy==1),merge(dzero,z,iz==1))*(-v)+zt(k)
+          endif
+
+          !
+          ! Center coefficient.
+          !
+          val(icoeff)=2.0_psb_dpk_*(c11+c22+c33)/sqdeltah
+          call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim)
+          irow(icoeff)=glob_row
+          icoeff=icoeff+1
+
+        end do
+
+        call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info)
+        if(info /= psb_success_) cycle
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info)
+        if(info /= psb_success_) cycle
+        zt(:)=dzero
+        call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
+        if(info /= psb_success_) cycle
+      end do
+
+      deallocate(val,irow,icol)
+    end block
+    tgen = psb_wtime()-t1
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='insert rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+
+
+    call psb_barrier(ctxt)
+    t1 = psb_wtime()
+    call psb_cdasb(desc_a,info)
+    tcdasb = psb_wtime()-t1
+    call psb_barrier(ctxt)
+    t1 = psb_wtime()
+    if (info == psb_success_) then
+      if (present(amold)) then
+        call psb_spasb(a,desc_a,info,mold=amold)
+      else
+        call psb_spasb(a,desc_a,info,afmt=afmt)
+      end if
+    end if
+    call psb_barrier(ctxt)
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='asb rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+    if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold)
+    if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold)
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='asb rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+    tasb = psb_wtime()-t1
+    call psb_barrier(ctxt)
+    ttot = psb_wtime() - t0
+
+    call psb_amx(ctxt,talc)
+    call psb_amx(ctxt,tgen)
+    call psb_amx(ctxt,tasb)
+    call psb_amx(ctxt,ttot)
+    if(iam == psb_root_) then
+      tmpfmt = a%get_fmt()
+      write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')&
+           &   tmpfmt
+      write(psb_out_unit,'("-allocation  time : ",es12.5)') talc
+      write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen
+      write(psb_out_unit,'("-desc asbly  time : ",es12.5)') tcdasb
+      write(psb_out_unit,'("- mat asbly  time : ",es12.5)') tasb
+      write(psb_out_unit,'("-total       time : ",es12.5)') ttot
+
+    end if
+    call psb_erractionrestore(err_act)
+    return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+    return
+  end subroutine amg_d_gen_aniso_poisson3d
+
 end module amg_d_genpde_mod

samples/advanced/pdegen/amg_d_pde3d.F90

@@ -74,6 +74,7 @@ program amg_d_pde3d
   use amg_d_pde3d_exp_mod
   use amg_d_pde3d_box_mod
   use amg_d_pde3d_gauss_mod
+  use amg_d_pde3d_aniso_mod
   use amg_d_genpde_mod
 #if defined(PSB_OPENMP)
   use omp_lib
@@ -90,6 +91,9 @@ program amg_d_pde3d
   ! miscellaneous
   real(psb_dpk_) :: t1, t2, tprec, thier, tslv, tsmth, tpgen
 
+  ! anisotropy parameters
+  real(psb_dpk_) :: epsilonaniso, thetaaniso, anisovec(2)
+
   ! sparse matrix and preconditioner
   type(psb_dspmat_type) :: a
   type(amg_dprec_type)  :: prec
@@ -160,6 +164,12 @@ program amg_d_pde3d
     integer(psb_ipk_)  :: fill         ! fill-in for incomplete LU factorization
     integer(psb_ipk_)  :: invfill      ! Inverse fill-in for INVK
     real(psb_dpk_)      :: thr          ! threshold for ILUT factorization
+    integer(psb_ipk_)  :: mumps_blr_icntl35 ! BLR activation/format (ICNTL(35))
+    integer(psb_ipk_)  :: mumps_blr_icntl36 ! BLR variant (ICNTL(36))
+    integer(psb_ipk_)  :: mumps_blr_icntl37 ! BLR block size (ICNTL(37))
+    integer(psb_ipk_)  :: mumps_blr_icntl38 ! BLR compression option (ICNTL(38))
+    real(psb_dpk_)     :: mumps_blr_cntl7   ! BLR drop parameter (CNTL(7))
+    integer(psb_ipk_)  :: mumps_sym_value ! Symmetry value for MUMPS 
 
     ! AMG post-smoother; ignored by 1-lev preconditioner
     character(len=32)  :: smther2      ! post-smoother type: BJAC, AS
@@ -253,7 +263,9 @@ program amg_d_pde3d
   !
   !  get parameters
   !
-  call get_parms(ctxt,afmt,idim,s_choice,p_choice,pdecoeff)
+  call get_parms(ctxt,afmt,idim,s_choice,p_choice,pdecoeff,epsilonaniso,thetaaniso)
+  anisovec(1) = epsilonaniso
+  anisovec(2) = thetaaniso
 
   !
   !  allocate and fill in the coefficient matrix, rhs and initial guess
@@ -278,6 +290,13 @@ program amg_d_pde3d
     call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
          & a1_gauss,a2_gauss,a3_gauss,&
          & b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
+  case("ANISO")
+    call pde_set_parm3d_aniso(anisovec)
+    if (iam == psb_root_) then
+      write(psb_out_unit,'("Anisotropy \epsilon = ",es12.5," \theta = ",es12.5)') anisovec(1), anisovec(2)
+    end if
+    call amg_d_gen_aniso_poisson3d(ctxt,idim,a,b,x,desc_a,afmt, &
+         & k11,k22,k33,k12,k13,k23,g_aniso,info)
   case default
     info=psb_err_from_subroutine_
     ch_err='amg_gen_pdecoeff'
@@ -319,8 +338,14 @@ program amg_d_pde3d
     call prec%set('solver_sweeps',   p_choice%ssweeps,   info)
     call prec%set('poly_degree',     p_choice%degree,    info)
     call prec%set('poly_variant',    p_choice%pvariant,  info)
-    if (psb_toupper(p_choice%solve)=='MUMPS') &
-         & call prec%set('mumps_loc_glob','local_solver',info)
+    if (psb_toupper(p_choice%solve)=='MUMPS') then
+      call prec%set('mumps_loc_glob','local_solver',info)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_)
+      call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_)
+    end if
     call prec%set('sub_fillin',      p_choice%fill,    info)
     call prec%set('sub_iluthrs',     p_choice%thr,     info)
 
@@ -331,8 +356,15 @@ program amg_d_pde3d
     call prec%set('sub_prol',        p_choice%prol,    info)
     call prec%set('sub_solve',       p_choice%solve,   info)
     call prec%set('solver_sweeps',   p_choice%ssweeps,   info)
-    if (psb_toupper(p_choice%solve)=='MUMPS') &
-         & call prec%set('mumps_loc_glob','local_solver',info)
+    if (psb_toupper(p_choice%solve)=='MUMPS') then
+      call prec%set('mumps_loc_glob','local_solver',info)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_)
+      call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_)
+      call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_)
+      call prec%set('MUMPS_SYM',p_choice%mumps_sym_value,info)
+    end if
     call prec%set('sub_fillin',      p_choice%fill,    info)
     call prec%set('sub_iluthrs',     p_choice%thr,     info)
 
@@ -391,8 +423,20 @@ program amg_d_pde3d
         call prec%set('ainv_alg', p_choice%variant,   info)
       case default
         call prec%set('sub_solve',       p_choice%solve,   info)
-        if (psb_toupper(p_choice%solve)=='MUMPS') &
-             & call prec%set('mumps_loc_glob','local_solver',info)
+        if (psb_toupper(p_choice%solve)=='MUMPS') then
+          if (psb_toupper(p_choice%smther) == 'RICHARDS' .or. &
+               & psb_toupper(p_choice%smther) == 'RICHARDSON') then
+            call prec%set('mumps_loc_glob','global_solver',info)
+          else
+            call prec%set('mumps_loc_glob','local_solver',info)
+          end if
+          call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_)
+          call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_)
+          call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_)
+          call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_)
+          call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_)
+          call prec%set('MUMPS_SYM',p_choice%mumps_sym_value,info)
+        end if
       end select
       call prec%set('solver_sweeps',   p_choice%ssweeps,   info)
       call prec%set('sub_fillin',      p_choice%fill,       info)
@@ -428,8 +472,20 @@ program amg_d_pde3d
           call prec%set('ainv_alg', p_choice%variant2,   info)
         case default
           call prec%set('sub_solve',       p_choice%solve2,   info, pos='post')
-          if (psb_toupper(p_choice%solve2)=='MUMPS') &
-               & call prec%set('mumps_loc_glob','local_solver',info)
+          if (psb_toupper(p_choice%solve2)=='MUMPS') then
+            if (psb_toupper(p_choice%smther2) == 'RICHARDS' .or. &
+                 & psb_toupper(p_choice%smther2) == 'RICHARDSON') then
+              call prec%set('mumps_loc_glob','global_solver',info)
+            else
+              call prec%set('mumps_loc_glob','local_solver',info)
+            end if
+            call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl35,info,idx=35_psb_ipk_)
+            call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl36,info,idx=36_psb_ipk_)
+            call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl37,info,idx=37_psb_ipk_)
+            call prec%set('mumps_ipar_entry',p_choice%mumps_blr_icntl38,info,idx=38_psb_ipk_)
+            call prec%set('mumps_rpar_entry',p_choice%mumps_blr_cntl7,info,idx=7_psb_ipk_)
+            call prec%set('MUMPS_SYM',p_choice%mumps_sym_value,info)
+          end if
         end select
         call prec%set('solver_sweeps',   p_choice%ssweeps2,   info,pos='post')
         call prec%set('sub_fillin',      p_choice%fill2,     info,pos='post')
@@ -610,7 +666,7 @@ contains
   !
   ! get iteration parameters from standard input
   !
-  subroutine get_parms(ctxt,afmt,idim,solve,prec,pdecoeff)
+  subroutine get_parms(ctxt,afmt,idim,solve,prec,pdecoeff,epsilonaniso,thetaaniso)
 
     implicit none
 
@@ -620,6 +676,7 @@ contains
     type(solverdata)    :: solve
     type(precdata)      :: prec
     character(len=*)    :: pdecoeff
+    real(psb_dpk_)      :: epsilonaniso, thetaaniso
     integer(psb_ipk_)   :: iam, nm, np, inp_unit
     character(len=1024)   :: filename
 
@@ -647,6 +704,8 @@ contains
       call read_data(afmt,inp_unit)            ! matrix storage format
       call read_data(idim,inp_unit)            ! Discretization grid size
       call read_data(pdecoeff,inp_unit)        ! PDE Coefficients
+      call read_data(epsilonaniso,inp_unit)    ! Anysotropy epsilon parameter
+      call read_data(thetaaniso,inp_unit)      ! Anysotropy angle
       ! Krylov solver data
       call read_data(solve%kmethd,inp_unit)    ! Krylov solver
       call read_data(solve%istopc,inp_unit)    ! stopping criterion
@@ -674,6 +733,12 @@ contains
       call read_data(prec%fill,inp_unit)       ! fill-in for incomplete LU
       call read_data(prec%invfill,inp_unit)    !Inverse fill-in for INVK
       call read_data(prec%thr,inp_unit)        ! threshold for ILUT
+      call read_data(prec%mumps_blr_icntl35,inp_unit) ! MUMPS ICNTL(35)
+      call read_data(prec%mumps_blr_icntl36,inp_unit) ! MUMPS ICNTL(36)
+      call read_data(prec%mumps_blr_icntl37,inp_unit) ! MUMPS ICNTL(37)
+      call read_data(prec%mumps_blr_icntl38,inp_unit) ! MUMPS ICNTL(38)
+      call read_data(prec%mumps_blr_cntl7,inp_unit)   ! MUMPS CNTL(7)
+      call read_data(prec%mumps_sym_value,inp_unit)   ! MUMPS SYMMETRY
       ! Second smoother/ AMG post-smoother (if NONE ignored in main)
       call read_data(prec%smther2,inp_unit)     ! smoother type
       call read_data(prec%jsweeps2,inp_unit)    ! (post-)smoother sweeps
@@ -748,6 +813,8 @@ contains
     call psb_bcast(ctxt,afmt)
     call psb_bcast(ctxt,idim)
     call psb_bcast(ctxt,pdecoeff)
+    call psb_bcast(ctxt,epsilonaniso)
+    call psb_bcast(ctxt,thetaaniso)
 
     call psb_bcast(ctxt,solve%kmethd)
     call psb_bcast(ctxt,solve%istopc)
@@ -775,6 +842,12 @@ contains
     call psb_bcast(ctxt,prec%fill)
     call psb_bcast(ctxt,prec%invfill)
     call psb_bcast(ctxt,prec%thr)
+    call psb_bcast(ctxt,prec%mumps_blr_icntl35)
+    call psb_bcast(ctxt,prec%mumps_blr_icntl36)
+    call psb_bcast(ctxt,prec%mumps_blr_icntl37)
+    call psb_bcast(ctxt,prec%mumps_blr_icntl38)
+    call psb_bcast(ctxt,prec%mumps_blr_cntl7)
+    call psb_bcast(ctxt,prec%mumps_sym_value)
     ! broadcast second (post-)smoother
     call psb_bcast(ctxt,prec%smther2)
     call psb_bcast(ctxt,prec%jsweeps2)

samples/advanced/pdegen/amg_d_pde3d_aniso_mod.f90

@@ -0,0 +1,161 @@
+!
+!
+!                             AMG4PSBLAS version 1.0
+!    Algebraic Multigrid Package
+!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Fabio Durastante
+!
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the AMG4PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific prior written permission.
+!
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!
+module amg_d_pde3d_aniso_mod
+
+  use psb_base_mod, only : psb_dpk_, done, dzero
+
+  real(psb_dpk_), save, private :: epsilon = done/80
+  real(psb_dpk_), save, private :: theta   = dzero
+
+contains
+
+  subroutine pde_set_parm3d_aniso(dat)
+
+    real(psb_dpk_), intent(in) :: dat(:)
+
+    epsilon = dat(1)
+    theta   = dat(2)
+
+  end subroutine pde_set_parm3d_aniso
+
+  !
+
+  ! Diffusion tensor coefficients
+  !
+  !        [ k11  k12  k13 ]
+  !  K  =  [ k12  k22  k23 ]
+  !        [ k13  k23  k33 ]
+  !
+  ! corresponding to a rotation by theta in the x-y plane
+  ! of the diagonal tensor diag(epsilon,1,1).
+
+  !
+
+  function k11(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: k11
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    k11 = epsilon*cos(theta)**2 + sin(theta)**2
+
+  end function k11
+
+  function k22(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: k22
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    k22 = epsilon*sin(theta)**2 + cos(theta)**2
+
+  end function k22
+
+  function k33(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: k33
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    k33 = done
+
+  end function k33
+
+  function k12(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: k12
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    k12 = (epsilon-done)*sin(theta)*cos(theta)
+
+  end function k12
+
+  function k13(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: k13
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    k13 = dzero
+
+  end function k13
+
+  function k23(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: k23
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    k23 = dzero
+
+  end function k23
+
+  !
+
+  ! Right-hand side / boundary data
+
+  !
+
+  function g_aniso(x,y,z)
+
+    implicit none
+
+    real(psb_dpk_) :: g_aniso
+    real(psb_dpk_), intent(in) :: x,y,z
+
+    g_aniso = dzero
+
+    if (x == done) then
+
+      g_aniso = done
+
+    else if (x == dzero) then
+
+      g_aniso = done
+
+    end if
+
+  end function g_aniso
+
+end module amg_d_pde3d_aniso_mod

samples/advanced/pdegen/runs/amg_pde3d.inp

@@ -1,12 +1,14 @@
 %%%%%%%%%%%  General  arguments % Lines starting with % are ignored.
 CSR                         ! Storage format CSR COO JAD
-0150                        ! IDIM; domain size. Linear system size is IDIM**3
-POISSON                     ! PDECOEFF: POISSON, EXP, BOX, GAUSS
+00070                       ! IDIM; domain size. Linear system size is IDIM**3
+ANISO                       ! PDECOEFF: POISSON, EXP, BOX, GAUSS, ANISO
+100                         ! Epsilon for anistropic PDE (ANISO only)
+0.785398163397448           ! Theta for anistropic PDE (ANISO only) needs to be in radians
 %                           ! Coefficients of the PDE
 CG                          ! Iterative method:
 %                           !    BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES RICHARDSON
 2                           ! ISTOPC
-00500                       ! ITMAX
+00100                       ! ITMAX
 1                           ! ITRACE
 30                          ! IRST (restart for RGMRES and BiCGSTABL)
 1.d-6                       ! EPS
@@ -15,8 +17,8 @@ ML-VBM-VCYC-POLY-L1-D-BJAC       ! Longer descriptive name for preconditioner (u
 ML                        ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY
 %
 %%%%%%%%%%%  First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
-POLY                        ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY  r 1-level, repeats previous.
-6                           ! Number of sweeps for smoother
+RICHARDS                    ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY RICHARDS  r 1-level, repeats previous.
+1                           ! Number of sweeps for smoother
 6                           ! degree for polynomial smoother
 CHEB_4_OPT            ! Polynomial variant
 %  Fields to be added for POLY
@@ -28,12 +30,18 @@ POLY_RHO_EST_POWER     !    POLY_RHO_ESTIMATE             Currently only POLY_RH
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
 NONE                        ! AS prolongation operator: NONE SUM AVG
-L1-JACOBI                     ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF
+MUMPS                       ! Subdomain solver for BJAC/AS: JACOBI GS BGS ILU ILUT MILU MUMPS SLU UMF
 1                           ! Inner solver sweeps (GS and JACOBI)
 LLK                         ! AINV variant
 0                           ! Fill level P for ILU(P) and ILU(T,P)
 1                           ! Inverse Fill level P for INVK
 1.d-4                       ! Threshold T for ILU(T,P)
+1                           ! MUMPS ICNTL(35) BLR activation/format
+0                           ! MUMPS ICNTL(36) BLR variant
+256                         ! MUMPS ICNTL(37) BLR block size
+333                         ! MUMPS ICNTL(38) BLR compression option
+1.d-2                       ! MUMPS CNTL(7) BLR dropping parameter
+2                           ! Force symmetry into MUMPS 
 %%%%%%%%%%%  Second smoother, always ignored for non-ML  %%%%%%%%%%%%%%%%
 NONE                        ! Second (post) smoother, ignored if NONE
 6                           ! Number of sweeps for (post) smoother
@@ -68,7 +76,7 @@ FILTER                      ! Filtering of matrix:  FILTER NOFILTER
 -2                           ! Number of thresholds in vector, next line ignored if <= 0
 0.05 0.025                  ! Thresholds
 %%%%%%%%%%%  Coarse level solver  %%%%%%%%%%%%%%%%
-BJAC                        ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC KRM
+MUMPS                       ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC KRM
 ILU                         ! Coarsest-level subsolver for BJAC: ILU ILUT MILU UMF MUMPS SLU
 DIST                        ! Coarsest-level matrix distribution: DIST  REPL
 1                           ! Coarsest-level fillin P for ILU(P) and ILU(T,P)