mld2p4-2:

mlprec/Makefile
 mlprec/mld_c_as_smoother.f90
 mlprec/mld_c_base_smoother_mod.f90
 mlprec/mld_c_base_solver_mod.f90
 mlprec/mld_c_diag_solver.f90
 mlprec/mld_c_id_solver.f90
 mlprec/mld_c_ilu_solver.f90
 mlprec/mld_c_jac_smoother.f90
 mlprec/mld_c_onelev_mod.f90
 mlprec/mld_c_prec_type.f90
 mlprec/mld_c_slu_solver.f90
 mlprec/mld_c_sludist_solver.f90
 mlprec/mld_c_umf_solver.f90
 mlprec/mld_d_as_smoother.f90
 mlprec/mld_d_base_smoother_mod.f90
 mlprec/mld_d_base_solver_mod.f90
 mlprec/mld_d_diag_solver.f90
 mlprec/mld_d_id_solver.f90
 mlprec/mld_d_ilu_solver.f90
 mlprec/mld_d_jac_smoother.f90
 mlprec/mld_d_onelev_mod.f90
 mlprec/mld_d_prec_type.f90
 mlprec/mld_d_slu_solver.f90
 mlprec/mld_d_sludist_solver.f90
 mlprec/mld_d_umf_solver.f90
 mlprec/mld_s_as_smoother.f90
 mlprec/mld_s_base_smoother_mod.f90
 mlprec/mld_s_base_solver_mod.f90
 mlprec/mld_s_diag_solver.f90
 mlprec/mld_s_id_solver.f90
 mlprec/mld_s_ilu_solver.f90
 mlprec/mld_s_jac_smoother.f90
 mlprec/mld_s_onelev_mod.f90
 mlprec/mld_s_prec_type.f90
 mlprec/mld_s_slu_solver.f90
 mlprec/mld_s_umf_solver.f90
 mlprec/mld_z_as_smoother.f90
 mlprec/mld_z_base_smoother_mod.f90
 mlprec/mld_z_base_solver_mod.f90
 mlprec/mld_z_diag_solver.f90
 mlprec/mld_z_id_solver.f90
 mlprec/mld_z_ilu_solver.f90
 mlprec/mld_z_jac_smoother.f90
 mlprec/mld_z_onelev_mod.f90
 mlprec/mld_z_prec_type.f90
 mlprec/mld_z_slu_solver.f90
 mlprec/mld_z_sludist_solver.f90
 mlprec/mld_z_umf_solver.f90
 tests/fileread/data_input.f90
 tests/fileread/df_sample.f90
 tests/fileread/runs/dfs.inp


Split mld_X_prec_type to reduce memory pressure during compilation. 
Possibly further split will be needed.

mlprec/Makefile

@@ -9,19 +9,23 @@ FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBINCDIR) $(FMFLAG)$(PSBLIBD
 
 DMODOBJS=mld_d_prec_type.o mld_d_prec_mod.o mld_d_move_alloc_mod.o mld_d_ilu_fact_mod.o \
 	mld_d_inner_mod.o mld_d_ilu_solver.o mld_d_diag_solver.o mld_d_jac_smoother.o mld_d_as_smoother.o \
-	mld_d_umf_solver.o mld_d_slu_solver.o mld_d_sludist_solver.o mld_d_id_solver.o 
+	mld_d_umf_solver.o mld_d_slu_solver.o mld_d_sludist_solver.o mld_d_id_solver.o\
+	mld_d_base_solver_mod.o mld_d_base_smoother_mod.o mld_d_onelev_mod.o
 
 SMODOBJS=mld_s_prec_type.o mld_s_prec_mod.o mld_s_move_alloc_mod.o mld_s_ilu_fact_mod.o \
 	mld_s_inner_mod.o mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_jac_smoother.o mld_s_as_smoother.o \
-	mld_s_umf_solver.o mld_s_slu_solver.o mld_s_sludist_solver.o mld_s_id_solver.o 
+	mld_s_umf_solver.o mld_s_slu_solver.o mld_s_sludist_solver.o mld_s_id_solver.o\
+	mld_s_base_solver_mod.o mld_s_base_smoother_mod.o mld_s_onelev_mod.o
 
 ZMODOBJS=mld_z_prec_type.o mld_z_prec_mod.o mld_z_move_alloc_mod.o mld_z_ilu_fact_mod.o \
 	mld_z_inner_mod.o mld_z_ilu_solver.o mld_z_diag_solver.o mld_z_jac_smoother.o mld_z_as_smoother.o \
-	mld_z_umf_solver.o mld_z_slu_solver.o mld_z_sludist_solver.o mld_z_id_solver.o 
+	mld_z_umf_solver.o mld_z_slu_solver.o mld_z_sludist_solver.o mld_z_id_solver.o\
+	mld_z_base_solver_mod.o mld_z_base_smoother_mod.o mld_z_onelev_mod.o
 
 CMODOBJS=mld_c_prec_type.o mld_c_prec_mod.o mld_c_move_alloc_mod.o mld_c_ilu_fact_mod.o \
 	mld_c_inner_mod.o mld_c_ilu_solver.o mld_c_diag_solver.o mld_c_jac_smoother.o mld_c_as_smoother.o \
-	mld_c_umf_solver.o mld_c_slu_solver.o mld_c_sludist_solver.o mld_c_id_solver.o 
+	mld_c_umf_solver.o mld_c_slu_solver.o mld_c_sludist_solver.o mld_c_id_solver.o\
+	mld_c_base_solver_mod.o mld_c_base_smoother_mod.o mld_c_onelev_mod.o
 
 
 
@@ -132,38 +136,63 @@ mld_c_prec_mod.o: mld_c_move_alloc_mod.o
 mld_z_prec_mod.o: mld_z_move_alloc_mod.o 
 
 
+mld_s_prec_type.o: mld_s_onelev_mod.o
+mld_d_prec_type.o: mld_d_onelev_mod.o
+mld_c_prec_type.o: mld_c_onelev_mod.o
+mld_z_prec_type.o: mld_z_onelev_mod.o
+
+mld_s_onelev_mod.o: mld_s_base_smoother_mod.o
+mld_d_onelev_mod.o: mld_d_base_smoother_mod.o
+mld_c_onelev_mod.o: mld_c_base_smoother_mod.o
+mld_z_onelev_mod.o: mld_z_base_smoother_mod.o
+
+mld_s_base_smoother_mod.o: mld_s_base_solver_mod.o
+mld_d_base_smoother_mod.o: mld_d_base_solver_mod.o
+mld_c_base_smoother_mod.o: mld_c_base_solver_mod.o
+mld_z_base_smoother_mod.o: mld_z_base_solver_mod.o
+
+
+mld_s_base_solver_mod.o mld_d_base_solver_mod.o mld_c_base_solver_mod.o mld_z_base_solver_mod.o: mld_base_prec_type.o
+
+
+
+
+mld_d_id_solver.o mld_d_sludist_solver.o mld_d_slu_solver.o \
+mld_d_umf_solver.o mld_d_diag_solver.o mld_d_ilu_solver.o: mld_d_base_solver_mod.o
 
-mld_d_id_solver.o mld_d_sludist_solver.o mld_d_slu_solver.o mld_d_umf_solver.o mld_d_diag_solver.o mld_d_ilu_solver.o: mld_d_prec_type.o 
 mld_d_ilu_fact_mod.o: mld_base_prec_type.o	
 mld_d_ilu_solver.o mld_d_iluk_fact.o: mld_d_ilu_fact_mod.o
-mld_d_as_smoother.o mld_d_jac_smoother.o: mld_d_prec_type.o 
+mld_d_as_smoother.o mld_d_jac_smoother.o: mld_d_base_smoother_mod.o
 mld_d_jac_smoother.o: mld_d_diag_solver.o
 mld_dprecinit.o mld_dprecset.o: mld_d_diag_solver.o mld_d_ilu_solver.o \
 		mld_d_umf_solver.o mld_d_as_smoother.o mld_d_jac_smoother.o \
 	        mld_d_id_solver.o mld_d_slu_solver.o mld_d_sludist_solver.o
 
-mld_s_id_solver.o mld_s_sludist_solver.o mld_s_slu_solver.o mld_s_umf_solver.o mld_s_diag_solver.o mld_s_ilu_solver.o: mld_s_prec_type.o 
+mld_s_id_solver.o mld_s_sludist_solver.o mld_s_slu_solver.o \
+mld_s_umf_solver.o mld_s_diag_solver.o mld_s_ilu_solver.o: mld_s_base_solver_mod.o 
 mld_s_ilu_fact_mod.o: mld_base_prec_type.o	
 mld_s_ilu_solver.o mld_s_iluk_fact.o: mld_s_ilu_fact_mod.o
-mld_s_as_smoother.o mld_s_jac_smoother.o: mld_s_prec_type.o 
+mld_s_as_smoother.o mld_s_jac_smoother.o: mld_s_base_smoother_mod.o
 mld_s_jac_smoother.o: mld_s_diag_solver.o
 mld_sprecinit.o mld_sprecset.o: mld_s_diag_solver.o mld_s_ilu_solver.o \
 		mld_s_umf_solver.o mld_s_as_smoother.o mld_s_jac_smoother.o \
 	        mld_s_id_solver.o mld_s_slu_solver.o mld_s_sludist_solver.o
 
-mld_z_id_solver.o mld_z_sludist_solver.o mld_z_slu_solver.o mld_z_umf_solver.o mld_z_diag_solver.o mld_z_ilu_solver.o: mld_z_prec_type.o 
+mld_z_id_solver.o mld_z_sludist_solver.o mld_z_slu_solver.o \
+mld_z_umf_solver.o mld_z_diag_solver.o mld_z_ilu_solver.o: mld_z_base_solver_mod.o
 mld_z_ilu_fact_mod.o: mld_base_prec_type.o	
 mld_z_ilu_solver.o mld_z_iluk_fact.o: mld_z_ilu_fact_mod.o
-mld_z_as_smoother.o mld_z_jac_smoother.o: mld_z_prec_type.o 
+mld_z_as_smoother.o mld_z_jac_smoother.o: mld_z_base_smoother_mod.o
 mld_z_jac_smoother.o: mld_z_diag_solver.o
 mld_zprecinit.o mld_zprecset.o: mld_z_diag_solver.o mld_z_ilu_solver.o \
 		mld_z_umf_solver.o mld_z_as_smoother.o mld_z_jac_smoother.o \
 	        mld_z_id_solver.o mld_z_slu_solver.o mld_z_sludist_solver.o
 
-mld_c_id_solver.o mld_c_sludist_solver.o mld_c_slu_solver.o mld_c_umf_solver.o mld_c_diag_solver.o mld_c_ilu_solver.o: mld_c_prec_type.o 
+mld_c_id_solver.o mld_c_sludist_solver.o mld_c_slu_solver.o \
+mld_c_umf_solver.o mld_c_diag_solver.o mld_c_ilu_solver.o: mld_c_base_solver_mod.o
 mld_c_ilu_fact_mod.o: mld_base_prec_type.o	
 mld_c_ilu_solver.o mld_c_iluk_fact.o: mld_c_ilu_fact_mod.o
-mld_c_as_smoother.o mld_c_jac_smoother.o: mld_c_prec_type.o 
+mld_c_as_smoother.o mld_c_jac_smoother.o: mld_c_base_smoother_mod.o
 mld_c_jac_smoother.o: mld_c_diag_solver.o
 mld_cprecinit.o mld_cprecset.o: mld_c_diag_solver.o mld_c_ilu_solver.o \
 		mld_c_umf_solver.o mld_c_as_smoother.o mld_c_jac_smoother.o \

mlprec/mld_c_as_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_c_as_smoother
 
-  use mld_c_prec_type
+  use mld_c_base_smoother_mod
   
   type, extends(mld_c_base_smoother_type) :: mld_c_as_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_c_base_smoother_mod.f90

@@ -0,0 +1,612 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_c_base_smoother_mod.f90
+!
+! Module: mld_c_base_smoother_mod
+!
+!  This module defines: 
+!  - the mld_c_base_smoother_type data structure containing the
+!    smoother  and related data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the smoother is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_c_base_smoother_mod
+
+  use mld_c_base_solver_mod
+  !
+  !
+  ! 
+  ! Type: mld_T_base_smoother_type.
+  ! 
+  !  It holds the smoother a single level. Its only mandatory component is a solver
+  !  object which holds a local solver; this decoupling allows to have the same solver
+  !  e.g ILU to work with Jacobi with multiple sweeps as well as with any AS variant.
+  !
+  !  type  mld_T_base_smoother_type
+  !    class(mld_T_base_solver_type), allocatable :: sv
+  !  end type mld_T_base_smoother_type
+  !
+  !   Methods:  
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the solver object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  ! 
+
+  type  mld_c_base_smoother_type
+    class(mld_c_base_solver_type), allocatable :: sv
+  contains
+    procedure, pass(sm) :: check => c_base_smoother_check
+    procedure, pass(sm) :: dump  => c_base_smoother_dmp
+    procedure, pass(sm) :: build => c_base_smoother_bld
+    procedure, pass(sm) :: apply_v => c_base_smoother_apply_vect
+    procedure, pass(sm) :: apply_a => c_base_smoother_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sm) :: free  => c_base_smoother_free
+    procedure, pass(sm) :: seti  => c_base_smoother_seti
+    procedure, pass(sm) :: setc  => c_base_smoother_setc
+    procedure, pass(sm) :: setr  => c_base_smoother_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sm) :: default => c_base_smoother_default
+    procedure, pass(sm) :: descr =>   c_base_smoother_descr
+    procedure, pass(sm) :: sizeof =>  c_base_smoother_sizeof
+    procedure, pass(sm) :: get_nzeros => c_base_smoother_get_nzeros
+  end type mld_c_base_smoother_type
+
+
+  private :: c_base_smoother_bld,   c_base_smoother_apply, &
+       &  c_base_smoother_free,  c_base_smoother_seti, &
+       &  c_base_smoother_setc,  c_base_smoother_setr,&
+       &  c_base_smoother_descr, c_base_smoother_sizeof, &
+       &  c_base_smoother_default, c_base_smoother_check, &
+       &  c_base_smoother_dmp, c_base_smoother_apply_vect, &
+       &  c_base_smoother_get_nzeros
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function c_base_smoother_get_nzeros(sm) result(val)
+    implicit none 
+    class(mld_c_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(sm%sv)) &
+         &  val =  sm%sv%get_nzeros()
+  end function c_base_smoother_get_nzeros
+
+  function c_base_smoother_sizeof(sm) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_)                    :: val
+    integer             :: i
+    
+    val = 0
+    if (allocated(sm%sv)) then 
+      val = sm%sv%sizeof()
+    end if
+
+    return
+  end function c_base_smoother_sizeof
+
+
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! This basic version just applies the local solver, whatever that
+  ! is. 
+  !
+
+  subroutine c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)             :: desc_data
+    class(mld_c_base_smoother_type), intent(in) :: sm
+    complex(psb_spk_),intent(inout)                :: x(:)
+    complex(psb_spk_),intent(inout)                :: y(:)
+    complex(psb_spk_),intent(in)                   :: alpha,beta
+    character(len=1),intent(in)                 :: trans
+    integer, intent(in)                         :: sweeps
+    complex(psb_spk_),target, intent(inout)        :: work(:)
+    integer, intent(out)                        :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine c_base_smoother_apply
+
+  subroutine c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
+       &  trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)                :: desc_data
+    class(mld_c_base_smoother_type), intent(inout) :: sm
+    type(psb_c_vect_type),intent(inout)            :: x
+    type(psb_c_vect_type),intent(inout)            :: y
+    complex(psb_spk_),intent(in)                      :: alpha,beta
+    character(len=1),intent(in)                    :: trans
+    integer, intent(in)                            :: sweeps
+    complex(psb_spk_),target, intent(inout)           :: work(:)
+    integer, intent(out)                           :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine c_base_smoother_apply_vect
+
+  !
+  ! Check:
+  ! 1. Check that we do have a solver object
+  ! 2. Call its check method
+  !
+
+  subroutine c_base_smoother_check(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(inout) :: sm 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_check
+
+  !
+  ! Set methods: the come in multiple versions according
+  ! to whether we are setting with integer, real or character
+  ! input.
+  ! The basic rule is: if the input refers to a parameter
+  ! of the smoother, use it, otherwise pass it to the
+  ! solver object for further processing.
+  ! Since there are no parameters in the base smoother
+  ! we just pass everything to the solver object. 
+  !
+  subroutine c_base_smoother_seti(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    integer, intent(in)                            :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_seti
+
+  subroutine c_base_smoother_setc(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_setc
+  
+  subroutine c_base_smoother_setr(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    real(psb_spk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_smoother_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_setr
+
+
+  !
+  ! Build method.
+  ! At base level we only have to pass data to the inner solver. 
+  ! AMOLD/VMOLD allow to have any relevant sparse matrix or vector
+  ! to be stored in a given format. This is essential e.g.
+  ! when dealing  with GPUs. 
+  !
+  subroutine c_base_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_cspmat_type), intent(in), target      :: a
+    Type(psb_desc_type), Intent(in)                :: desc_a 
+    class(mld_c_base_smoother_type), intent(inout) :: sm 
+    character, intent(in)                          :: upd
+    integer, intent(out)                           :: info
+    class(psb_c_base_sparse_mat), intent(in), optional :: amold
+    class(psb_c_base_vect_type), intent(in), optional  :: vmold
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold)
+    else
+      info = 1121
+      call psb_errpush(info,name)
+    endif
+    if (info /= psb_success_) goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_bld
+
+  !
+  ! Free method (aka destructor).
+  ! In most cases we could do without; however
+  ! for cases where there are data objects allocated outside
+  ! of the Fortran RTE we need to free them explicitly.
+  !
+  ! Even in that case, we could do without this if FINAL
+  ! subroutines were supported, which is not the case
+  ! in GNU Fortran up to 4.7. 
+  !
+  subroutine c_base_smoother_free(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(inout) :: sm
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_free'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    
+    if (allocated(sm%sv)) then 
+      call sm%sv%free(info)
+    end if
+    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
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_free
+
+  !
+  ! Print a description
+  !
+
+  subroutine c_base_smoother_descr(sm,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(in) :: sm
+    integer, intent(out)                        :: info
+    integer, intent(in), optional               :: iout
+    logical, intent(in), optional               :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_c_base_smoother_descr'
+    integer :: iout_
+    logical      :: coarse_
+
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (present(coarse)) then 
+      coarse_ = coarse
+    else
+      coarse_ = .false.
+    end if
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+
+    if (.not.coarse_) &
+         &  write(iout_,*) 'Base smoother with local solver'
+    if (allocated(sm%sv)) then 
+      call sm%sv%descr(info,iout,coarse)
+      if (info /= psb_success_) then 
+        info = psb_err_from_subroutine_ 
+        call psb_errpush(info,name,a_err='Local solver')
+        goto 9999
+      end if
+    end if
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_smoother_descr
+
+  !
+  ! Dump 
+  ! to file, for debugging purposes.
+  !
+  subroutine c_base_smoother_dmp(sm,ictxt,level,info,prefix,head,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_c_base_smoother_type), intent(in) :: sm
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: smoother_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_smth_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(smoother)) then 
+      smoother_ = smoother
+    else
+      smoother_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the smoother
+    if (allocated(sm%sv)) &
+         & call sm%sv%dump(ictxt,level,info,solver=solver)
+
+  end subroutine c_base_smoother_dmp
+
+  !
+  ! Set sensible defaults.
+  ! To be called immediately after allocation
+  !
+  subroutine c_base_smoother_default(sm) 
+    implicit none 
+    ! Arguments
+    class(mld_c_base_smoother_type), intent(inout) :: sm
+    ! Do nothing for base version
+
+    if (allocated(sm%sv)) call sm%sv%default()
+
+    return
+  end subroutine c_base_smoother_default
+
+  
+
+end module mld_c_base_smoother_mod

mlprec/mld_c_base_solver_mod.f90

@@ -0,0 +1,497 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_c_base_solver_mod.f90
+!
+! Module: mld_c_base_solver_mod
+!
+!  This module defines: 
+!  - the mld_c_base_solver_type data structure containing the
+!    basic solver type acting on a subdomain
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the solver is correctly defined;
+!  - printing a	description of the solver;
+!  - deallocating the data structure.  
+!
+
+module mld_c_base_solver_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_c_vect_type, psb_c_base_vect_type
+  !
+  ! 
+  ! Type: mld_T_base_solver_type.
+  ! 
+  !  It holds the local solver; it has no mandatory components. 
+  !
+  !  type  mld_T_base_solver_type
+  !  end type mld_T_base_solver_type
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  !
+
+  type mld_c_base_solver_type
+  contains
+    procedure, pass(sv) :: check => c_base_solver_check
+    procedure, pass(sv) :: dump  => c_base_solver_dmp
+    procedure, pass(sv) :: build => c_base_solver_bld
+    procedure, pass(sv) :: apply_v => c_base_solver_apply_vect
+    procedure, pass(sv) :: apply_a => c_base_solver_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sv) :: free  => c_base_solver_free
+    procedure, pass(sv) :: seti  => c_base_solver_seti
+    procedure, pass(sv) :: setc  => c_base_solver_setc
+    procedure, pass(sv) :: setr  => c_base_solver_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sv) :: default => c_base_solver_default
+    procedure, pass(sv) :: descr   => c_base_solver_descr
+    procedure, pass(sv) :: sizeof  => c_base_solver_sizeof
+    procedure, pass(sv) :: get_nzeros => c_base_solver_get_nzeros
+  end type mld_c_base_solver_type
+
+  private :: c_base_solver_bld,  c_base_solver_apply, &
+       &  c_base_solver_free,    c_base_solver_seti, &
+       &  c_base_solver_setc,    c_base_solver_setr, &
+       &  c_base_solver_descr,   c_base_solver_sizeof, &
+       &  c_base_solver_default, c_base_solver_check,&
+       &  c_base_solver_dmp, c_base_solver_apply_vect, &
+       &  c_base_solver_get_nzeros
+
+
+
+contains
+  !
+  ! Function returning the size of the data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function c_base_solver_sizeof(sv) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_c_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_)                  :: val
+    integer             :: i
+    val = 0
+
+    return
+  end function c_base_solver_sizeof
+
+  function c_base_solver_get_nzeros(sv) result(val)
+    implicit none 
+    class(mld_c_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+  end function c_base_solver_get_nzeros
+
+  
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  ! Question: would it make sense to transform the base version into
+  ! the ID version, i.e. "base_solver" is the identity operator? 
+  ! 
+
+  subroutine c_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)           :: desc_data
+    class(mld_c_base_solver_type), intent(in) :: sv
+    complex(psb_spk_),intent(inout)              :: x(:)
+    complex(psb_spk_),intent(inout)              :: y(:)
+    complex(psb_spk_),intent(in)                 :: alpha,beta
+    character(len=1),intent(in)               :: trans
+    complex(psb_spk_),target, intent(inout)      :: work(:)
+    integer, intent(out)                      :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine c_base_solver_apply
+
+  subroutine c_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)              :: desc_data
+    class(mld_c_base_solver_type), intent(inout) :: sv
+    type(psb_c_vect_type),intent(inout)          :: x
+    type(psb_c_vect_type),intent(inout)          :: y
+    complex(psb_spk_),intent(in)                    :: alpha,beta
+    character(len=1),intent(in)                  :: trans
+    complex(psb_spk_),target, intent(inout)         :: work(:)
+    integer, intent(out)                         :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine c_base_solver_apply_vect
+
+
+  !
+  ! Build
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  !
+  subroutine c_base_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_cspmat_type), intent(in), target           :: a
+    Type(psb_desc_type), Intent(in)                     :: desc_a 
+    class(mld_c_base_solver_type), intent(inout)        :: sv
+    character, intent(in)                               :: upd
+    integer, intent(out)                                :: info
+    type(psb_cspmat_type), intent(in), target, optional :: b
+    class(psb_c_base_sparse_mat), intent(in), optional  :: amold
+    class(psb_c_base_vect_type), intent(in), optional   :: vmold
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_solver_bld
+
+  subroutine c_base_solver_check(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_solver_check
+
+  !
+  ! Set.
+  ! The base version does nothing; the principle is that
+  ! SET acts on what is known, and delegates what is unknown.
+  ! Since we are at the bottom of the hierarchy, there's no one
+  ! to delegate, so we do nothing. 
+  !
+  subroutine c_base_solver_seti(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_seti'
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine c_base_solver_seti
+
+  subroutine c_base_solver_setc(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv
+    integer, intent(in)                          :: what 
+    character(len=*), intent(in)                 :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act, ival 
+    character(len=20) :: name='d_base_solver_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call sv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_solver_setc
+  
+  subroutine c_base_solver_setr(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    real(psb_spk_), intent(in)                   :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_setr'
+
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine c_base_solver_setr
+
+  !
+  ! Free
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. IS THIS CORRECT? I suspect it would be better
+  ! to be silent here, to cover reallocation. 
+  !
+  subroutine c_base_solver_free(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_free'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_solver_free
+
+  subroutine c_base_solver_descr(sv,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_c_base_solver_type), intent(in) :: sv
+    integer, intent(out)                      :: info
+    integer, intent(in), optional             :: iout
+    logical, intent(in), optional             :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_c_base_solver_descr'
+    integer      :: iout_
+
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_solver_descr
+
+  !
+  ! Dump. For debugging purposes. 
+  !
+  subroutine c_base_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_c_base_solver_type), intent(in) :: sv
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: solver_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_slv_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(solver)) then 
+      solver_ = solver
+    else
+      solver_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the solver
+
+  end subroutine c_base_solver_dmp
+
+  subroutine c_base_solver_default(sv) 
+    implicit none 
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv
+    ! Do nothing for base version
+
+    return
+  end subroutine c_base_solver_default
+
+
+
+end module mld_c_base_solver_mod

mlprec/mld_c_diag_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_c_diag_solver
 
-  use mld_c_prec_type
+  use mld_c_base_solver_mod
 
   type, extends(mld_c_base_solver_type) :: mld_c_diag_solver_type
     type(psb_c_vect_type), allocatable :: dv

mlprec/mld_c_id_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_c_id_solver
 
-  use mld_c_prec_type
+  use mld_c_base_solver_mod
 
   type, extends(mld_c_base_solver_type) :: mld_c_id_solver_type
   contains

mlprec/mld_c_ilu_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_c_ilu_solver
 
-  use mld_c_prec_type
+  use mld_c_base_solver_mod
   use mld_c_ilu_fact_mod
 
   type, extends(mld_c_base_solver_type) :: mld_c_ilu_solver_type

mlprec/mld_c_jac_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_c_jac_smoother
 
-  use mld_c_prec_type
+  use mld_c_base_smoother_mod
 
   type, extends(mld_c_base_smoother_type) :: mld_c_jac_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_c_onelev_mod.f90

@@ -0,0 +1,666 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_c_onelev_mod.f90
+!
+! Module: mld_c_onelev_mod
+!
+!  This module defines: 
+!  - the mld_c_onelev_type data structure containing one level
+!    of a multilevel  preconditioner and related
+!    data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the preconditioner is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_c_onelev_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_c_vect_type, psb_c_base_vect_type
+  use mld_c_base_smoother_mod
+  !
+  !
+  ! Type: mld_Tonelev_type.
+  !
+  !  It is the data type containing the necessary items for the	current
+  !  level (essentially, the smoother, the current-level matrix
+  !  and the restriction and prolongation operators).
+  !
+  !  type mld_Tonelev_type
+  !    class(mld_T_base_smoother_type), allocatable :: sm
+  !    type(mld_RTml_parms)            :: parms 
+  !    type(psb_Tspmat_type)           :: ac
+  !    type(psb_Tesc_type)             :: desc_ac
+  !    type(psb_Tspmat_type), pointer  :: base_a    => null() 
+  !    type(psb_Tesc_type), pointer    :: base_desc => null() 
+  !    type(psb_Tlinmap_type)          :: map
+  !  end type mld_Tonelev_type
+  !
+  !  Note that psb_Tpk denotes the kind of the real data type to be chosen
+  !  according to single/double precision version of MLD2P4.
+  !
+  !   sm           -  class(mld_T_base_smoother_type), allocatable
+  !                   The current level preconditioner (aka smoother).
+  !   parms        -  type(mld_RTml_parms)
+  !                   The parameters defining the multilevel strategy.
+  !   ac           -  The local part of the current-level matrix, built by
+  !                   coarsening the previous-level matrix.
+  !   desc_ac      -  type(psb_desc_type).
+  !                   The communication descriptor associated to the matrix
+  !                   stored in ac.
+  !   base_a       -  type(psb_Tspmat_type), pointer.
+  !                   Pointer (really a pointer!) to the local part of the current 
+  !                   matrix (so we have a unified treatment of residuals).
+  !                   We need this to avoid passing explicitly the current matrix
+  !                   to the routine which applies the preconditioner.
+  !   base_desc    -  type(psb_desc_type), pointer.
+  !                   Pointer to the communication descriptor associated to the
+  !                   matrix pointed by base_a.
+  !   map          -  Stores the maps (restriction and prolongation) between the
+  !                   vector spaces associated to the index spaces of the previous
+  !                   and current levels.
+  !
+  !   Methods:  
+  !     Most methods follow the encapsulation hierarchy: they take whatever action
+  !     is appropriate for the current object, then call the corresponding method for
+  !     the contained object.
+  !     As an example: the descr() method prints out a description of the
+  !     level. It starts by invoking the descr() method of the parms object,
+  !     then calls the descr() method of the smoother object. 
+  !
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  type mld_conelev_type
+    class(mld_c_base_smoother_type), allocatable :: sm
+    type(mld_sml_parms)             :: parms 
+    type(psb_cspmat_type)           :: ac
+    type(psb_desc_type)             :: desc_ac
+    type(psb_cspmat_type), pointer  :: base_a    => null() 
+    type(psb_desc_type), pointer    :: base_desc => null() 
+    type(psb_clinmap_type)          :: map
+  contains
+    procedure, pass(lv) :: descr   => c_base_onelev_descr
+    procedure, pass(lv) :: default => c_base_onelev_default
+    procedure, pass(lv) :: free    => c_base_onelev_free
+    procedure, pass(lv) :: nullify => c_base_onelev_nullify
+    procedure, pass(lv) :: check => c_base_onelev_check
+    procedure, pass(lv) :: dump  => c_base_onelev_dump
+    procedure, pass(lv) :: seti  => c_base_onelev_seti
+    procedure, pass(lv) :: setr  => c_base_onelev_setr
+    procedure, pass(lv) :: setc  => c_base_onelev_setc
+    generic, public     :: set   => seti, setr, setc
+    procedure, pass(lv) :: sizeof => c_base_onelev_sizeof
+    procedure, pass(lv) :: get_nzeros => c_base_onelev_get_nzeros
+  end type mld_conelev_type
+
+  private :: c_base_onelev_seti, c_base_onelev_setc, &
+       &  c_base_onelev_setr, c_base_onelev_check, &
+       &  c_base_onelev_default, c_base_onelev_dump, &
+       &  c_base_onelev_descr,  c_base_onelev_sizeof, &
+       &  c_base_onelev_free,  c_base_onelev_nullify,&
+       &  c_base_onelev_get_nzeros
+
+
+  interface mld_nullify_onelevprec
+    module procedure  mld_nullify_d_onelevprec
+  end interface
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function c_base_onelev_get_nzeros(lv) result(val)
+    implicit none 
+    class(mld_conelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(lv%sm)) &
+         &  val =  lv%sm%get_nzeros()
+  end function c_base_onelev_get_nzeros
+
+  function c_base_onelev_sizeof(lv) result(val)
+    implicit none 
+    class(mld_conelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    
+    val = 0
+    val = val + lv%desc_ac%sizeof()
+    val = val + lv%ac%sizeof()
+    val = val + lv%map%sizeof() 
+    if (allocated(lv%sm))  val = val + lv%sm%sizeof()
+  end function c_base_onelev_sizeof
+
+
+  !
+  ! Subroutine: mld_file_onelev_descr
+  ! Version: complex
+  !
+  !  This routine prints a description of the preconditioner to the standard 
+  !  output or to a file. It must be called after the preconditioner has been
+  !  built by mld_precbld.
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tprec_type), input.
+  !             The preconditioner data structure to be printed out.
+  !  info    -  integer, output.
+  !             error code.
+  !  iout    -  integer, input, optional.
+  !             The id of the file where the preconditioner description
+  !             will be printed. If iout is not present, then the standard
+  !             output is condidered.
+  !
+  subroutine c_base_onelev_descr(lv,il,nl,info,iout)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_conelev_type), intent(in) :: lv
+    integer, intent(in)                 :: il,nl
+    integer, intent(out)                :: info
+    integer, intent(in), optional       :: iout
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_c_base_onelev_descr'
+    integer      :: iout_
+    logical      :: coarse
+
+
+    call psb_erractionsave(err_act)
+
+
+    coarse = (il==nl)
+
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+    
+    write(iout_,*) 
+    if (il == 2) then 
+      call lv%parms%mldescr(iout_,info)
+      write(iout_,*) 
+    end if
+
+    if (coarse)  then 
+      write(iout_,*) ' Level ',il,' (coarsest)'
+    else
+      write(iout_,*) ' Level ',il
+    end if
+
+    call lv%parms%descr(iout_,info,coarse=coarse)
+        
+    if (nl > 1) then 
+      if (allocated(lv%map%naggr)) then
+        write(iout_,*) '  Size of coarse matrix: ', &
+             &  sum(lv%map%naggr(:))
+        write(iout_,*) '  Sizes of aggregates: ', &
+             &  lv%map%naggr(:)
+      end if
+    end if
+
+    if (coarse.and.allocated(lv%sm)) &
+         & call lv%sm%descr(info,iout=iout_,coarse=coarse)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_onelev_descr
+
+
+  !
+  ! Subroutines: mld_T_onelev_precfree
+  ! Version: complex
+  !
+  !  These routines deallocate the mld_Tonelev_type 
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tonelev_type), input.
+  !             The data structure to be deallocated.
+  !  info    -  integer, output.
+  !             error code.
+  !
+  subroutine c_base_onelev_free(lv,info)
+    use psb_base_mod
+    implicit none 
+
+    class(mld_conelev_type), intent(inout) :: lv
+    integer, intent(out)                :: info
+    integer :: i
+
+    info = psb_success_
+
+    ! We might just deallocate the top level array, except 
+    ! that there may be inner objects containing C pointers,
+    ! e.g.  UMFPACK, SLU or CUDA stuff.
+    ! We really need FINALs. 
+    call lv%sm%free(info)
+    
+    call lv%ac%free()
+    if (psb_is_ok_desc(lv%desc_ac)) &
+         & call psb_cdfree(lv%desc_ac,info)
+    call lv%map%free(info)
+
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_a) 
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_desc) 
+
+    call lv%nullify()
+
+  end subroutine c_base_onelev_free
+
+
+  subroutine c_base_onelev_nullify(lv)
+    implicit none 
+
+    class(mld_conelev_type), intent(inout) :: lv
+
+    nullify(lv%base_a) 
+    nullify(lv%base_desc) 
+
+  end subroutine c_base_onelev_nullify
+
+
+  subroutine mld_nullify_d_onelevprec(p)
+    implicit none 
+
+    type(mld_conelev_type), intent(inout) :: p
+
+    nullify(p%base_a) 
+    nullify(p%base_desc) 
+
+  end subroutine mld_nullify_d_onelevprec
+
+  !
+  ! Onelevel checks. 
+  ! The number of Jacobi sweeps to be applied is not
+  ! tied to the Jacobi smoother: logically, you have
+  ! a smoother and you can choose to apply it any number
+  ! of times you like. 
+  !
+  subroutine c_base_onelev_check(lv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_conelev_type), intent(inout) :: lv 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    call mld_check_def(lv%parms%sweeps,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_pre,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_post,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+
+    
+    if (allocated(lv%sm)) then 
+      call lv%sm%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_onelev_check
+
+  !
+  ! Multilevel defaults: 
+  !  multiplicative vs. additive ML framework;
+  !  Smoothed decoupled aggregation with zero threshold; 
+  !  distributed coarse matrix;
+  !  damping omega  computed with the max-norm estimate of the
+  !  dominant eigenvalue;
+  !  two-sided smoothing (i.e. V-cycle) with 1 smoothing sweep;
+  !
+
+  subroutine c_base_onelev_default(lv)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_conelev_type), intent(inout) :: lv 
+
+    lv%parms%sweeps          = 1
+    lv%parms%sweeps_pre      = 1
+    lv%parms%sweeps_post     = 1
+    lv%parms%ml_type         = mld_mult_ml_
+    lv%parms%aggr_alg        = mld_dec_aggr_
+    lv%parms%aggr_kind       = mld_smooth_prol_
+    lv%parms%coarse_mat      = mld_distr_mat_
+    lv%parms%smoother_pos    = mld_twoside_smooth_
+    lv%parms%aggr_omega_alg  = mld_eig_est_
+    lv%parms%aggr_eig        = mld_max_norm_
+    lv%parms%aggr_filter     = mld_no_filter_mat_
+    lv%parms%aggr_omega_val  = szero
+    lv%parms%aggr_thresh     = szero
+    
+    if (allocated(lv%sm)) call lv%sm%default()
+
+    return
+
+  end subroutine c_base_onelev_default
+
+  ! 
+  ! Set routines:
+  ! Parameters belonging here are:
+  !   Number of smoothing sweeps;
+  !   Smoother position;
+  !   Aggregation related parameters
+  !   Record request on coarse level solver,
+  !    for checks on solver vs. smoother nomenclature
+  !    reconciliation. 
+  ! 
+  subroutine c_base_onelev_seti(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_conelev_type), intent(inout) :: lv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    select case (what) 
+
+    case (mld_smoother_sweeps_)
+      lv%parms%sweeps      = val
+      lv%parms%sweeps_pre  = val
+      lv%parms%sweeps_post = val
+
+    case (mld_smoother_sweeps_pre_)
+      lv%parms%sweeps_pre  = val
+
+    case (mld_smoother_sweeps_post_)
+      lv%parms%sweeps_post = val
+
+    case (mld_ml_type_)
+      lv%parms%ml_type       = val
+
+    case (mld_aggr_alg_)
+      lv%parms%aggr_alg      = val
+
+    case (mld_aggr_kind_)
+      lv%parms%aggr_kind     = val
+
+    case (mld_coarse_mat_)
+      lv%parms%coarse_mat    = val
+
+    case (mld_smoother_pos_)
+      lv%parms%smoother_pos  = val
+
+    case (mld_aggr_omega_alg_)
+      lv%parms%aggr_omega_alg= val
+
+    case (mld_aggr_eig_)
+      lv%parms%aggr_eig      = val
+
+    case (mld_aggr_filter_)
+      lv%parms%aggr_filter   = val
+
+    case (mld_coarse_solve_)
+      lv%parms%coarse_solve    = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_onelev_seti
+
+  subroutine c_base_onelev_setc(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_conelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_setc'
+    integer :: ival 
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call lv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_onelev_setc
+  
+  subroutine c_base_onelev_setr(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_conelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    real(psb_spk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_onelev_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+    
+    select case (what) 
+
+    case (mld_aggr_omega_val_)
+      lv%parms%aggr_omega_val= val
+   
+    case (mld_aggr_thresh_)
+      lv%parms%aggr_thresh   = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine c_base_onelev_setr
+
+  !
+  ! Dump on file: can be fine-tuned to include the (aggregated) matrix
+  ! as well as smoother and solver. 
+  !
+  subroutine c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_conelev_type), intent(in) :: lv
+    integer, intent(in)              :: level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: ac, rp, smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: ac_, rp_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_lev_d"
+    end if
+
+    if (associated(lv%base_desc)) then 
+      icontxt = lv%base_desc%get_context()
+      call psb_info(icontxt,iam,np)
+    else 
+      icontxt = -1 
+      iam     = -1
+    end if
+    if (present(ac)) then 
+      ac_ = ac
+    else
+      ac_ = .false. 
+    end if
+    if (present(rp)) then 
+      rp_ = rp
+    else
+      rp_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    if (level >= 2) then 
+      if (ac_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%ac%print(fname,head=head)
+      end if
+      if (rp_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_X2Y%print(fname,head=head)
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_Y2X%print(fname,head=head)
+      end if
+    end if
+    if (allocated(lv%sm)) &
+         & call lv%sm%dump(icontxt,level,info,smoother=smoother,solver=solver)
+
+  end subroutine c_base_onelev_dump
+
+
+end module mld_c_onelev_mod

mlprec/mld_c_slu_solver.f90

@@ -46,7 +46,7 @@
 module mld_c_slu_solver
 
   use iso_c_binding
-  use mld_c_prec_type
+  use mld_c_base_solver_mod
 
   type, extends(mld_c_base_solver_type) :: mld_c_slu_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_c_sludist_solver.f90

@@ -46,7 +46,7 @@
 module mld_c_sludist_solver
 
   use iso_c_binding
-  use mld_c_prec_type
+  use mld_c_base_solver_mod
 
   type, extends(mld_c_base_solver_type) :: mld_c_sludist_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_c_umf_solver.f90

@@ -46,7 +46,7 @@
 module mld_c_umf_solver
 
   use iso_c_binding
-  use mld_c_prec_type
+  use mld_c_base_solver_mod
 
   type, extends(mld_c_base_solver_type) :: mld_c_umf_solver_type
     type(c_ptr)                 :: symbolic=c_null_ptr, numeric=c_null_ptr

mlprec/mld_d_as_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_d_as_smoother
 
-  use mld_d_prec_type
+  use mld_d_base_smoother_mod
   
   type, extends(mld_d_base_smoother_type) :: mld_d_as_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_d_base_smoother_mod.f90

@@ -0,0 +1,612 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_d_base_smoother_mod.f90
+!
+! Module: mld_d_base_smoother_mod
+!
+!  This module defines: 
+!  - the mld_d_base_smoother_type data structure containing the
+!    smoother  and related data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the smoother is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_d_base_smoother_mod
+
+  use mld_d_base_solver_mod
+  !
+  !
+  ! 
+  ! Type: mld_T_base_smoother_type.
+  ! 
+  !  It holds the smoother a single level. Its only mandatory component is a solver
+  !  object which holds a local solver; this decoupling allows to have the same solver
+  !  e.g ILU to work with Jacobi with multiple sweeps as well as with any AS variant.
+  !
+  !  type  mld_T_base_smoother_type
+  !    class(mld_T_base_solver_type), allocatable :: sv
+  !  end type mld_T_base_smoother_type
+  !
+  !   Methods:  
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the solver object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  ! 
+
+  type  mld_d_base_smoother_type
+    class(mld_d_base_solver_type), allocatable :: sv
+  contains
+    procedure, pass(sm) :: check => d_base_smoother_check
+    procedure, pass(sm) :: dump  => d_base_smoother_dmp
+    procedure, pass(sm) :: build => d_base_smoother_bld
+    procedure, pass(sm) :: apply_v => d_base_smoother_apply_vect
+    procedure, pass(sm) :: apply_a => d_base_smoother_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sm) :: free  => d_base_smoother_free
+    procedure, pass(sm) :: seti  => d_base_smoother_seti
+    procedure, pass(sm) :: setc  => d_base_smoother_setc
+    procedure, pass(sm) :: setr  => d_base_smoother_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sm) :: default => d_base_smoother_default
+    procedure, pass(sm) :: descr =>   d_base_smoother_descr
+    procedure, pass(sm) :: sizeof =>  d_base_smoother_sizeof
+    procedure, pass(sm) :: get_nzeros => d_base_smoother_get_nzeros
+  end type mld_d_base_smoother_type
+
+
+  private :: d_base_smoother_bld,   d_base_smoother_apply, &
+       &  d_base_smoother_free,  d_base_smoother_seti, &
+       &  d_base_smoother_setc,  d_base_smoother_setr,&
+       &  d_base_smoother_descr, d_base_smoother_sizeof, &
+       &  d_base_smoother_default, d_base_smoother_check, &
+       &  d_base_smoother_dmp, d_base_smoother_apply_vect, &
+       &  d_base_smoother_get_nzeros
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function d_base_smoother_get_nzeros(sm) result(val)
+    implicit none 
+    class(mld_d_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(sm%sv)) &
+         &  val =  sm%sv%get_nzeros()
+  end function d_base_smoother_get_nzeros
+
+  function d_base_smoother_sizeof(sm) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_)                    :: val
+    integer             :: i
+    
+    val = 0
+    if (allocated(sm%sv)) then 
+      val = sm%sv%sizeof()
+    end if
+
+    return
+  end function d_base_smoother_sizeof
+
+
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! This basic version just applies the local solver, whatever that
+  ! is. 
+  !
+
+  subroutine d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)             :: desc_data
+    class(mld_d_base_smoother_type), intent(in) :: 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, intent(in)                         :: sweeps
+    real(psb_dpk_),target, intent(inout)        :: work(:)
+    integer, intent(out)                        :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine d_base_smoother_apply
+
+  subroutine d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
+       &  trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)                :: desc_data
+    class(mld_d_base_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, intent(in)                            :: sweeps
+    real(psb_dpk_),target, intent(inout)           :: work(:)
+    integer, intent(out)                           :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine d_base_smoother_apply_vect
+
+  !
+  ! Check:
+  ! 1. Check that we do have a solver object
+  ! 2. Call its check method
+  !
+
+  subroutine d_base_smoother_check(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(inout) :: sm 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_check
+
+  !
+  ! Set methods: the come in multiple versions according
+  ! to whether we are setting with integer, real or character
+  ! input.
+  ! The basic rule is: if the input refers to a parameter
+  ! of the smoother, use it, otherwise pass it to the
+  ! solver object for further processing.
+  ! Since there are no parameters in the base smoother
+  ! we just pass everything to the solver object. 
+  !
+  subroutine d_base_smoother_seti(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    integer, intent(in)                            :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_seti
+
+  subroutine d_base_smoother_setc(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_setc
+  
+  subroutine d_base_smoother_setr(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    real(psb_dpk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_smoother_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_setr
+
+
+  !
+  ! Build method.
+  ! At base level we only have to pass data to the inner solver. 
+  ! AMOLD/VMOLD allow to have any relevant sparse matrix or vector
+  ! to be stored in a given format. This is essential e.g.
+  ! when dealing  with GPUs. 
+  !
+  subroutine d_base_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_dspmat_type), intent(in), target      :: a
+    Type(psb_desc_type), Intent(in)                :: desc_a 
+    class(mld_d_base_smoother_type), intent(inout) :: sm 
+    character, intent(in)                          :: upd
+    integer, intent(out)                           :: info
+    class(psb_d_base_sparse_mat), intent(in), optional :: amold
+    class(psb_d_base_vect_type), intent(in), optional  :: vmold
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold)
+    else
+      info = 1121
+      call psb_errpush(info,name)
+    endif
+    if (info /= psb_success_) goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_bld
+
+  !
+  ! Free method (aka destructor).
+  ! In most cases we could do without; however
+  ! for cases where there are data objects allocated outside
+  ! of the Fortran RTE we need to free them explicitly.
+  !
+  ! Even in that case, we could do without this if FINAL
+  ! subroutines were supported, which is not the case
+  ! in GNU Fortran up to 4.7. 
+  !
+  subroutine d_base_smoother_free(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(inout) :: sm
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_free'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    
+    if (allocated(sm%sv)) then 
+      call sm%sv%free(info)
+    end if
+    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
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_free
+
+  !
+  ! Print a description
+  !
+
+  subroutine d_base_smoother_descr(sm,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(in) :: sm
+    integer, intent(out)                        :: info
+    integer, intent(in), optional               :: iout
+    logical, intent(in), optional               :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_d_base_smoother_descr'
+    integer :: iout_
+    logical      :: coarse_
+
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (present(coarse)) then 
+      coarse_ = coarse
+    else
+      coarse_ = .false.
+    end if
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+
+    if (.not.coarse_) &
+         &  write(iout_,*) 'Base smoother with local solver'
+    if (allocated(sm%sv)) then 
+      call sm%sv%descr(info,iout,coarse)
+      if (info /= psb_success_) then 
+        info = psb_err_from_subroutine_ 
+        call psb_errpush(info,name,a_err='Local solver')
+        goto 9999
+      end if
+    end if
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_smoother_descr
+
+  !
+  ! Dump 
+  ! to file, for debugging purposes.
+  !
+  subroutine d_base_smoother_dmp(sm,ictxt,level,info,prefix,head,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_d_base_smoother_type), intent(in) :: sm
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: smoother_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_smth_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(smoother)) then 
+      smoother_ = smoother
+    else
+      smoother_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the smoother
+    if (allocated(sm%sv)) &
+         & call sm%sv%dump(ictxt,level,info,solver=solver)
+
+  end subroutine d_base_smoother_dmp
+
+  !
+  ! Set sensible defaults.
+  ! To be called immediately after allocation
+  !
+  subroutine d_base_smoother_default(sm) 
+    implicit none 
+    ! Arguments
+    class(mld_d_base_smoother_type), intent(inout) :: sm
+    ! Do nothing for base version
+
+    if (allocated(sm%sv)) call sm%sv%default()
+
+    return
+  end subroutine d_base_smoother_default
+
+  
+
+end module mld_d_base_smoother_mod

mlprec/mld_d_base_solver_mod.f90

@@ -0,0 +1,497 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_d_base_solver_mod.f90
+!
+! Module: mld_d_base_solver_mod
+!
+!  This module defines: 
+!  - the mld_d_base_solver_type data structure containing the
+!    basic solver type acting on a subdomain
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the solver is correctly defined;
+!  - printing a	description of the solver;
+!  - deallocating the data structure.  
+!
+
+module mld_d_base_solver_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_d_vect_type, psb_d_base_vect_type
+  !
+  ! 
+  ! Type: mld_T_base_solver_type.
+  ! 
+  !  It holds the local solver; it has no mandatory components. 
+  !
+  !  type  mld_T_base_solver_type
+  !  end type mld_T_base_solver_type
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  !
+
+  type mld_d_base_solver_type
+  contains
+    procedure, pass(sv) :: check => d_base_solver_check
+    procedure, pass(sv) :: dump  => d_base_solver_dmp
+    procedure, pass(sv) :: build => d_base_solver_bld
+    procedure, pass(sv) :: apply_v => d_base_solver_apply_vect
+    procedure, pass(sv) :: apply_a => d_base_solver_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sv) :: free  => d_base_solver_free
+    procedure, pass(sv) :: seti  => d_base_solver_seti
+    procedure, pass(sv) :: setc  => d_base_solver_setc
+    procedure, pass(sv) :: setr  => d_base_solver_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sv) :: default => d_base_solver_default
+    procedure, pass(sv) :: descr   => d_base_solver_descr
+    procedure, pass(sv) :: sizeof  => d_base_solver_sizeof
+    procedure, pass(sv) :: get_nzeros => d_base_solver_get_nzeros
+  end type mld_d_base_solver_type
+
+  private :: d_base_solver_bld,  d_base_solver_apply, &
+       &  d_base_solver_free,    d_base_solver_seti, &
+       &  d_base_solver_setc,    d_base_solver_setr, &
+       &  d_base_solver_descr,   d_base_solver_sizeof, &
+       &  d_base_solver_default, d_base_solver_check,&
+       &  d_base_solver_dmp, d_base_solver_apply_vect, &
+       &  d_base_solver_get_nzeros
+
+
+
+contains
+  !
+  ! Function returning the size of the data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function d_base_solver_sizeof(sv) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_d_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_)                  :: val
+    integer             :: i
+    val = 0
+
+    return
+  end function d_base_solver_sizeof
+
+  function d_base_solver_get_nzeros(sv) result(val)
+    implicit none 
+    class(mld_d_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+  end function d_base_solver_get_nzeros
+
+  
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  ! Question: would it make sense to transform the base version into
+  ! the ID version, i.e. "base_solver" is the identity operator? 
+  ! 
+
+  subroutine d_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)           :: desc_data
+    class(mld_d_base_solver_type), intent(in) :: sv
+    real(psb_dpk_),intent(inout)              :: x(:)
+    real(psb_dpk_),intent(inout)              :: y(:)
+    real(psb_dpk_),intent(in)                 :: alpha,beta
+    character(len=1),intent(in)               :: trans
+    real(psb_dpk_),target, intent(inout)      :: work(:)
+    integer, intent(out)                      :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine d_base_solver_apply
+
+  subroutine d_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)              :: desc_data
+    class(mld_d_base_solver_type), intent(inout) :: sv
+    type(psb_d_vect_type),intent(inout)          :: x
+    type(psb_d_vect_type),intent(inout)          :: y
+    real(psb_dpk_),intent(in)                    :: alpha,beta
+    character(len=1),intent(in)                  :: trans
+    real(psb_dpk_),target, intent(inout)         :: work(:)
+    integer, intent(out)                         :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine d_base_solver_apply_vect
+
+
+  !
+  ! Build
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  !
+  subroutine d_base_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_dspmat_type), intent(in), target           :: a
+    Type(psb_desc_type), Intent(in)                     :: desc_a 
+    class(mld_d_base_solver_type), intent(inout)        :: sv
+    character, intent(in)                               :: upd
+    integer, intent(out)                                :: info
+    type(psb_dspmat_type), intent(in), target, optional :: b
+    class(psb_d_base_sparse_mat), intent(in), optional  :: amold
+    class(psb_d_base_vect_type), intent(in), optional   :: vmold
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_solver_bld
+
+  subroutine d_base_solver_check(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_solver_check
+
+  !
+  ! Set.
+  ! The base version does nothing; the principle is that
+  ! SET acts on what is known, and delegates what is unknown.
+  ! Since we are at the bottom of the hierarchy, there's no one
+  ! to delegate, so we do nothing. 
+  !
+  subroutine d_base_solver_seti(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_seti'
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine d_base_solver_seti
+
+  subroutine d_base_solver_setc(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv
+    integer, intent(in)                          :: what 
+    character(len=*), intent(in)                 :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act, ival 
+    character(len=20) :: name='d_base_solver_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call sv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_solver_setc
+  
+  subroutine d_base_solver_setr(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    real(psb_dpk_), intent(in)                   :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_setr'
+
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine d_base_solver_setr
+
+  !
+  ! Free
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. IS THIS CORRECT? I suspect it would be better
+  ! to be silent here, to cover reallocation. 
+  !
+  subroutine d_base_solver_free(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_free'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_solver_free
+
+  subroutine d_base_solver_descr(sv,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_d_base_solver_type), intent(in) :: sv
+    integer, intent(out)                      :: info
+    integer, intent(in), optional             :: iout
+    logical, intent(in), optional             :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_d_base_solver_descr'
+    integer      :: iout_
+
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_solver_descr
+
+  !
+  ! Dump. For debugging purposes. 
+  !
+  subroutine d_base_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_d_base_solver_type), intent(in) :: sv
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: solver_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_slv_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(solver)) then 
+      solver_ = solver
+    else
+      solver_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the solver
+
+  end subroutine d_base_solver_dmp
+
+  subroutine d_base_solver_default(sv) 
+    implicit none 
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv
+    ! Do nothing for base version
+
+    return
+  end subroutine d_base_solver_default
+
+
+
+end module mld_d_base_solver_mod

mlprec/mld_d_diag_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_d_diag_solver
 
-  use mld_d_prec_type
+  use mld_d_base_solver_mod
 
   type, extends(mld_d_base_solver_type) :: mld_d_diag_solver_type
     type(psb_d_vect_type), allocatable :: dv

mlprec/mld_d_id_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_d_id_solver
 
-  use mld_d_prec_type
+  use mld_d_base_solver_mod
 
   type, extends(mld_d_base_solver_type) :: mld_d_id_solver_type
   contains

mlprec/mld_d_ilu_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_d_ilu_solver
 
-  use mld_d_prec_type
+  use mld_d_base_solver_mod
   use mld_d_ilu_fact_mod
 
   type, extends(mld_d_base_solver_type) :: mld_d_ilu_solver_type

mlprec/mld_d_jac_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_d_jac_smoother
 
-  use mld_d_prec_type
+  use mld_d_base_smoother_mod
 
   type, extends(mld_d_base_smoother_type) :: mld_d_jac_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_d_onelev_mod.f90

@@ -0,0 +1,666 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_d_onelev_mod.f90
+!
+! Module: mld_d_onelev_mod
+!
+!  This module defines: 
+!  - the mld_d_onelev_type data structure containing one level
+!    of a multilevel  preconditioner and related
+!    data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the preconditioner is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_d_onelev_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_d_vect_type, psb_d_base_vect_type
+  use mld_d_base_smoother_mod
+  !
+  !
+  ! Type: mld_Tonelev_type.
+  !
+  !  It is the data type containing the necessary items for the	current
+  !  level (essentially, the smoother, the current-level matrix
+  !  and the restriction and prolongation operators).
+  !
+  !  type mld_Tonelev_type
+  !    class(mld_T_base_smoother_type), allocatable :: sm
+  !    type(mld_RTml_parms)            :: parms 
+  !    type(psb_Tspmat_type)           :: ac
+  !    type(psb_Tesc_type)             :: desc_ac
+  !    type(psb_Tspmat_type), pointer  :: base_a    => null() 
+  !    type(psb_Tesc_type), pointer    :: base_desc => null() 
+  !    type(psb_Tlinmap_type)          :: map
+  !  end type mld_Tonelev_type
+  !
+  !  Note that psb_Tpk denotes the kind of the real data type to be chosen
+  !  according to single/double precision version of MLD2P4.
+  !
+  !   sm           -  class(mld_T_base_smoother_type), allocatable
+  !                   The current level preconditioner (aka smoother).
+  !   parms        -  type(mld_RTml_parms)
+  !                   The parameters defining the multilevel strategy.
+  !   ac           -  The local part of the current-level matrix, built by
+  !                   coarsening the previous-level matrix.
+  !   desc_ac      -  type(psb_desc_type).
+  !                   The communication descriptor associated to the matrix
+  !                   stored in ac.
+  !   base_a       -  type(psb_Tspmat_type), pointer.
+  !                   Pointer (really a pointer!) to the local part of the current 
+  !                   matrix (so we have a unified treatment of residuals).
+  !                   We need this to avoid passing explicitly the current matrix
+  !                   to the routine which applies the preconditioner.
+  !   base_desc    -  type(psb_desc_type), pointer.
+  !                   Pointer to the communication descriptor associated to the
+  !                   matrix pointed by base_a.
+  !   map          -  Stores the maps (restriction and prolongation) between the
+  !                   vector spaces associated to the index spaces of the previous
+  !                   and current levels.
+  !
+  !   Methods:  
+  !     Most methods follow the encapsulation hierarchy: they take whatever action
+  !     is appropriate for the current object, then call the corresponding method for
+  !     the contained object.
+  !     As an example: the descr() method prints out a description of the
+  !     level. It starts by invoking the descr() method of the parms object,
+  !     then calls the descr() method of the smoother object. 
+  !
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  type mld_donelev_type
+    class(mld_d_base_smoother_type), allocatable :: sm
+    type(mld_dml_parms)             :: parms 
+    type(psb_dspmat_type)           :: ac
+    type(psb_desc_type)             :: desc_ac
+    type(psb_dspmat_type), pointer  :: base_a    => null() 
+    type(psb_desc_type), pointer    :: base_desc => null() 
+    type(psb_dlinmap_type)          :: map
+  contains
+    procedure, pass(lv) :: descr   => d_base_onelev_descr
+    procedure, pass(lv) :: default => d_base_onelev_default
+    procedure, pass(lv) :: free    => d_base_onelev_free
+    procedure, pass(lv) :: nullify => d_base_onelev_nullify
+    procedure, pass(lv) :: check => d_base_onelev_check
+    procedure, pass(lv) :: dump  => d_base_onelev_dump
+    procedure, pass(lv) :: seti  => d_base_onelev_seti
+    procedure, pass(lv) :: setr  => d_base_onelev_setr
+    procedure, pass(lv) :: setc  => d_base_onelev_setc
+    generic, public     :: set   => seti, setr, setc
+    procedure, pass(lv) :: sizeof => d_base_onelev_sizeof
+    procedure, pass(lv) :: get_nzeros => d_base_onelev_get_nzeros
+  end type mld_donelev_type
+
+  private :: d_base_onelev_seti, d_base_onelev_setc, &
+       &  d_base_onelev_setr, d_base_onelev_check, &
+       &  d_base_onelev_default, d_base_onelev_dump, &
+       &  d_base_onelev_descr,  d_base_onelev_sizeof, &
+       &  d_base_onelev_free,  d_base_onelev_nullify,&
+       &  d_base_onelev_get_nzeros
+
+
+  interface mld_nullify_onelevprec
+    module procedure  mld_nullify_d_onelevprec
+  end interface
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function d_base_onelev_get_nzeros(lv) result(val)
+    implicit none 
+    class(mld_donelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(lv%sm)) &
+         &  val =  lv%sm%get_nzeros()
+  end function d_base_onelev_get_nzeros
+
+  function d_base_onelev_sizeof(lv) result(val)
+    implicit none 
+    class(mld_donelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    
+    val = 0
+    val = val + lv%desc_ac%sizeof()
+    val = val + lv%ac%sizeof()
+    val = val + lv%map%sizeof() 
+    if (allocated(lv%sm))  val = val + lv%sm%sizeof()
+  end function d_base_onelev_sizeof
+
+
+  !
+  ! Subroutine: mld_file_onelev_descr
+  ! Version: real
+  !
+  !  This routine prints a description of the preconditioner to the standard 
+  !  output or to a file. It must be called after the preconditioner has been
+  !  built by mld_precbld.
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tprec_type), input.
+  !             The preconditioner data structure to be printed out.
+  !  info    -  integer, output.
+  !             error code.
+  !  iout    -  integer, input, optional.
+  !             The id of the file where the preconditioner description
+  !             will be printed. If iout is not present, then the standard
+  !             output is condidered.
+  !
+  subroutine d_base_onelev_descr(lv,il,nl,info,iout)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_donelev_type), intent(in) :: lv
+    integer, intent(in)                 :: il,nl
+    integer, intent(out)                :: info
+    integer, intent(in), optional       :: iout
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_d_base_onelev_descr'
+    integer      :: iout_
+    logical      :: coarse
+
+
+    call psb_erractionsave(err_act)
+
+
+    coarse = (il==nl)
+
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+    
+    write(iout_,*) 
+    if (il == 2) then 
+      call lv%parms%mldescr(iout_,info)
+      write(iout_,*) 
+    end if
+
+    if (coarse)  then 
+      write(iout_,*) ' Level ',il,' (coarsest)'
+    else
+      write(iout_,*) ' Level ',il
+    end if
+
+    call lv%parms%descr(iout_,info,coarse=coarse)
+        
+    if (nl > 1) then 
+      if (allocated(lv%map%naggr)) then
+        write(iout_,*) '  Size of coarse matrix: ', &
+             &  sum(lv%map%naggr(:))
+        write(iout_,*) '  Sizes of aggregates: ', &
+             &  lv%map%naggr(:)
+      end if
+    end if
+
+    if (coarse.and.allocated(lv%sm)) &
+         & call lv%sm%descr(info,iout=iout_,coarse=coarse)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_onelev_descr
+
+
+  !
+  ! Subroutines: mld_T_onelev_precfree
+  ! Version: real
+  !
+  !  These routines deallocate the mld_Tonelev_type 
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tonelev_type), input.
+  !             The data structure to be deallocated.
+  !  info    -  integer, output.
+  !             error code.
+  !
+  subroutine d_base_onelev_free(lv,info)
+    use psb_base_mod
+    implicit none 
+
+    class(mld_donelev_type), intent(inout) :: lv
+    integer, intent(out)                :: info
+    integer :: i
+
+    info = psb_success_
+
+    ! We might just deallocate the top level array, except 
+    ! that there may be inner objects containing C pointers,
+    ! e.g.  UMFPACK, SLU or CUDA stuff.
+    ! We really need FINALs. 
+    call lv%sm%free(info)
+    
+    call lv%ac%free()
+    if (psb_is_ok_desc(lv%desc_ac)) &
+         & call psb_cdfree(lv%desc_ac,info)
+    call lv%map%free(info)
+
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_a) 
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_desc) 
+
+    call lv%nullify()
+
+  end subroutine d_base_onelev_free
+
+
+  subroutine d_base_onelev_nullify(lv)
+    implicit none 
+
+    class(mld_donelev_type), intent(inout) :: lv
+
+    nullify(lv%base_a) 
+    nullify(lv%base_desc) 
+
+  end subroutine d_base_onelev_nullify
+
+
+  subroutine mld_nullify_d_onelevprec(p)
+    implicit none 
+
+    type(mld_donelev_type), intent(inout) :: p
+
+    nullify(p%base_a) 
+    nullify(p%base_desc) 
+
+  end subroutine mld_nullify_d_onelevprec
+
+  !
+  ! Onelevel checks. 
+  ! The number of Jacobi sweeps to be applied is not
+  ! tied to the Jacobi smoother: logically, you have
+  ! a smoother and you can choose to apply it any number
+  ! of times you like. 
+  !
+  subroutine d_base_onelev_check(lv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_donelev_type), intent(inout) :: lv 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    call mld_check_def(lv%parms%sweeps,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_pre,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_post,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+
+    
+    if (allocated(lv%sm)) then 
+      call lv%sm%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_onelev_check
+
+  !
+  ! Multilevel defaults: 
+  !  multiplicative vs. additive ML framework;
+  !  Smoothed decoupled aggregation with zero threshold; 
+  !  distributed coarse matrix;
+  !  damping omega  computed with the max-norm estimate of the
+  !  dominant eigenvalue;
+  !  two-sided smoothing (i.e. V-cycle) with 1 smoothing sweep;
+  !
+
+  subroutine d_base_onelev_default(lv)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_donelev_type), intent(inout) :: lv 
+
+    lv%parms%sweeps          = 1
+    lv%parms%sweeps_pre      = 1
+    lv%parms%sweeps_post     = 1
+    lv%parms%ml_type         = mld_mult_ml_
+    lv%parms%aggr_alg        = mld_dec_aggr_
+    lv%parms%aggr_kind       = mld_smooth_prol_
+    lv%parms%coarse_mat      = mld_distr_mat_
+    lv%parms%smoother_pos    = mld_twoside_smooth_
+    lv%parms%aggr_omega_alg  = mld_eig_est_
+    lv%parms%aggr_eig        = mld_max_norm_
+    lv%parms%aggr_filter     = mld_no_filter_mat_
+    lv%parms%aggr_omega_val  = dzero
+    lv%parms%aggr_thresh     = dzero
+    
+    if (allocated(lv%sm)) call lv%sm%default()
+
+    return
+
+  end subroutine d_base_onelev_default
+
+  ! 
+  ! Set routines:
+  ! Parameters belonging here are:
+  !   Number of smoothing sweeps;
+  !   Smoother position;
+  !   Aggregation related parameters
+  !   Record request on coarse level solver,
+  !    for checks on solver vs. smoother nomenclature
+  !    reconciliation. 
+  ! 
+  subroutine d_base_onelev_seti(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_donelev_type), intent(inout) :: lv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    select case (what) 
+
+    case (mld_smoother_sweeps_)
+      lv%parms%sweeps      = val
+      lv%parms%sweeps_pre  = val
+      lv%parms%sweeps_post = val
+
+    case (mld_smoother_sweeps_pre_)
+      lv%parms%sweeps_pre  = val
+
+    case (mld_smoother_sweeps_post_)
+      lv%parms%sweeps_post = val
+
+    case (mld_ml_type_)
+      lv%parms%ml_type       = val
+
+    case (mld_aggr_alg_)
+      lv%parms%aggr_alg      = val
+
+    case (mld_aggr_kind_)
+      lv%parms%aggr_kind     = val
+
+    case (mld_coarse_mat_)
+      lv%parms%coarse_mat    = val
+
+    case (mld_smoother_pos_)
+      lv%parms%smoother_pos  = val
+
+    case (mld_aggr_omega_alg_)
+      lv%parms%aggr_omega_alg= val
+
+    case (mld_aggr_eig_)
+      lv%parms%aggr_eig      = val
+
+    case (mld_aggr_filter_)
+      lv%parms%aggr_filter   = val
+
+    case (mld_coarse_solve_)
+      lv%parms%coarse_solve    = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_onelev_seti
+
+  subroutine d_base_onelev_setc(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_donelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_setc'
+    integer :: ival 
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call lv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_onelev_setc
+  
+  subroutine d_base_onelev_setr(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_donelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    real(psb_dpk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_onelev_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+    
+    select case (what) 
+
+    case (mld_aggr_omega_val_)
+      lv%parms%aggr_omega_val= val
+   
+    case (mld_aggr_thresh_)
+      lv%parms%aggr_thresh   = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine d_base_onelev_setr
+
+  !
+  ! Dump on file: can be fine-tuned to include the (aggregated) matrix
+  ! as well as smoother and solver. 
+  !
+  subroutine d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_donelev_type), intent(in) :: lv
+    integer, intent(in)              :: level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: ac, rp, smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: ac_, rp_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_lev_d"
+    end if
+
+    if (associated(lv%base_desc)) then 
+      icontxt = lv%base_desc%get_context()
+      call psb_info(icontxt,iam,np)
+    else 
+      icontxt = -1 
+      iam     = -1
+    end if
+    if (present(ac)) then 
+      ac_ = ac
+    else
+      ac_ = .false. 
+    end if
+    if (present(rp)) then 
+      rp_ = rp
+    else
+      rp_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    if (level >= 2) then 
+      if (ac_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%ac%print(fname,head=head)
+      end if
+      if (rp_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_X2Y%print(fname,head=head)
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_Y2X%print(fname,head=head)
+      end if
+    end if
+    if (allocated(lv%sm)) &
+         & call lv%sm%dump(icontxt,level,info,smoother=smoother,solver=solver)
+
+  end subroutine d_base_onelev_dump
+
+
+end module mld_d_onelev_mod

mlprec/mld_d_slu_solver.f90

@@ -46,7 +46,7 @@
 module mld_d_slu_solver
 
   use iso_c_binding
-  use mld_d_prec_type
+  use mld_d_base_solver_mod
 
   type, extends(mld_d_base_solver_type) :: mld_d_slu_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_d_sludist_solver.f90

@@ -46,7 +46,7 @@
 module mld_d_sludist_solver
 
   use iso_c_binding
-  use mld_d_prec_type
+  use mld_d_base_solver_mod
 
   type, extends(mld_d_base_solver_type) :: mld_d_sludist_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_d_umf_solver.f90

@@ -46,7 +46,7 @@
 module mld_d_umf_solver
 
   use iso_c_binding
-  use mld_d_prec_type
+  use mld_d_base_solver_mod
 
   type, extends(mld_d_base_solver_type) :: mld_d_umf_solver_type
     type(c_ptr)                 :: symbolic=c_null_ptr, numeric=c_null_ptr

mlprec/mld_s_as_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_s_as_smoother
 
-  use mld_s_prec_type
+  use mld_s_base_smoother_mod  
   
   type, extends(mld_s_base_smoother_type) :: mld_s_as_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_s_base_smoother_mod.f90

@@ -0,0 +1,612 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_s_base_smoother_mod.f90
+!
+! Module: mld_s_base_smoother_mod
+!
+!  This module defines: 
+!  - the mld_s_base_smoother_type data structure containing the
+!    smoother  and related data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the smoother is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_s_base_smoother_mod
+
+  use mld_s_base_solver_mod
+  !
+  !
+  ! 
+  ! Type: mld_T_base_smoother_type.
+  ! 
+  !  It holds the smoother a single level. Its only mandatory component is a solver
+  !  object which holds a local solver; this decoupling allows to have the same solver
+  !  e.g ILU to work with Jacobi with multiple sweeps as well as with any AS variant.
+  !
+  !  type  mld_T_base_smoother_type
+  !    class(mld_T_base_solver_type), allocatable :: sv
+  !  end type mld_T_base_smoother_type
+  !
+  !   Methods:  
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the solver object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  ! 
+
+  type  mld_s_base_smoother_type
+    class(mld_s_base_solver_type), allocatable :: sv
+  contains
+    procedure, pass(sm) :: check => s_base_smoother_check
+    procedure, pass(sm) :: dump  => s_base_smoother_dmp
+    procedure, pass(sm) :: build => s_base_smoother_bld
+    procedure, pass(sm) :: apply_v => s_base_smoother_apply_vect
+    procedure, pass(sm) :: apply_a => s_base_smoother_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sm) :: free  => s_base_smoother_free
+    procedure, pass(sm) :: seti  => s_base_smoother_seti
+    procedure, pass(sm) :: setc  => s_base_smoother_setc
+    procedure, pass(sm) :: setr  => s_base_smoother_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sm) :: default => s_base_smoother_default
+    procedure, pass(sm) :: descr =>   s_base_smoother_descr
+    procedure, pass(sm) :: sizeof =>  s_base_smoother_sizeof
+    procedure, pass(sm) :: get_nzeros => s_base_smoother_get_nzeros
+  end type mld_s_base_smoother_type
+
+
+  private :: s_base_smoother_bld,   s_base_smoother_apply, &
+       &  s_base_smoother_free,  s_base_smoother_seti, &
+       &  s_base_smoother_setc,  s_base_smoother_setr,&
+       &  s_base_smoother_descr, s_base_smoother_sizeof, &
+       &  s_base_smoother_default, s_base_smoother_check, &
+       &  s_base_smoother_dmp, s_base_smoother_apply_vect, &
+       &  s_base_smoother_get_nzeros
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function s_base_smoother_get_nzeros(sm) result(val)
+    implicit none 
+    class(mld_s_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(sm%sv)) &
+         &  val =  sm%sv%get_nzeros()
+  end function s_base_smoother_get_nzeros
+
+  function s_base_smoother_sizeof(sm) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_)                    :: val
+    integer             :: i
+    
+    val = 0
+    if (allocated(sm%sv)) then 
+      val = sm%sv%sizeof()
+    end if
+
+    return
+  end function s_base_smoother_sizeof
+
+
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! This basic version just applies the local solver, whatever that
+  ! is. 
+  !
+
+  subroutine s_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)             :: desc_data
+    class(mld_s_base_smoother_type), intent(in) :: sm
+    real(psb_spk_),intent(inout)                :: x(:)
+    real(psb_spk_),intent(inout)                :: y(:)
+    real(psb_spk_),intent(in)                   :: alpha,beta
+    character(len=1),intent(in)                 :: trans
+    integer, intent(in)                         :: sweeps
+    real(psb_spk_),target, intent(inout)        :: work(:)
+    integer, intent(out)                        :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine s_base_smoother_apply
+
+  subroutine s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
+       &  trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)                :: desc_data
+    class(mld_s_base_smoother_type), intent(inout) :: sm
+    type(psb_s_vect_type),intent(inout)            :: x
+    type(psb_s_vect_type),intent(inout)            :: y
+    real(psb_spk_),intent(in)                      :: alpha,beta
+    character(len=1),intent(in)                    :: trans
+    integer, intent(in)                            :: sweeps
+    real(psb_spk_),target, intent(inout)           :: work(:)
+    integer, intent(out)                           :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine s_base_smoother_apply_vect
+
+  !
+  ! Check:
+  ! 1. Check that we do have a solver object
+  ! 2. Call its check method
+  !
+
+  subroutine s_base_smoother_check(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(inout) :: sm 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_check
+
+  !
+  ! Set methods: the come in multiple versions according
+  ! to whether we are setting with integer, real or character
+  ! input.
+  ! The basic rule is: if the input refers to a parameter
+  ! of the smoother, use it, otherwise pass it to the
+  ! solver object for further processing.
+  ! Since there are no parameters in the base smoother
+  ! we just pass everything to the solver object. 
+  !
+  subroutine s_base_smoother_seti(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    integer, intent(in)                            :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_seti
+
+  subroutine s_base_smoother_setc(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_setc
+  
+  subroutine s_base_smoother_setr(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    real(psb_spk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_smoother_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_setr
+
+
+  !
+  ! Build method.
+  ! At base level we only have to pass data to the inner solver. 
+  ! AMOLD/VMOLD allow to have any relevant sparse matrix or vector
+  ! to be stored in a given format. This is essential e.g.
+  ! when dealing  with GPUs. 
+  !
+  subroutine s_base_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_sspmat_type), intent(in), target      :: a
+    Type(psb_desc_type), Intent(in)                :: desc_a 
+    class(mld_s_base_smoother_type), intent(inout) :: sm 
+    character, intent(in)                          :: upd
+    integer, intent(out)                           :: info
+    class(psb_s_base_sparse_mat), intent(in), optional :: amold
+    class(psb_s_base_vect_type), intent(in), optional  :: vmold
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold)
+    else
+      info = 1121
+      call psb_errpush(info,name)
+    endif
+    if (info /= psb_success_) goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_bld
+
+  !
+  ! Free method (aka destructor).
+  ! In most cases we could do without; however
+  ! for cases where there are data objects allocated outside
+  ! of the Fortran RTE we need to free them explicitly.
+  !
+  ! Even in that case, we could do without this if FINAL
+  ! subroutines were supported, which is not the case
+  ! in GNU Fortran up to 4.7. 
+  !
+  subroutine s_base_smoother_free(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(inout) :: sm
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_free'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    
+    if (allocated(sm%sv)) then 
+      call sm%sv%free(info)
+    end if
+    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
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_free
+
+  !
+  ! Print a description
+  !
+
+  subroutine s_base_smoother_descr(sm,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(in) :: sm
+    integer, intent(out)                        :: info
+    integer, intent(in), optional               :: iout
+    logical, intent(in), optional               :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_s_base_smoother_descr'
+    integer :: iout_
+    logical      :: coarse_
+
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (present(coarse)) then 
+      coarse_ = coarse
+    else
+      coarse_ = .false.
+    end if
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+
+    if (.not.coarse_) &
+         &  write(iout_,*) 'Base smoother with local solver'
+    if (allocated(sm%sv)) then 
+      call sm%sv%descr(info,iout,coarse)
+      if (info /= psb_success_) then 
+        info = psb_err_from_subroutine_ 
+        call psb_errpush(info,name,a_err='Local solver')
+        goto 9999
+      end if
+    end if
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_smoother_descr
+
+  !
+  ! Dump 
+  ! to file, for debugging purposes.
+  !
+  subroutine s_base_smoother_dmp(sm,ictxt,level,info,prefix,head,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_s_base_smoother_type), intent(in) :: sm
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: smoother_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_smth_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(smoother)) then 
+      smoother_ = smoother
+    else
+      smoother_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the smoother
+    if (allocated(sm%sv)) &
+         & call sm%sv%dump(ictxt,level,info,solver=solver)
+
+  end subroutine s_base_smoother_dmp
+
+  !
+  ! Set sensible defaults.
+  ! To be called immediately after allocation
+  !
+  subroutine s_base_smoother_default(sm) 
+    implicit none 
+    ! Arguments
+    class(mld_s_base_smoother_type), intent(inout) :: sm
+    ! Do nothing for base version
+
+    if (allocated(sm%sv)) call sm%sv%default()
+
+    return
+  end subroutine s_base_smoother_default
+
+  
+
+end module mld_s_base_smoother_mod

mlprec/mld_s_base_solver_mod.f90

@@ -0,0 +1,497 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_s_base_solver_mod.f90
+!
+! Module: mld_s_base_solver_mod
+!
+!  This module defines: 
+!  - the mld_s_base_solver_type data structure containing the
+!    basic solver type acting on a subdomain
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the solver is correctly defined;
+!  - printing a	description of the solver;
+!  - deallocating the data structure.  
+!
+
+module mld_s_base_solver_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_s_vect_type, psb_s_base_vect_type
+  !
+  ! 
+  ! Type: mld_T_base_solver_type.
+  ! 
+  !  It holds the local solver; it has no mandatory components. 
+  !
+  !  type  mld_T_base_solver_type
+  !  end type mld_T_base_solver_type
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  !
+
+  type mld_s_base_solver_type
+  contains
+    procedure, pass(sv) :: check => s_base_solver_check
+    procedure, pass(sv) :: dump  => s_base_solver_dmp
+    procedure, pass(sv) :: build => s_base_solver_bld
+    procedure, pass(sv) :: apply_v => s_base_solver_apply_vect
+    procedure, pass(sv) :: apply_a => s_base_solver_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sv) :: free  => s_base_solver_free
+    procedure, pass(sv) :: seti  => s_base_solver_seti
+    procedure, pass(sv) :: setc  => s_base_solver_setc
+    procedure, pass(sv) :: setr  => s_base_solver_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sv) :: default => s_base_solver_default
+    procedure, pass(sv) :: descr   => s_base_solver_descr
+    procedure, pass(sv) :: sizeof  => s_base_solver_sizeof
+    procedure, pass(sv) :: get_nzeros => s_base_solver_get_nzeros
+  end type mld_s_base_solver_type
+
+  private :: s_base_solver_bld,  s_base_solver_apply, &
+       &  s_base_solver_free,    s_base_solver_seti, &
+       &  s_base_solver_setc,    s_base_solver_setr, &
+       &  s_base_solver_descr,   s_base_solver_sizeof, &
+       &  s_base_solver_default, s_base_solver_check,&
+       &  s_base_solver_dmp, s_base_solver_apply_vect, &
+       &  s_base_solver_get_nzeros
+
+
+
+contains
+  !
+  ! Function returning the size of the data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function s_base_solver_sizeof(sv) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_s_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_)                  :: val
+    integer             :: i
+    val = 0
+
+    return
+  end function s_base_solver_sizeof
+
+  function s_base_solver_get_nzeros(sv) result(val)
+    implicit none 
+    class(mld_s_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+  end function s_base_solver_get_nzeros
+
+  
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  ! Question: would it make sense to transform the base version into
+  ! the ID version, i.e. "base_solver" is the identity operator? 
+  ! 
+
+  subroutine s_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)           :: desc_data
+    class(mld_s_base_solver_type), intent(in) :: sv
+    real(psb_spk_),intent(inout)              :: x(:)
+    real(psb_spk_),intent(inout)              :: y(:)
+    real(psb_spk_),intent(in)                 :: alpha,beta
+    character(len=1),intent(in)               :: trans
+    real(psb_spk_),target, intent(inout)      :: work(:)
+    integer, intent(out)                      :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine s_base_solver_apply
+
+  subroutine s_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)              :: desc_data
+    class(mld_s_base_solver_type), intent(inout) :: sv
+    type(psb_s_vect_type),intent(inout)          :: x
+    type(psb_s_vect_type),intent(inout)          :: y
+    real(psb_spk_),intent(in)                    :: alpha,beta
+    character(len=1),intent(in)                  :: trans
+    real(psb_spk_),target, intent(inout)         :: work(:)
+    integer, intent(out)                         :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine s_base_solver_apply_vect
+
+
+  !
+  ! Build
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  !
+  subroutine s_base_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_sspmat_type), intent(in), target           :: a
+    Type(psb_desc_type), Intent(in)                     :: desc_a 
+    class(mld_s_base_solver_type), intent(inout)        :: sv
+    character, intent(in)                               :: upd
+    integer, intent(out)                                :: info
+    type(psb_sspmat_type), intent(in), target, optional :: b
+    class(psb_s_base_sparse_mat), intent(in), optional  :: amold
+    class(psb_s_base_vect_type), intent(in), optional   :: vmold
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_solver_bld
+
+  subroutine s_base_solver_check(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_solver_check
+
+  !
+  ! Set.
+  ! The base version does nothing; the principle is that
+  ! SET acts on what is known, and delegates what is unknown.
+  ! Since we are at the bottom of the hierarchy, there's no one
+  ! to delegate, so we do nothing. 
+  !
+  subroutine s_base_solver_seti(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_seti'
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine s_base_solver_seti
+
+  subroutine s_base_solver_setc(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv
+    integer, intent(in)                          :: what 
+    character(len=*), intent(in)                 :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act, ival 
+    character(len=20) :: name='d_base_solver_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call sv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_solver_setc
+  
+  subroutine s_base_solver_setr(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    real(psb_spk_), intent(in)                   :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_setr'
+
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine s_base_solver_setr
+
+  !
+  ! Free
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. IS THIS CORRECT? I suspect it would be better
+  ! to be silent here, to cover reallocation. 
+  !
+  subroutine s_base_solver_free(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_free'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_solver_free
+
+  subroutine s_base_solver_descr(sv,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_s_base_solver_type), intent(in) :: sv
+    integer, intent(out)                      :: info
+    integer, intent(in), optional             :: iout
+    logical, intent(in), optional             :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_s_base_solver_descr'
+    integer      :: iout_
+
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_solver_descr
+
+  !
+  ! Dump. For debugging purposes. 
+  !
+  subroutine s_base_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_s_base_solver_type), intent(in) :: sv
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: solver_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_slv_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(solver)) then 
+      solver_ = solver
+    else
+      solver_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the solver
+
+  end subroutine s_base_solver_dmp
+
+  subroutine s_base_solver_default(sv) 
+    implicit none 
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv
+    ! Do nothing for base version
+
+    return
+  end subroutine s_base_solver_default
+
+
+
+end module mld_s_base_solver_mod

mlprec/mld_s_diag_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_s_diag_solver
 
-  use mld_s_prec_type
+  use mld_s_base_solver_mod
 
   type, extends(mld_s_base_solver_type) :: mld_s_diag_solver_type
     type(psb_s_vect_type), allocatable :: dv

mlprec/mld_s_id_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_s_id_solver
 
-  use mld_s_prec_type
+  use mld_s_base_solver_mod
 
   type, extends(mld_s_base_solver_type) :: mld_s_id_solver_type
   contains

mlprec/mld_s_ilu_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_s_ilu_solver
 
-  use mld_s_prec_type
+  use mld_s_base_solver_mod
   use mld_s_ilu_fact_mod
 
   type, extends(mld_s_base_solver_type) :: mld_s_ilu_solver_type

mlprec/mld_s_jac_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_s_jac_smoother
 
-  use mld_s_prec_type
+  use mld_s_base_smoother_mod
 
   type, extends(mld_s_base_smoother_type) :: mld_s_jac_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_s_onelev_mod.f90

@@ -0,0 +1,666 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_s_onelev_mod.f90
+!
+! Module: mld_s_onelev_mod
+!
+!  This module defines: 
+!  - the mld_s_onelev_type data structure containing one level
+!    of a multilevel  preconditioner and related
+!    data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the preconditioner is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_s_onelev_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_s_vect_type, psb_s_base_vect_type
+  use mld_s_base_smoother_mod
+  !
+  !
+  ! Type: mld_Tonelev_type.
+  !
+  !  It is the data type containing the necessary items for the	current
+  !  level (essentially, the smoother, the current-level matrix
+  !  and the restriction and prolongation operators).
+  !
+  !  type mld_Tonelev_type
+  !    class(mld_T_base_smoother_type), allocatable :: sm
+  !    type(mld_RTml_parms)            :: parms 
+  !    type(psb_Tspmat_type)           :: ac
+  !    type(psb_Tesc_type)             :: desc_ac
+  !    type(psb_Tspmat_type), pointer  :: base_a    => null() 
+  !    type(psb_Tesc_type), pointer    :: base_desc => null() 
+  !    type(psb_Tlinmap_type)          :: map
+  !  end type mld_Tonelev_type
+  !
+  !  Note that psb_Tpk denotes the kind of the real data type to be chosen
+  !  according to single/double precision version of MLD2P4.
+  !
+  !   sm           -  class(mld_T_base_smoother_type), allocatable
+  !                   The current level preconditioner (aka smoother).
+  !   parms        -  type(mld_RTml_parms)
+  !                   The parameters defining the multilevel strategy.
+  !   ac           -  The local part of the current-level matrix, built by
+  !                   coarsening the previous-level matrix.
+  !   desc_ac      -  type(psb_desc_type).
+  !                   The communication descriptor associated to the matrix
+  !                   stored in ac.
+  !   base_a       -  type(psb_Tspmat_type), pointer.
+  !                   Pointer (really a pointer!) to the local part of the current 
+  !                   matrix (so we have a unified treatment of residuals).
+  !                   We need this to avoid passing explicitly the current matrix
+  !                   to the routine which applies the preconditioner.
+  !   base_desc    -  type(psb_desc_type), pointer.
+  !                   Pointer to the communication descriptor associated to the
+  !                   matrix pointed by base_a.
+  !   map          -  Stores the maps (restriction and prolongation) between the
+  !                   vector spaces associated to the index spaces of the previous
+  !                   and current levels.
+  !
+  !   Methods:  
+  !     Most methods follow the encapsulation hierarchy: they take whatever action
+  !     is appropriate for the current object, then call the corresponding method for
+  !     the contained object.
+  !     As an example: the descr() method prints out a description of the
+  !     level. It starts by invoking the descr() method of the parms object,
+  !     then calls the descr() method of the smoother object. 
+  !
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  type mld_sonelev_type
+    class(mld_s_base_smoother_type), allocatable :: sm
+    type(mld_sml_parms)             :: parms 
+    type(psb_sspmat_type)           :: ac
+    type(psb_desc_type)             :: desc_ac
+    type(psb_sspmat_type), pointer  :: base_a    => null() 
+    type(psb_desc_type), pointer    :: base_desc => null() 
+    type(psb_slinmap_type)          :: map
+  contains
+    procedure, pass(lv) :: descr   => s_base_onelev_descr
+    procedure, pass(lv) :: default => s_base_onelev_default
+    procedure, pass(lv) :: free    => s_base_onelev_free
+    procedure, pass(lv) :: nullify => s_base_onelev_nullify
+    procedure, pass(lv) :: check => s_base_onelev_check
+    procedure, pass(lv) :: dump  => s_base_onelev_dump
+    procedure, pass(lv) :: seti  => s_base_onelev_seti
+    procedure, pass(lv) :: setr  => s_base_onelev_setr
+    procedure, pass(lv) :: setc  => s_base_onelev_setc
+    generic, public     :: set   => seti, setr, setc
+    procedure, pass(lv) :: sizeof => s_base_onelev_sizeof
+    procedure, pass(lv) :: get_nzeros => s_base_onelev_get_nzeros
+  end type mld_sonelev_type
+
+  private :: s_base_onelev_seti, s_base_onelev_setc, &
+       &  s_base_onelev_setr, s_base_onelev_check, &
+       &  s_base_onelev_default, s_base_onelev_dump, &
+       &  s_base_onelev_descr,  s_base_onelev_sizeof, &
+       &  s_base_onelev_free,  s_base_onelev_nullify,&
+       &  s_base_onelev_get_nzeros
+
+
+  interface mld_nullify_onelevprec
+    module procedure  mld_nullify_d_onelevprec
+  end interface
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function s_base_onelev_get_nzeros(lv) result(val)
+    implicit none 
+    class(mld_sonelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(lv%sm)) &
+         &  val =  lv%sm%get_nzeros()
+  end function s_base_onelev_get_nzeros
+
+  function s_base_onelev_sizeof(lv) result(val)
+    implicit none 
+    class(mld_sonelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    
+    val = 0
+    val = val + lv%desc_ac%sizeof()
+    val = val + lv%ac%sizeof()
+    val = val + lv%map%sizeof() 
+    if (allocated(lv%sm))  val = val + lv%sm%sizeof()
+  end function s_base_onelev_sizeof
+
+
+  !
+  ! Subroutine: mld_file_onelev_descr
+  ! Version: real
+  !
+  !  This routine prints a description of the preconditioner to the standard 
+  !  output or to a file. It must be called after the preconditioner has been
+  !  built by mld_precbld.
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tprec_type), input.
+  !             The preconditioner data structure to be printed out.
+  !  info    -  integer, output.
+  !             error code.
+  !  iout    -  integer, input, optional.
+  !             The id of the file where the preconditioner description
+  !             will be printed. If iout is not present, then the standard
+  !             output is condidered.
+  !
+  subroutine s_base_onelev_descr(lv,il,nl,info,iout)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_sonelev_type), intent(in) :: lv
+    integer, intent(in)                 :: il,nl
+    integer, intent(out)                :: info
+    integer, intent(in), optional       :: iout
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_s_base_onelev_descr'
+    integer      :: iout_
+    logical      :: coarse
+
+
+    call psb_erractionsave(err_act)
+
+
+    coarse = (il==nl)
+
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+    
+    write(iout_,*) 
+    if (il == 2) then 
+      call lv%parms%mldescr(iout_,info)
+      write(iout_,*) 
+    end if
+
+    if (coarse)  then 
+      write(iout_,*) ' Level ',il,' (coarsest)'
+    else
+      write(iout_,*) ' Level ',il
+    end if
+
+    call lv%parms%descr(iout_,info,coarse=coarse)
+        
+    if (nl > 1) then 
+      if (allocated(lv%map%naggr)) then
+        write(iout_,*) '  Size of coarse matrix: ', &
+             &  sum(lv%map%naggr(:))
+        write(iout_,*) '  Sizes of aggregates: ', &
+             &  lv%map%naggr(:)
+      end if
+    end if
+
+    if (coarse.and.allocated(lv%sm)) &
+         & call lv%sm%descr(info,iout=iout_,coarse=coarse)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_onelev_descr
+
+
+  !
+  ! Subroutines: mld_T_onelev_precfree
+  ! Version: real
+  !
+  !  These routines deallocate the mld_Tonelev_type 
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tonelev_type), input.
+  !             The data structure to be deallocated.
+  !  info    -  integer, output.
+  !             error code.
+  !
+  subroutine s_base_onelev_free(lv,info)
+    use psb_base_mod
+    implicit none 
+
+    class(mld_sonelev_type), intent(inout) :: lv
+    integer, intent(out)                :: info
+    integer :: i
+
+    info = psb_success_
+
+    ! We might just deallocate the top level array, except 
+    ! that there may be inner objects containing C pointers,
+    ! e.g.  UMFPACK, SLU or CUDA stuff.
+    ! We really need FINALs. 
+    call lv%sm%free(info)
+    
+    call lv%ac%free()
+    if (psb_is_ok_desc(lv%desc_ac)) &
+         & call psb_cdfree(lv%desc_ac,info)
+    call lv%map%free(info)
+
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_a) 
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_desc) 
+
+    call lv%nullify()
+
+  end subroutine s_base_onelev_free
+
+
+  subroutine s_base_onelev_nullify(lv)
+    implicit none 
+
+    class(mld_sonelev_type), intent(inout) :: lv
+
+    nullify(lv%base_a) 
+    nullify(lv%base_desc) 
+
+  end subroutine s_base_onelev_nullify
+
+
+  subroutine mld_nullify_d_onelevprec(p)
+    implicit none 
+
+    type(mld_sonelev_type), intent(inout) :: p
+
+    nullify(p%base_a) 
+    nullify(p%base_desc) 
+
+  end subroutine mld_nullify_d_onelevprec
+
+  !
+  ! Onelevel checks. 
+  ! The number of Jacobi sweeps to be applied is not
+  ! tied to the Jacobi smoother: logically, you have
+  ! a smoother and you can choose to apply it any number
+  ! of times you like. 
+  !
+  subroutine s_base_onelev_check(lv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_sonelev_type), intent(inout) :: lv 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    call mld_check_def(lv%parms%sweeps,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_pre,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_post,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+
+    
+    if (allocated(lv%sm)) then 
+      call lv%sm%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_onelev_check
+
+  !
+  ! Multilevel defaults: 
+  !  multiplicative vs. additive ML framework;
+  !  Smoothed decoupled aggregation with zero threshold; 
+  !  distributed coarse matrix;
+  !  damping omega  computed with the max-norm estimate of the
+  !  dominant eigenvalue;
+  !  two-sided smoothing (i.e. V-cycle) with 1 smoothing sweep;
+  !
+
+  subroutine s_base_onelev_default(lv)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_sonelev_type), intent(inout) :: lv 
+
+    lv%parms%sweeps          = 1
+    lv%parms%sweeps_pre      = 1
+    lv%parms%sweeps_post     = 1
+    lv%parms%ml_type         = mld_mult_ml_
+    lv%parms%aggr_alg        = mld_dec_aggr_
+    lv%parms%aggr_kind       = mld_smooth_prol_
+    lv%parms%coarse_mat      = mld_distr_mat_
+    lv%parms%smoother_pos    = mld_twoside_smooth_
+    lv%parms%aggr_omega_alg  = mld_eig_est_
+    lv%parms%aggr_eig        = mld_max_norm_
+    lv%parms%aggr_filter     = mld_no_filter_mat_
+    lv%parms%aggr_omega_val  = szero
+    lv%parms%aggr_thresh     = szero
+    
+    if (allocated(lv%sm)) call lv%sm%default()
+
+    return
+
+  end subroutine s_base_onelev_default
+
+  ! 
+  ! Set routines:
+  ! Parameters belonging here are:
+  !   Number of smoothing sweeps;
+  !   Smoother position;
+  !   Aggregation related parameters
+  !   Record request on coarse level solver,
+  !    for checks on solver vs. smoother nomenclature
+  !    reconciliation. 
+  ! 
+  subroutine s_base_onelev_seti(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_sonelev_type), intent(inout) :: lv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    select case (what) 
+
+    case (mld_smoother_sweeps_)
+      lv%parms%sweeps      = val
+      lv%parms%sweeps_pre  = val
+      lv%parms%sweeps_post = val
+
+    case (mld_smoother_sweeps_pre_)
+      lv%parms%sweeps_pre  = val
+
+    case (mld_smoother_sweeps_post_)
+      lv%parms%sweeps_post = val
+
+    case (mld_ml_type_)
+      lv%parms%ml_type       = val
+
+    case (mld_aggr_alg_)
+      lv%parms%aggr_alg      = val
+
+    case (mld_aggr_kind_)
+      lv%parms%aggr_kind     = val
+
+    case (mld_coarse_mat_)
+      lv%parms%coarse_mat    = val
+
+    case (mld_smoother_pos_)
+      lv%parms%smoother_pos  = val
+
+    case (mld_aggr_omega_alg_)
+      lv%parms%aggr_omega_alg= val
+
+    case (mld_aggr_eig_)
+      lv%parms%aggr_eig      = val
+
+    case (mld_aggr_filter_)
+      lv%parms%aggr_filter   = val
+
+    case (mld_coarse_solve_)
+      lv%parms%coarse_solve    = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_onelev_seti
+
+  subroutine s_base_onelev_setc(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_sonelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_setc'
+    integer :: ival 
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call lv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_onelev_setc
+  
+  subroutine s_base_onelev_setr(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_sonelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    real(psb_spk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_onelev_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+    
+    select case (what) 
+
+    case (mld_aggr_omega_val_)
+      lv%parms%aggr_omega_val= val
+   
+    case (mld_aggr_thresh_)
+      lv%parms%aggr_thresh   = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine s_base_onelev_setr
+
+  !
+  ! Dump on file: can be fine-tuned to include the (aggregated) matrix
+  ! as well as smoother and solver. 
+  !
+  subroutine s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_sonelev_type), intent(in) :: lv
+    integer, intent(in)              :: level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: ac, rp, smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: ac_, rp_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_lev_d"
+    end if
+
+    if (associated(lv%base_desc)) then 
+      icontxt = lv%base_desc%get_context()
+      call psb_info(icontxt,iam,np)
+    else 
+      icontxt = -1 
+      iam     = -1
+    end if
+    if (present(ac)) then 
+      ac_ = ac
+    else
+      ac_ = .false. 
+    end if
+    if (present(rp)) then 
+      rp_ = rp
+    else
+      rp_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    if (level >= 2) then 
+      if (ac_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%ac%print(fname,head=head)
+      end if
+      if (rp_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_X2Y%print(fname,head=head)
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_Y2X%print(fname,head=head)
+      end if
+    end if
+    if (allocated(lv%sm)) &
+         & call lv%sm%dump(icontxt,level,info,smoother=smoother,solver=solver)
+
+  end subroutine s_base_onelev_dump
+
+
+end module mld_s_onelev_mod

mlprec/mld_s_slu_solver.f90

@@ -46,7 +46,7 @@
 module mld_s_slu_solver
 
   use iso_c_binding
-  use mld_s_prec_type
+  use mld_s_base_solver_mod
 
   type, extends(mld_s_base_solver_type) :: mld_s_slu_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_s_umf_solver.f90

@@ -46,7 +46,7 @@
 module mld_s_umf_solver
 
   use iso_c_binding
-  use mld_s_prec_type
+  use mld_s_base_solver_mod
 
   type, extends(mld_s_base_solver_type) :: mld_s_umf_solver_type
     type(c_ptr)                 :: symbolic=c_null_ptr, numeric=c_null_ptr

mlprec/mld_z_as_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_z_as_smoother
 
-  use mld_z_prec_type
+  use mld_z_base_smoother_mod
   
   type, extends(mld_z_base_smoother_type) :: mld_z_as_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_z_base_smoother_mod.f90

@@ -0,0 +1,612 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_z_base_smoother_mod.f90
+!
+! Module: mld_z_base_smoother_mod
+!
+!  This module defines: 
+!  - the mld_z_base_smoother_type data structure containing the
+!    smoother  and related data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the smoother is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_z_base_smoother_mod
+
+  use mld_z_base_solver_mod
+  !
+  !
+  ! 
+  ! Type: mld_T_base_smoother_type.
+  ! 
+  !  It holds the smoother a single level. Its only mandatory component is a solver
+  !  object which holds a local solver; this decoupling allows to have the same solver
+  !  e.g ILU to work with Jacobi with multiple sweeps as well as with any AS variant.
+  !
+  !  type  mld_T_base_smoother_type
+  !    class(mld_T_base_solver_type), allocatable :: sv
+  !  end type mld_T_base_smoother_type
+  !
+  !   Methods:  
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the solver object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  ! 
+
+  type  mld_z_base_smoother_type
+    class(mld_z_base_solver_type), allocatable :: sv
+  contains
+    procedure, pass(sm) :: check => z_base_smoother_check
+    procedure, pass(sm) :: dump  => z_base_smoother_dmp
+    procedure, pass(sm) :: build => z_base_smoother_bld
+    procedure, pass(sm) :: apply_v => z_base_smoother_apply_vect
+    procedure, pass(sm) :: apply_a => z_base_smoother_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sm) :: free  => z_base_smoother_free
+    procedure, pass(sm) :: seti  => z_base_smoother_seti
+    procedure, pass(sm) :: setc  => z_base_smoother_setc
+    procedure, pass(sm) :: setr  => z_base_smoother_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sm) :: default => z_base_smoother_default
+    procedure, pass(sm) :: descr =>   z_base_smoother_descr
+    procedure, pass(sm) :: sizeof =>  z_base_smoother_sizeof
+    procedure, pass(sm) :: get_nzeros => z_base_smoother_get_nzeros
+  end type mld_z_base_smoother_type
+
+
+  private :: z_base_smoother_bld,   z_base_smoother_apply, &
+       &  z_base_smoother_free,  z_base_smoother_seti, &
+       &  z_base_smoother_setc,  z_base_smoother_setr,&
+       &  z_base_smoother_descr, z_base_smoother_sizeof, &
+       &  z_base_smoother_default, z_base_smoother_check, &
+       &  z_base_smoother_dmp, z_base_smoother_apply_vect, &
+       &  z_base_smoother_get_nzeros
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function z_base_smoother_get_nzeros(sm) result(val)
+    implicit none 
+    class(mld_z_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(sm%sv)) &
+         &  val =  sm%sv%get_nzeros()
+  end function z_base_smoother_get_nzeros
+
+  function z_base_smoother_sizeof(sm) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(in) :: sm
+    integer(psb_long_int_k_)                    :: val
+    integer             :: i
+    
+    val = 0
+    if (allocated(sm%sv)) then 
+      val = sm%sv%sizeof()
+    end if
+
+    return
+  end function z_base_smoother_sizeof
+
+
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! This basic version just applies the local solver, whatever that
+  ! is. 
+  !
+
+  subroutine z_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)             :: desc_data
+    class(mld_z_base_smoother_type), intent(in) :: sm
+    complex(psb_dpk_),intent(inout)                :: x(:)
+    complex(psb_dpk_),intent(inout)                :: y(:)
+    complex(psb_dpk_),intent(in)                   :: alpha,beta
+    character(len=1),intent(in)                 :: trans
+    integer, intent(in)                         :: sweeps
+    complex(psb_dpk_),target, intent(inout)        :: work(:)
+    integer, intent(out)                        :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine z_base_smoother_apply
+
+  subroutine z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
+       &  trans,sweeps,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)                :: desc_data
+    class(mld_z_base_smoother_type), intent(inout) :: sm
+    type(psb_z_vect_type),intent(inout)            :: x
+    type(psb_z_vect_type),intent(inout)            :: y
+    complex(psb_dpk_),intent(in)                      :: alpha,beta
+    character(len=1),intent(in)                    :: trans
+    integer, intent(in)                            :: sweeps
+    complex(psb_dpk_),target, intent(inout)           :: work(:)
+    integer, intent(out)                           :: info
+    
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_apply'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+    else
+      info = 1121
+    endif
+    if (info /= psb_success_) then 
+      call psb_errpush(info,name)
+      goto 9999 
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine z_base_smoother_apply_vect
+
+  !
+  ! Check:
+  ! 1. Check that we do have a solver object
+  ! 2. Call its check method
+  !
+
+  subroutine z_base_smoother_check(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(inout) :: sm 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_check
+
+  !
+  ! Set methods: the come in multiple versions according
+  ! to whether we are setting with integer, real or character
+  ! input.
+  ! The basic rule is: if the input refers to a parameter
+  ! of the smoother, use it, otherwise pass it to the
+  ! solver object for further processing.
+  ! Since there are no parameters in the base smoother
+  ! we just pass everything to the solver object. 
+  !
+  subroutine z_base_smoother_seti(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    integer, intent(in)                            :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_seti
+
+  subroutine z_base_smoother_setc(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_setc
+  
+  subroutine z_base_smoother_setr(sm,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(inout) :: sm 
+    integer, intent(in)                            :: what 
+    real(psb_dpk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_smoother_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+
+    if (allocated(sm%sv)) then 
+      call sm%sv%set(what,val,info)
+    end if
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_setr
+
+
+  !
+  ! Build method.
+  ! At base level we only have to pass data to the inner solver. 
+  ! AMOLD/VMOLD allow to have any relevant sparse matrix or vector
+  ! to be stored in a given format. This is essential e.g.
+  ! when dealing  with GPUs. 
+  !
+  subroutine z_base_smoother_bld(a,desc_a,sm,upd,info,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_zspmat_type), intent(in), target      :: a
+    Type(psb_desc_type), Intent(in)                :: desc_a 
+    class(mld_z_base_smoother_type), intent(inout) :: sm 
+    character, intent(in)                          :: upd
+    integer, intent(out)                           :: info
+    class(psb_z_base_sparse_mat), intent(in), optional :: amold
+    class(psb_z_base_vect_type), intent(in), optional  :: vmold
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+    if (allocated(sm%sv)) then 
+      call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold)
+    else
+      info = 1121
+      call psb_errpush(info,name)
+    endif
+    if (info /= psb_success_) goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_bld
+
+  !
+  ! Free method (aka destructor).
+  ! In most cases we could do without; however
+  ! for cases where there are data objects allocated outside
+  ! of the Fortran RTE we need to free them explicitly.
+  !
+  ! Even in that case, we could do without this if FINAL
+  ! subroutines were supported, which is not the case
+  ! in GNU Fortran up to 4.7. 
+  !
+  subroutine z_base_smoother_free(sm,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(inout) :: sm
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_smoother_free'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+    
+    if (allocated(sm%sv)) then 
+      call sm%sv%free(info)
+    end if
+    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
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_free
+
+  !
+  ! Print a description
+  !
+
+  subroutine z_base_smoother_descr(sm,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(in) :: sm
+    integer, intent(out)                        :: info
+    integer, intent(in), optional               :: iout
+    logical, intent(in), optional               :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_z_base_smoother_descr'
+    integer :: iout_
+    logical      :: coarse_
+
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    if (present(coarse)) then 
+      coarse_ = coarse
+    else
+      coarse_ = .false.
+    end if
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+
+    if (.not.coarse_) &
+         &  write(iout_,*) 'Base smoother with local solver'
+    if (allocated(sm%sv)) then 
+      call sm%sv%descr(info,iout,coarse)
+      if (info /= psb_success_) then 
+        info = psb_err_from_subroutine_ 
+        call psb_errpush(info,name,a_err='Local solver')
+        goto 9999
+      end if
+    end if
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_smoother_descr
+
+  !
+  ! Dump 
+  ! to file, for debugging purposes.
+  !
+  subroutine z_base_smoother_dmp(sm,ictxt,level,info,prefix,head,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_z_base_smoother_type), intent(in) :: sm
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: smoother_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_smth_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(smoother)) then 
+      smoother_ = smoother
+    else
+      smoother_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the smoother
+    if (allocated(sm%sv)) &
+         & call sm%sv%dump(ictxt,level,info,solver=solver)
+
+  end subroutine z_base_smoother_dmp
+
+  !
+  ! Set sensible defaults.
+  ! To be called immediately after allocation
+  !
+  subroutine z_base_smoother_default(sm) 
+    implicit none 
+    ! Arguments
+    class(mld_z_base_smoother_type), intent(inout) :: sm
+    ! Do nothing for base version
+
+    if (allocated(sm%sv)) call sm%sv%default()
+
+    return
+  end subroutine z_base_smoother_default
+
+  
+
+end module mld_z_base_smoother_mod

mlprec/mld_z_base_solver_mod.f90

@@ -0,0 +1,497 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_z_base_solver_mod.f90
+!
+! Module: mld_z_base_solver_mod
+!
+!  This module defines: 
+!  - the mld_z_base_solver_type data structure containing the
+!    basic solver type acting on a subdomain
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the solver is correctly defined;
+!  - printing a	description of the solver;
+!  - deallocating the data structure.  
+!
+
+module mld_z_base_solver_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_z_vect_type, psb_z_base_vect_type
+  !
+  ! 
+  ! Type: mld_T_base_solver_type.
+  ! 
+  !  It holds the local solver; it has no mandatory components. 
+  !
+  !  type  mld_T_base_solver_type
+  !  end type mld_T_base_solver_type
+  !
+  !    build      -   Compute the actual contents of the smoother; includes
+  !                   invocation of the build method on the solver component. 
+  !    free       -   Release memory
+  !    apply      -   Apply the smoother to a vector (or to an array); includes
+  !                   invocation of the apply method on the solver component. 
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  !
+
+  type mld_z_base_solver_type
+  contains
+    procedure, pass(sv) :: check => z_base_solver_check
+    procedure, pass(sv) :: dump  => z_base_solver_dmp
+    procedure, pass(sv) :: build => z_base_solver_bld
+    procedure, pass(sv) :: apply_v => z_base_solver_apply_vect
+    procedure, pass(sv) :: apply_a => z_base_solver_apply
+    generic, public     :: apply => apply_a, apply_v
+    procedure, pass(sv) :: free  => z_base_solver_free
+    procedure, pass(sv) :: seti  => z_base_solver_seti
+    procedure, pass(sv) :: setc  => z_base_solver_setc
+    procedure, pass(sv) :: setr  => z_base_solver_setr
+    generic, public     :: set   => seti, setc, setr
+    procedure, pass(sv) :: default => z_base_solver_default
+    procedure, pass(sv) :: descr   => z_base_solver_descr
+    procedure, pass(sv) :: sizeof  => z_base_solver_sizeof
+    procedure, pass(sv) :: get_nzeros => z_base_solver_get_nzeros
+  end type mld_z_base_solver_type
+
+  private :: z_base_solver_bld,  z_base_solver_apply, &
+       &  z_base_solver_free,    z_base_solver_seti, &
+       &  z_base_solver_setc,    z_base_solver_setr, &
+       &  z_base_solver_descr,   z_base_solver_sizeof, &
+       &  z_base_solver_default, z_base_solver_check,&
+       &  z_base_solver_dmp, z_base_solver_apply_vect, &
+       &  z_base_solver_get_nzeros
+
+
+
+contains
+  !
+  ! Function returning the size of the data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function z_base_solver_sizeof(sv) result(val)
+    implicit none 
+    ! Arguments
+    class(mld_z_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_)                  :: val
+    integer             :: i
+    val = 0
+
+    return
+  end function z_base_solver_sizeof
+
+  function z_base_solver_get_nzeros(sv) result(val)
+    implicit none 
+    class(mld_z_base_solver_type), intent(in) :: sv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+  end function z_base_solver_get_nzeros
+
+  
+  !
+  ! Apply: comes in two versions, on plain arrays or on encapsulated
+  ! vectors.
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  ! Question: would it make sense to transform the base version into
+  ! the ID version, i.e. "base_solver" is the identity operator? 
+  ! 
+
+  subroutine z_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)           :: desc_data
+    class(mld_z_base_solver_type), intent(in) :: sv
+    complex(psb_dpk_),intent(inout)              :: x(:)
+    complex(psb_dpk_),intent(inout)              :: y(:)
+    complex(psb_dpk_),intent(in)                 :: alpha,beta
+    character(len=1),intent(in)               :: trans
+    complex(psb_dpk_),target, intent(inout)      :: work(:)
+    integer, intent(out)                      :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine z_base_solver_apply
+
+  subroutine z_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    use psb_base_mod
+    type(psb_desc_type), intent(in)              :: desc_data
+    class(mld_z_base_solver_type), intent(inout) :: sv
+    type(psb_z_vect_type),intent(inout)          :: x
+    type(psb_z_vect_type),intent(inout)          :: y
+    complex(psb_dpk_),intent(in)                    :: alpha,beta
+    character(len=1),intent(in)                  :: trans
+    complex(psb_dpk_),target, intent(inout)         :: work(:)
+    integer, intent(out)                         :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_apply'
+
+    call psb_erractionsave(err_act)
+    
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+    
+  end subroutine z_base_solver_apply_vect
+
+
+  !
+  ! Build
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. 
+  !
+  subroutine z_base_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    type(psb_zspmat_type), intent(in), target           :: a
+    Type(psb_desc_type), Intent(in)                     :: desc_a 
+    class(mld_z_base_solver_type), intent(inout)        :: sv
+    character, intent(in)                               :: upd
+    integer, intent(out)                                :: info
+    type(psb_zspmat_type), intent(in), target, optional :: b
+    class(psb_z_base_sparse_mat), intent(in), optional  :: amold
+    class(psb_z_base_vect_type), intent(in), optional   :: vmold
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_solver_bld'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_solver_bld
+
+  subroutine z_base_solver_check(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_solver_check
+
+  !
+  ! Set.
+  ! The base version does nothing; the principle is that
+  ! SET acts on what is known, and delegates what is unknown.
+  ! Since we are at the bottom of the hierarchy, there's no one
+  ! to delegate, so we do nothing. 
+  !
+  subroutine z_base_solver_seti(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_seti'
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine z_base_solver_seti
+
+  subroutine z_base_solver_setc(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv
+    integer, intent(in)                          :: what 
+    character(len=*), intent(in)                 :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act, ival 
+    character(len=20) :: name='d_base_solver_setc'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call sv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_solver_setc
+  
+  subroutine z_base_solver_setr(sv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv 
+    integer, intent(in)                          :: what 
+    real(psb_dpk_), intent(in)                   :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_setr'
+
+    
+    ! Correct action here is doing nothing. 
+    info = 0
+    
+    return
+  end subroutine z_base_solver_setr
+
+  !
+  ! Free
+  ! The base version throws an error, since it means that no explicit
+  ! choice was made. IS THIS CORRECT? I suspect it would be better
+  ! to be silent here, to cover reallocation. 
+  !
+  subroutine z_base_solver_free(sv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_solver_free'
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_solver_free
+
+  subroutine z_base_solver_descr(sv,info,iout,coarse)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_z_base_solver_type), intent(in) :: sv
+    integer, intent(out)                      :: info
+    integer, intent(in), optional             :: iout
+    logical, intent(in), optional             :: coarse
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_z_base_solver_descr'
+    integer      :: iout_
+
+
+    call psb_erractionsave(err_act)
+
+    info = psb_err_missing_override_method_
+    call psb_errpush(info,name)
+    goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_solver_descr
+
+  !
+  ! Dump. For debugging purposes. 
+  !
+  subroutine z_base_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_z_base_solver_type), intent(in) :: sv
+    integer, intent(in)              :: ictxt,level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: solver_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_slv_d"
+    end if
+
+    call psb_info(ictxt,iam,np)
+
+    if (present(solver)) then 
+      solver_ = solver
+    else
+      solver_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    ! At base level do nothing for the solver
+
+  end subroutine z_base_solver_dmp
+
+  subroutine z_base_solver_default(sv) 
+    implicit none 
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv
+    ! Do nothing for base version
+
+    return
+  end subroutine z_base_solver_default
+
+
+
+end module mld_z_base_solver_mod

mlprec/mld_z_diag_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_z_diag_solver
 
-  use mld_z_prec_type
+  use mld_z_base_solver_mod
 
   type, extends(mld_z_base_solver_type) :: mld_z_diag_solver_type
     type(psb_z_vect_type), allocatable :: dv

mlprec/mld_z_id_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_z_id_solver
 
-  use mld_z_prec_type
+  use mld_z_base_solver_mod
 
   type, extends(mld_z_base_solver_type) :: mld_z_id_solver_type
   contains

mlprec/mld_z_ilu_solver.f90

@@ -45,7 +45,7 @@
 
 module mld_z_ilu_solver
 
-  use mld_z_prec_type
+  use mld_z_base_solver_mod
   use mld_z_ilu_fact_mod
 
   type, extends(mld_z_base_solver_type) :: mld_z_ilu_solver_type

mlprec/mld_z_jac_smoother.f90

@@ -44,7 +44,7 @@
 !
 module mld_z_jac_smoother
 
-  use mld_z_prec_type
+  use mld_z_base_smoother_mod
 
   type, extends(mld_z_base_smoother_type) :: mld_z_jac_smoother_type
     ! The local solver component is inherited from the

mlprec/mld_z_onelev_mod.f90

@@ -0,0 +1,666 @@
+!!$
+!!$ 
+!!$                           MLD2P4  version 2.0
+!!$  MultiLevel Domain Decomposition Parallel Preconditioners Package
+!!$             based on PSBLAS (Parallel Sparse BLAS version 3.0)
+!!$  
+!!$  (C) Copyright 2008,2009,2010
+!!$
+!!$                      Salvatore Filippone  University of Rome Tor Vergata
+!!$                      Alfredo Buttari      CNRS-IRIT, Toulouse
+!!$                      Pasqua D'Ambra       ICAR-CNR, Naples
+!!$                      Daniela di Serafino  Second University of Naples
+!!$ 
+!!$  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 MLD2P4 group or the names of its contributors may
+!!$       not be used to endorse or promote products derived from this
+!!$       software without specific 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 MLD2P4 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: mld_z_onelev_mod.f90
+!
+! Module: mld_z_onelev_mod
+!
+!  This module defines: 
+!  - the mld_z_onelev_type data structure containing one level
+!    of a multilevel  preconditioner and related
+!    data structures;
+!
+!  It contains routines for
+!  - Building and applying; 
+!  - checking if the preconditioner is correctly defined;
+!  - printing a	description of the preconditioner;
+!  - deallocating the preconditioner data structure.  
+!
+
+module mld_z_onelev_mod
+
+  use mld_base_prec_type
+  use psb_base_mod, only : psb_z_vect_type, psb_z_base_vect_type
+  use mld_z_base_smoother_mod
+  !
+  !
+  ! Type: mld_Tonelev_type.
+  !
+  !  It is the data type containing the necessary items for the	current
+  !  level (essentially, the smoother, the current-level matrix
+  !  and the restriction and prolongation operators).
+  !
+  !  type mld_Tonelev_type
+  !    class(mld_T_base_smoother_type), allocatable :: sm
+  !    type(mld_RTml_parms)            :: parms 
+  !    type(psb_Tspmat_type)           :: ac
+  !    type(psb_Tesc_type)             :: desc_ac
+  !    type(psb_Tspmat_type), pointer  :: base_a    => null() 
+  !    type(psb_Tesc_type), pointer    :: base_desc => null() 
+  !    type(psb_Tlinmap_type)          :: map
+  !  end type mld_Tonelev_type
+  !
+  !  Note that psb_Tpk denotes the kind of the real data type to be chosen
+  !  according to single/double precision version of MLD2P4.
+  !
+  !   sm           -  class(mld_T_base_smoother_type), allocatable
+  !                   The current level preconditioner (aka smoother).
+  !   parms        -  type(mld_RTml_parms)
+  !                   The parameters defining the multilevel strategy.
+  !   ac           -  The local part of the current-level matrix, built by
+  !                   coarsening the previous-level matrix.
+  !   desc_ac      -  type(psb_desc_type).
+  !                   The communication descriptor associated to the matrix
+  !                   stored in ac.
+  !   base_a       -  type(psb_Tspmat_type), pointer.
+  !                   Pointer (really a pointer!) to the local part of the current 
+  !                   matrix (so we have a unified treatment of residuals).
+  !                   We need this to avoid passing explicitly the current matrix
+  !                   to the routine which applies the preconditioner.
+  !   base_desc    -  type(psb_desc_type), pointer.
+  !                   Pointer to the communication descriptor associated to the
+  !                   matrix pointed by base_a.
+  !   map          -  Stores the maps (restriction and prolongation) between the
+  !                   vector spaces associated to the index spaces of the previous
+  !                   and current levels.
+  !
+  !   Methods:  
+  !     Most methods follow the encapsulation hierarchy: they take whatever action
+  !     is appropriate for the current object, then call the corresponding method for
+  !     the contained object.
+  !     As an example: the descr() method prints out a description of the
+  !     level. It starts by invoking the descr() method of the parms object,
+  !     then calls the descr() method of the smoother object. 
+  !
+  !    descr      -   Prints a description of the object.
+  !    default    -   Set default values
+  !    dump       -   Dump to file object contents
+  !    set        -   Sets various parameters; when a request is unknown
+  !                   it is passed to the smoother object for further processing.
+  !    check      -   Sanity checks.
+  !    sizeof     -   Total memory occupation in bytes
+  !    get_nzeros -   Number of nonzeros 
+  !
+  !
+  type mld_zonelev_type
+    class(mld_z_base_smoother_type), allocatable :: sm
+    type(mld_dml_parms)             :: parms 
+    type(psb_zspmat_type)           :: ac
+    type(psb_desc_type)             :: desc_ac
+    type(psb_zspmat_type), pointer  :: base_a    => null() 
+    type(psb_desc_type), pointer    :: base_desc => null() 
+    type(psb_zlinmap_type)          :: map
+  contains
+    procedure, pass(lv) :: descr   => z_base_onelev_descr
+    procedure, pass(lv) :: default => z_base_onelev_default
+    procedure, pass(lv) :: free    => z_base_onelev_free
+    procedure, pass(lv) :: nullify => z_base_onelev_nullify
+    procedure, pass(lv) :: check => z_base_onelev_check
+    procedure, pass(lv) :: dump  => z_base_onelev_dump
+    procedure, pass(lv) :: seti  => z_base_onelev_seti
+    procedure, pass(lv) :: setr  => z_base_onelev_setr
+    procedure, pass(lv) :: setc  => z_base_onelev_setc
+    generic, public     :: set   => seti, setr, setc
+    procedure, pass(lv) :: sizeof => z_base_onelev_sizeof
+    procedure, pass(lv) :: get_nzeros => z_base_onelev_get_nzeros
+  end type mld_zonelev_type
+
+  private :: z_base_onelev_seti, z_base_onelev_setc, &
+       &  z_base_onelev_setr, z_base_onelev_check, &
+       &  z_base_onelev_default, z_base_onelev_dump, &
+       &  z_base_onelev_descr,  z_base_onelev_sizeof, &
+       &  z_base_onelev_free,  z_base_onelev_nullify,&
+       &  z_base_onelev_get_nzeros
+
+
+  interface mld_nullify_onelevprec
+    module procedure  mld_nullify_d_onelevprec
+  end interface
+
+
+contains
+  !
+  ! Function returning the size of the mld_prec_type data structure
+  ! in bytes or in number of nonzeros of the operator(s) involved. 
+  !
+
+  function z_base_onelev_get_nzeros(lv) result(val)
+    implicit none 
+    class(mld_zonelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    val = 0
+    if (allocated(lv%sm)) &
+         &  val =  lv%sm%get_nzeros()
+  end function z_base_onelev_get_nzeros
+
+  function z_base_onelev_sizeof(lv) result(val)
+    implicit none 
+    class(mld_zonelev_type), intent(in) :: lv
+    integer(psb_long_int_k_) :: val
+    integer             :: i
+    
+    val = 0
+    val = val + lv%desc_ac%sizeof()
+    val = val + lv%ac%sizeof()
+    val = val + lv%map%sizeof() 
+    if (allocated(lv%sm))  val = val + lv%sm%sizeof()
+  end function z_base_onelev_sizeof
+
+
+  !
+  ! Subroutine: mld_file_onelev_descr
+  ! Version: complex
+  !
+  !  This routine prints a description of the preconditioner to the standard 
+  !  output or to a file. It must be called after the preconditioner has been
+  !  built by mld_precbld.
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tprec_type), input.
+  !             The preconditioner data structure to be printed out.
+  !  info    -  integer, output.
+  !             error code.
+  !  iout    -  integer, input, optional.
+  !             The id of the file where the preconditioner description
+  !             will be printed. If iout is not present, then the standard
+  !             output is condidered.
+  !
+  subroutine z_base_onelev_descr(lv,il,nl,info,iout)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_zonelev_type), intent(in) :: lv
+    integer, intent(in)                 :: il,nl
+    integer, intent(out)                :: info
+    integer, intent(in), optional       :: iout
+
+    ! Local variables
+    integer      :: err_act
+    integer      :: ictxt, me, np
+    character(len=20), parameter :: name='mld_z_base_onelev_descr'
+    integer      :: iout_
+    logical      :: coarse
+
+
+    call psb_erractionsave(err_act)
+
+
+    coarse = (il==nl)
+
+    if (present(iout)) then 
+      iout_ = iout
+    else 
+      iout_ = 6
+    end if
+    
+    write(iout_,*) 
+    if (il == 2) then 
+      call lv%parms%mldescr(iout_,info)
+      write(iout_,*) 
+    end if
+
+    if (coarse)  then 
+      write(iout_,*) ' Level ',il,' (coarsest)'
+    else
+      write(iout_,*) ' Level ',il
+    end if
+
+    call lv%parms%descr(iout_,info,coarse=coarse)
+        
+    if (nl > 1) then 
+      if (allocated(lv%map%naggr)) then
+        write(iout_,*) '  Size of coarse matrix: ', &
+             &  sum(lv%map%naggr(:))
+        write(iout_,*) '  Sizes of aggregates: ', &
+             &  lv%map%naggr(:)
+      end if
+    end if
+
+    if (coarse.and.allocated(lv%sm)) &
+         & call lv%sm%descr(info,iout=iout_,coarse=coarse)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_onelev_descr
+
+
+  !
+  ! Subroutines: mld_T_onelev_precfree
+  ! Version: complex
+  !
+  !  These routines deallocate the mld_Tonelev_type 
+  !
+  ! Arguments:
+  !  p       -  type(mld_Tonelev_type), input.
+  !             The data structure to be deallocated.
+  !  info    -  integer, output.
+  !             error code.
+  !
+  subroutine z_base_onelev_free(lv,info)
+    use psb_base_mod
+    implicit none 
+
+    class(mld_zonelev_type), intent(inout) :: lv
+    integer, intent(out)                :: info
+    integer :: i
+
+    info = psb_success_
+
+    ! We might just deallocate the top level array, except 
+    ! that there may be inner objects containing C pointers,
+    ! e.g.  UMFPACK, SLU or CUDA stuff.
+    ! We really need FINALs. 
+    call lv%sm%free(info)
+    
+    call lv%ac%free()
+    if (psb_is_ok_desc(lv%desc_ac)) &
+         & call psb_cdfree(lv%desc_ac,info)
+    call lv%map%free(info)
+
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_a) 
+    ! This is a pointer to something else, must not free it here. 
+    nullify(lv%base_desc) 
+
+    call lv%nullify()
+
+  end subroutine z_base_onelev_free
+
+
+  subroutine z_base_onelev_nullify(lv)
+    implicit none 
+
+    class(mld_zonelev_type), intent(inout) :: lv
+
+    nullify(lv%base_a) 
+    nullify(lv%base_desc) 
+
+  end subroutine z_base_onelev_nullify
+
+
+  subroutine mld_nullify_d_onelevprec(p)
+    implicit none 
+
+    type(mld_zonelev_type), intent(inout) :: p
+
+    nullify(p%base_a) 
+    nullify(p%base_desc) 
+
+  end subroutine mld_nullify_d_onelevprec
+
+  !
+  ! Onelevel checks. 
+  ! The number of Jacobi sweeps to be applied is not
+  ! tied to the Jacobi smoother: logically, you have
+  ! a smoother and you can choose to apply it any number
+  ! of times you like. 
+  !
+  subroutine z_base_onelev_check(lv,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_zonelev_type), intent(inout) :: lv 
+    integer, intent(out)                   :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_check'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    call mld_check_def(lv%parms%sweeps,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_pre,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+    call mld_check_def(lv%parms%sweeps_post,&
+         & 'Jacobi sweeps',1,is_legal_jac_sweeps)
+
+    
+    if (allocated(lv%sm)) then 
+      call lv%sm%check(info)
+    else 
+      info=3111
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    if (info /= psb_success_) goto 9999
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_onelev_check
+
+  !
+  ! Multilevel defaults: 
+  !  multiplicative vs. additive ML framework;
+  !  Smoothed decoupled aggregation with zero threshold; 
+  !  distributed coarse matrix;
+  !  damping omega  computed with the max-norm estimate of the
+  !  dominant eigenvalue;
+  !  two-sided smoothing (i.e. V-cycle) with 1 smoothing sweep;
+  !
+
+  subroutine z_base_onelev_default(lv)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_zonelev_type), intent(inout) :: lv 
+
+    lv%parms%sweeps          = 1
+    lv%parms%sweeps_pre      = 1
+    lv%parms%sweeps_post     = 1
+    lv%parms%ml_type         = mld_mult_ml_
+    lv%parms%aggr_alg        = mld_dec_aggr_
+    lv%parms%aggr_kind       = mld_smooth_prol_
+    lv%parms%coarse_mat      = mld_distr_mat_
+    lv%parms%smoother_pos    = mld_twoside_smooth_
+    lv%parms%aggr_omega_alg  = mld_eig_est_
+    lv%parms%aggr_eig        = mld_max_norm_
+    lv%parms%aggr_filter     = mld_no_filter_mat_
+    lv%parms%aggr_omega_val  = dzero
+    lv%parms%aggr_thresh     = dzero
+    
+    if (allocated(lv%sm)) call lv%sm%default()
+
+    return
+
+  end subroutine z_base_onelev_default
+
+  ! 
+  ! Set routines:
+  ! Parameters belonging here are:
+  !   Number of smoothing sweeps;
+  !   Smoother position;
+  !   Aggregation related parameters
+  !   Record request on coarse level solver,
+  !    for checks on solver vs. smoother nomenclature
+  !    reconciliation. 
+  ! 
+  subroutine z_base_onelev_seti(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_zonelev_type), intent(inout) :: lv 
+    integer, intent(in)                          :: what 
+    integer, intent(in)                          :: val
+    integer, intent(out)                         :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_seti'
+
+    call psb_erractionsave(err_act)
+    info = psb_success_
+
+    select case (what) 
+
+    case (mld_smoother_sweeps_)
+      lv%parms%sweeps      = val
+      lv%parms%sweeps_pre  = val
+      lv%parms%sweeps_post = val
+
+    case (mld_smoother_sweeps_pre_)
+      lv%parms%sweeps_pre  = val
+
+    case (mld_smoother_sweeps_post_)
+      lv%parms%sweeps_post = val
+
+    case (mld_ml_type_)
+      lv%parms%ml_type       = val
+
+    case (mld_aggr_alg_)
+      lv%parms%aggr_alg      = val
+
+    case (mld_aggr_kind_)
+      lv%parms%aggr_kind     = val
+
+    case (mld_coarse_mat_)
+      lv%parms%coarse_mat    = val
+
+    case (mld_smoother_pos_)
+      lv%parms%smoother_pos  = val
+
+    case (mld_aggr_omega_alg_)
+      lv%parms%aggr_omega_alg= val
+
+    case (mld_aggr_eig_)
+      lv%parms%aggr_eig      = val
+
+    case (mld_aggr_filter_)
+      lv%parms%aggr_filter   = val
+
+    case (mld_coarse_solve_)
+      lv%parms%coarse_solve    = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_onelev_seti
+
+  subroutine z_base_onelev_setc(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_zonelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    character(len=*), intent(in)                   :: val
+    integer, intent(out)                           :: info
+    Integer           :: err_act
+    character(len=20) :: name='d_base_onelev_setc'
+    integer :: ival 
+
+    call psb_erractionsave(err_act)
+
+    info = psb_success_
+
+    call mld_stringval(val,ival,info)
+    if (info == psb_success_) call lv%set(what,ival,info)
+
+    if (info /= psb_success_) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_onelev_setc
+  
+  subroutine z_base_onelev_setr(lv,what,val,info)
+
+    use psb_base_mod
+
+    Implicit None
+
+    ! Arguments
+    class(mld_zonelev_type), intent(inout) :: lv 
+    integer, intent(in)                            :: what 
+    real(psb_dpk_), intent(in)                     :: val
+    integer, intent(out)                           :: info
+    Integer :: err_act
+    character(len=20)  :: name='d_base_onelev_setr'
+
+    call psb_erractionsave(err_act)
+
+
+    info = psb_success_
+    
+    select case (what) 
+
+    case (mld_aggr_omega_val_)
+      lv%parms%aggr_omega_val= val
+   
+    case (mld_aggr_thresh_)
+      lv%parms%aggr_thresh   = val
+
+    case default
+      if (allocated(lv%sm)) then 
+        call lv%sm%set(what,val,info)
+      end if
+      if (info /= psb_success_) goto 9999
+    end select
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine z_base_onelev_setr
+
+  !
+  ! Dump on file: can be fine-tuned to include the (aggregated) matrix
+  ! as well as smoother and solver. 
+  !
+  subroutine z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,smoother,solver)
+    use psb_base_mod
+    implicit none 
+    class(mld_zonelev_type), intent(in) :: lv
+    integer, intent(in)              :: level
+    integer, intent(out)             :: info
+    character(len=*), intent(in), optional :: prefix, head
+    logical, optional, intent(in)    :: ac, rp, smoother, solver
+    integer :: i, j, il1, iln, lname, lev
+    integer :: icontxt,iam, np
+    character(len=80)  :: prefix_
+    character(len=120) :: fname ! len should be at least 20 more than
+    logical :: ac_, rp_
+    !  len of prefix_ 
+
+    info = 0
+
+    if (present(prefix)) then 
+      prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
+    else
+      prefix_ = "dump_lev_d"
+    end if
+
+    if (associated(lv%base_desc)) then 
+      icontxt = lv%base_desc%get_context()
+      call psb_info(icontxt,iam,np)
+    else 
+      icontxt = -1 
+      iam     = -1
+    end if
+    if (present(ac)) then 
+      ac_ = ac
+    else
+      ac_ = .false. 
+    end if
+    if (present(rp)) then 
+      rp_ = rp
+    else
+      rp_ = .false. 
+    end if
+    lname = len_trim(prefix_)
+    fname = trim(prefix_)
+    write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
+    lname = lname + 5
+
+    if (level >= 2) then 
+      if (ac_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_ac.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%ac%print(fname,head=head)
+      end if
+      if (rp_) then 
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_r.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_X2Y%print(fname,head=head)
+        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_p.mtx'
+        write(0,*) 'Filename ',fname
+        call lv%map%map_Y2X%print(fname,head=head)
+      end if
+    end if
+    if (allocated(lv%sm)) &
+         & call lv%sm%dump(icontxt,level,info,smoother=smoother,solver=solver)
+
+  end subroutine z_base_onelev_dump
+
+
+end module mld_z_onelev_mod

mlprec/mld_z_slu_solver.f90

@@ -46,7 +46,7 @@
 module mld_z_slu_solver
 
   use iso_c_binding
-  use mld_z_prec_type
+  use mld_z_base_solver_mod
 
   type, extends(mld_z_base_solver_type) :: mld_z_slu_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_z_sludist_solver.f90

@@ -46,7 +46,7 @@
 module mld_z_sludist_solver
 
   use iso_c_binding
-  use mld_z_prec_type
+  use mld_z_base_solver_mod
 
   type, extends(mld_z_base_solver_type) :: mld_z_sludist_solver_type
     type(c_ptr)                 :: lufactors=c_null_ptr

mlprec/mld_z_umf_solver.f90

@@ -46,7 +46,7 @@
 module mld_z_umf_solver
 
   use iso_c_binding
-  use mld_z_prec_type
+  use mld_z_base_solver_mod
 
   type, extends(mld_z_base_solver_type) :: mld_z_umf_solver_type
     type(c_ptr)                 :: symbolic=c_null_ptr, numeric=c_null_ptr

tests/fileread/data_input.f90

@@ -40,9 +40,10 @@ module data_input
   
   interface read_data
     module procedure read_char, read_int,&
-         & read_double, read_single,&
+         & read_double, read_single, read_logical,&
          & string_read_char, string_read_int,&
-         & string_read_double, string_read_single
+         & string_read_double, string_read_single, &
+         & string_read_logical
   end interface read_data
   interface trim_string
     module procedure trim_string
@@ -53,6 +54,16 @@ module data_input
 
 contains
 
+  subroutine read_logical(val,file,marker)
+    logical, intent(out) :: val
+    integer, intent(in)           :: file
+    character(len=1), optional, intent(in) :: marker
+
+    read(file,'(a)')charbuf
+    call read_data(val,charbuf,marker)
+
+  end subroutine read_logical
+
   subroutine read_char(val,file,marker)
     character(len=*), intent(out) :: val
     integer, intent(in)           :: file
@@ -63,6 +74,7 @@ contains
 
   end subroutine read_char
 
+
   subroutine read_int(val,file,marker)
     integer, intent(out) :: val
     integer, intent(in)  :: file
@@ -130,6 +142,7 @@ contains
     if (idx == 0) idx = len(charbuf)+1
     read(charbuf(1:idx-1),*) val
   end subroutine string_read_int
+
   subroutine string_read_single(val,file,marker)
     use psb_base_mod
     real(psb_spk_), intent(out) :: val
@@ -149,6 +162,7 @@ contains
     if (idx == 0) idx = len(charbuf)+1
     read(charbuf(1:idx-1),*) val
   end subroutine string_read_single
+
   subroutine string_read_double(val,file,marker)
     use psb_base_mod
     real(psb_dpk_), intent(out) :: val
@@ -169,6 +183,26 @@ contains
     read(charbuf(1:idx-1),*) val
   end subroutine string_read_double
 
+  subroutine string_read_logical(val,file,marker)
+    use psb_base_mod
+    logical, intent(out) :: val
+    character(len=*), intent(in)         :: file
+    character(len=1), optional, intent(in) :: marker
+    character(len=1)    :: marker_
+    character(len=1024) :: charbuf
+    integer :: idx
+    if (present(marker)) then 
+      marker_ = marker
+    else
+      marker_ = def_marker
+    end if
+    read(file,'(a)')charbuf
+    charbuf = adjustl(charbuf)
+    idx=index(charbuf,marker_)
+    if (idx == 0) idx = len(charbuf)+1
+    read(charbuf(1:idx-1),*) val
+  end subroutine string_read_logical
+
   function  trim_string(string,marker)
     character(len=*), intent(in) :: string
     character(len=1), optional, intent(in) :: marker

tests/fileread/df_sample.f90

@@ -71,6 +71,7 @@ program df_sample
     real(psb_dpk_)     :: cthres      ! threshold for coarse fact.  ILU(T)
     integer            :: cjswp       ! block-Jacobi sweeps
     real(psb_dpk_)     :: athres      ! smoothed aggregation threshold
+    logical            :: dump        ! Dump preconditioner on file
   end type precdata
   type(precdata)        :: prec_choice
 
@@ -314,7 +315,10 @@ program df_sample
   call psb_spmm(-done,a,x_col,done,r_col,desc_a,info)
   resmx  = psb_genrm2(r_col,desc_a,info)
   resmxp = psb_geamax(r_col,desc_a,info)
-
+  if (prec_choice%dump) &
+       & call prec%dump(info,istart=1,prefix="out-"//trim(prec_choice%solve),&
+       &  solver=.true.)
+  
   amatsize = a%sizeof()
   descsize = desc_a%sizeof()
   precsize = mld_sizeof(prec)
@@ -409,6 +413,7 @@ contains
       call read_data(itrace,5)
       call read_data(irst,5)
       call read_data(eps,5)
+      call read_data(prec%dump,5)        ! dump prec on file
       call read_data(prec%descr,5)       ! verbose description of the prec
       call read_data(prec%prec,5)        ! overall prectype
       call read_data(prec%novr,5)        ! number of overlap layers
@@ -447,6 +452,7 @@ contains
     call psb_bcast(icontxt,irst)
     call psb_bcast(icontxt,eps)
 
+    call psb_bcast(icontxt,prec%dump)        ! 
     call psb_bcast(icontxt,prec%descr)       ! verbose description of the prec
     call psb_bcast(icontxt,prec%prec)        ! overall prectype
     call psb_bcast(icontxt,prec%novr)        ! number of overlap layers

tests/fileread/runs/dfs.inp

@@ -31,5 +31,5 @@ UMF                         ! Coarsest-level subsolver: ILU UMF SLU SLUDIST (DSC
 1.d-4                       ! Coarsest-level threshold T for ILU(T,P) 
 4                           ! Number of Jacobi sweeps for BJAC/PJAC coarsest-level solver
 0.01d0                      ! Smoothed aggregation threshold: >= 0.0
-F                           ! Read aggregation maps from file 
+F                           ! dump preconditioner data.
-mlaggr.out-FRT.lev          ! File name prefix 
+