amg4psblas/psb_prec_type.f90

!!$ 
!!$ 
!!$                    MD2P4
!!$    Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS
!!$                      for 
!!$              Parallel Sparse BLAS  v2.0
!!$    (C) Copyright 2006 Salvatore Filippone    University of Rome Tor Vergata
!!$                       Alfredo Buttari        University of Rome Tor Vergata
!!$                       Daniela di Serafino    Second University of Naples
!!$                       Pasqua D'Ambra         ICAR-CNR                      
!!$ 
!!$  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 MD2P4 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 MD2P4 GROUP OR ITS CONTRIBUTORS
!!$  BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$  POSSIBILITY OF SUCH DAMAGE.
!!$ 
module psb_prec_type
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!	Module to   define PREC_DATA,           !!
!!      structure for preconditioning.          !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! Reduces size of .mod file. Without the ONLY clause compilation 
  ! blows up on some systems.
  use psb_base_mod, only : psb_dspmat_type, psb_zspmat_type, psb_desc_type,&
       & psb_sizeof

  integer, parameter :: min_prec_=0, noprec_=0, diag_=1, bjac_=2,&
       & ras_=3,asm_=4, ash_=5, rash_=6, ras2lv_=7, ras2lvm_=8,&
       & lv2mras_=9, lv2smth_=10, lv2lsm_=11, sl2sm_=12, superlu_=13,&
       & new_loc_smth_=14, new_glb_smth_=15, ag2lsm_=16,&
       & msy2l_=18, msy2g_=19, max_prec_=19

  ! Multilevel stuff.
  integer, parameter :: no_ml_=0, add_ml_prec_=1, mult_ml_prec_=2
  integer, parameter :: new_ml_prec_=3, max_ml_=new_ml_prec_
  integer, parameter :: pre_smooth_=1, post_smooth_=2, smooth_both_=3,&
       &    max_smooth_=smooth_both_
  integer, parameter :: loc_aggr_=0, glb_aggr_=1, new_loc_aggr_=2
  integer, parameter :: new_glb_aggr_=3, max_aggr_=new_glb_aggr_
  integer, parameter :: no_smth_=0, smth_omg_=1, smth_biz_=2
  integer, parameter :: lib_choice_=0, user_choice_=1
  integer, parameter :: mat_distr_=0, mat_repl_=1
  ! Entries in iprcparm: preconditioner type, factorization type,
  ! prolongation type, restriction type, renumbering algorithm,
  ! number of overlap layers, pointer to SuperLU factors, 
  ! levels of fill in for ILU(N), 
  integer, parameter :: p_type_=1, f_type_=2, restr_=3, prol_=4
  integer, parameter :: iren_=5, n_ovr_=6
  integer, parameter :: ilu_fill_in_=8, jac_sweeps_=9, ml_type_=10
  integer, parameter :: smth_pos_=11, aggr_alg_=12, smth_kind_=13
  integer, parameter :: om_choice_=14, glb_smth_=15, coarse_mat_=16
  !! 2 ints for 64 bit versions
  integer, parameter :: slu_ptr_=17, umf_symptr_=17, umf_numptr_=19
  integer, parameter :: slud_ptr_=21
  integer, parameter :: ifpsz=24, prec_status_=ifpsz
  integer, parameter :: prec_built=98765
  !Renumbering. SEE BELOW
  integer, parameter :: renum_none_=0, renum_glb_=1, renum_gps_=2
  ! Entries in dprcparm: ILU(E) epsilon, smoother omega
  integer, parameter :: fact_eps_=1, smooth_omega_=2
  integer, parameter :: dfpsz=4
  ! Factorization types: none, ILU(N), ILU(E), SuperLU, UMFPACK
  integer, parameter :: f_none_=0,f_ilu_n_=1,f_ilu_e_=2,f_slu_=3
  integer, parameter :: f_umf_=4, f_slud_=5
  ! Fields for sparse matrices ensembles: 
  integer, parameter :: l_pr_=1, u_pr_=2, bp_ilu_avsz=2
  integer, parameter :: ap_nd_=3, ac_=4, sm_pr_t_=5, sm_pr_=6
  integer, parameter :: smth_avsz=6, max_avsz=smth_avsz 


  type psb_dbaseprc_type

    type(psb_dspmat_type), allocatable :: av(:) 
    real(kind(1.d0)), allocatable      :: d(:)  
    type(psb_desc_type)                :: desc_data , desc_ac
    integer, allocatable               :: iprcparm(:) 
    real(kind(1.d0)), allocatable      :: dprcparm(:) 
    integer, allocatable               :: perm(:),  invperm(:) 
    integer, allocatable               :: mlia(:), nlaggr(:) 
    type(psb_dspmat_type), pointer     :: base_a    => null() !
    type(psb_desc_type), pointer       :: base_desc=> null() ! 
    real(kind(1.d0)), allocatable      :: dorig(:) 

  end type psb_dbaseprc_type


  !
  !  Multilevel preconditioning
  !
  !  To each level I there corresponds a matrix A(I) and a preconditioner K(I)
  !
  !  A notational difference: in the DD reference above the preconditioner for 
  !  a given level K(I) is written out as a sum over the subdomains
  !
  !  SUM_k(R_k^T A_k R_k) 
  !
  !  whereas in this code the sum is implicit in the parallelization, 
  !  i.e. each process takes care of one subdomain, and for each level we have 
  !  as many subdomains as there are processes (except for the coarsest level where 
  !  we might have a replicated index space). Thus the sum apparently disappears 
  !  from our code, but only apparently, because it is implicit in the call 
  !  to psb_baseprc_aply. 
  !
  !  A bit of description of the baseprecv(:) data structure:
  !   1. Number of levels = NLEV = size(baseprecv(:))
  !   2. baseprecv(ilev)%av(:)    sparse matrices needed for the current level. 
  !      Includes:
  !   2.1.:  baseprecv(ilev)%av(l_pr_)    L factor of ILU preconditioners
  !   2.2.:  baseprecv(ilev)%av(u_pr_)    U factor of ILU preconditioners
  !   2.3.:  baseprecv(ilev)%av(ap_nd_)   Off-diagonal part of A for Jacobi sweeps
  !   2.4.:  baseprecv(ilev)%av(ac_)      Aggregated matrix of level ILEV 
  !   2.5.:  baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors  
  !                                          (ilev-1) --->  (ilev) 
  !   2.6.:  baseprecv(ilev)%av(sm_pr_)   Smoother prolongator; maps vectors  
  !                                          (ilev)   --->  (ilev-1) 
  !   Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe 
  !
  !   3.    baseprecv(ilev)%desc_data     comm descriptor for level ILEV
  !   4.    baseprecv(ilev)%base_a        Pointer (really a pointer!) to the base matrix 
  !                                       of the current level, i.e.: if ILEV=1 then  A
  !                                       else the aggregated matrix av(ac_); so we have 
  !                                       a unified treatment of residuals. Need this to 
  !                                       avoid passing explicitly matrix A to the 
  !                                       outer prec. routine
  !   5.    baseprecv(ilev)%mlia          The aggregation map from (ilev-1)-->(ilev)
  !                                       if no smoother, it is used instead of sm_pr_
  !   6.    baseprecv(ilev)%nlaggr        Number of aggregates on the various procs. 
  !   
  type psb_dprec_type
    type(psb_dbaseprc_type), allocatable  :: baseprecv(:) 
    ! contain type of preconditioning to be performed
    integer                       :: prec, base_prec
  end type psb_dprec_type

  type psb_zbaseprc_type

    type(psb_zspmat_type), allocatable :: av(:) 
    complex(kind(1.d0)), allocatable   :: d(:)  
    type(psb_desc_type)                :: desc_data , desc_ac
    integer, allocatable               :: iprcparm(:) 
    real(kind(1.d0)), allocatable      :: dprcparm(:) 
    integer, allocatable               :: perm(:),  invperm(:) 
    integer, allocatable               :: mlia(:), nlaggr(:) 
    type(psb_zspmat_type), pointer     :: base_a    => null() !
    type(psb_desc_type), pointer       :: base_desc => null() !
    complex(kind(1.d0)), allocatable   :: dorig(:)

  end type psb_zbaseprc_type

  type psb_zprec_type
    type(psb_zbaseprc_type), allocatable  :: baseprecv(:) 
    ! contain type of preconditioning to be performed
    integer                       :: prec, base_prec
  end type psb_zprec_type


  character(len=15), parameter, private :: &
       &  smooth_names(1:3)=(/'Pre-smoothing ','Post-smoothing',&
       & 'Smooth both   '/)
  character(len=15), parameter, private :: &
       &  smooth_kinds(0:2)=(/'No  smoother  ','Omega smoother',&
       &           'Bizr. smoother'/)
  character(len=15), parameter, private :: &
       &  matrix_names(0:1)=(/'Distributed   ','Replicated    '/)
  character(len=18), parameter, private :: &
       &  aggr_names(0:3)=(/'Local aggregation ','Global aggregation',&
       &     'New local aggr.   ','New global aggr.  '/)
  character(len=6), parameter, private :: &
       &  restrict_names(0:4)=(/'None ','     ','     ','     ','Halo '/)
  character(len=12), parameter, private :: &
       &  prolong_names(0:3)=(/'None       ','Sum        ','Average    ','Square root'/)
  character(len=15), parameter, private :: &
       &  ml_names(0:3)=(/'None          ','Additive      ','Multiplicative',&
       & 'New ML        '/)
  character(len=15), parameter, private :: &
       &  fact_names(0:5)=(/'None          ','ILU(n)        ',&
       &  'ILU(eps)      ','Sparse SuperLU','UMFPACK Sp. LU',&
       &  'SuperLU_Dist  '/)

  interface psb_base_precfree
    module procedure psb_dbase_precfree, psb_zbase_precfree
  end interface

  interface psb_nullify_baseprec
    module procedure psb_nullify_dbaseprec, psb_nullify_zbaseprec
  end interface

  interface psb_check_def
    module procedure psb_icheck_def, psb_dcheck_def
  end interface

  interface psb_prec_descr
    module procedure psb_out_prec_descr, psb_file_prec_descr, &
         &  psb_zout_prec_descr, psb_zfile_prec_descr
  end interface

  interface psb_prec_short_descr
    module procedure psb_prec_short_descr, psb_zprec_short_descr
  end interface

  interface psb_sizeof
    module procedure psb_dprec_sizeof, psb_zprec_sizeof, &
         & psb_dbaseprc_sizeof, psb_zbaseprc_sizeof
  end interface

contains

  function psb_dprec_sizeof(prec)
    use psb_base_mod
    type(psb_dprec_type), intent(in) :: prec
    integer             :: psb_dprec_sizeof
    integer             :: val,i
    val = 8
    if (allocated(prec%baseprecv)) then 
      do i=1, size(prec%baseprecv)
        val = val + psb_prec_sizeof(prec%baseprecv(i))
      end do
    end if
    psb_dprec_sizeof = val
  end function psb_dprec_sizeof

  function psb_zprec_sizeof(prec)
    use psb_base_mod
    type(psb_zprec_type), intent(in) :: prec
    integer             :: psb_zprec_sizeof
    integer             :: val,i
    val = 8
    if (allocated(prec%baseprecv)) then 
      do i=1, size(prec%baseprecv)
        val = val + psb_prec_sizeof(prec%baseprecv(i))
      end do
    end if
    psb_zprec_sizeof = val
  end function psb_zprec_sizeof

  function psb_dbaseprc_sizeof(prec)
    use psb_base_mod
    type(psb_dbaseprc_type), intent(in) :: prec
    integer             :: psb_dbaseprc_sizeof
    integer             :: val,i
    
    val = 0
    if (allocated(prec%iprcparm)) then 
      val = val + 4 * size(prec%iprcparm)
      if (prec%iprcparm(prec_status_) == prec_built) then 
        select case(prec%iprcparm(f_type_)) 
        case(f_ilu_n_,f_ilu_e_)
          ! do nothing
        case(f_slu_)
          write(0,*) 'Should implement check for size of SuperLU data structs'
        case(f_umf_)
          write(0,*) 'Should implement check for size of UMFPACK data structs'
        case(f_slud_)
          write(0,*) 'Should implement check for size of SuperLUDist data structs'
        case default
        end select
        
      end if
    end if
    if (allocated(prec%dprcparm)) val = val + 8 * size(prec%dprcparm)
    if (allocated(prec%d))        val = val + 8 * size(prec%d)
    if (allocated(prec%perm))     val = val + 4 * size(prec%perm)
    if (allocated(prec%invperm))  val = val + 4 * size(prec%invperm)
                                  val = val + psb_sizeof(prec%desc_data)
    if (allocated(prec%av))  then 
      do i=1,size(prec%av)
        val = val + psb_sizeof(prec%av(i))
      end do
    end if

    psb_dbaseprc_sizeof = val 
    
  end function psb_dbaseprc_sizeof

  function psb_zbaseprc_sizeof(prec)
    use psb_base_mod
    type(psb_zbaseprc_type), intent(in) :: prec
    integer             :: psb_zbaseprc_sizeof
    integer             :: val,i
    
    val = 0
    if (allocated(prec%iprcparm)) then 
      val = val + 4 * size(prec%iprcparm)
      if (prec%iprcparm(prec_status_) == prec_built) then 
        select case(prec%iprcparm(f_type_)) 
        case(f_ilu_n_,f_ilu_e_)
          ! do nothing
        case(f_slu_)
          write(0,*) 'Should implement check for size of SuperLU data structs'
        case(f_umf_)
          write(0,*) 'Should implement check for size of UMFPACK data structs'
        case(f_slud_)
          write(0,*) 'Should implement check for size of SuperLUDist data structs'
        case default
        end select
        
      end if
    end if
    if (allocated(prec%dprcparm)) val = val + 8 * size(prec%dprcparm)
    if (allocated(prec%d))        val = val + 16 * size(prec%d)
    if (allocated(prec%perm))     val = val + 4 * size(prec%perm)
    if (allocated(prec%invperm))  val = val + 4 * size(prec%invperm)
                                  val = val + psb_sizeof(prec%desc_data)
    if (allocated(prec%av))  then 
      do i=1,size(prec%av)
        val = val + psb_sizeof(prec%av(i))
      end do
    end if
    
    psb_zbaseprc_sizeof = val 
    
  end function psb_zbaseprc_sizeof
    


  subroutine psb_out_prec_descr(p)
    use psb_base_mod
    type(psb_dprec_type), intent(in) :: p
    call psb_file_prec_descr(6,p)
  end subroutine psb_out_prec_descr

  subroutine psb_zout_prec_descr(p)
    use psb_base_mod
    type(psb_zprec_type), intent(in) :: p
    call psb_zfile_prec_descr(6,p)
  end subroutine psb_zout_prec_descr

  subroutine psb_file_prec_descr(iout,p)
    use psb_base_mod
    integer, intent(in)              :: iout
    type(psb_dprec_type), intent(in) :: p
    integer  :: ilev

    write(iout,*) 'Preconditioner description'
    if (allocated(p%baseprecv)) then 
      if (size(p%baseprecv)>=1) then 
        write(iout,*) 'Base preconditioner'
        select case(p%baseprecv(1)%iprcparm(p_type_))
        case(noprec_)
          write(iout,*) 'No preconditioning'
        case(diag_)
          write(iout,*) 'Diagonal scaling'
        case(bjac_)
          write(iout,*) 'Block Jacobi with: ',&
               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
        case(asm_,ras_,ash_,rash_)
          write(iout,*) 'Additive Schwarz with: ',&
               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
          write(iout,*) 'Overlap:',&
               &  p%baseprecv(1)%iprcparm(n_ovr_)
          write(iout,*) 'Restriction: ',&
               &  restrict_names(p%baseprecv(1)%iprcparm(restr_))
          write(iout,*) 'Prolongation: ',&
               &  prolong_names(p%baseprecv(1)%iprcparm(prol_))
        end select
      end if
      if (size(p%baseprecv)>=2) then 
        do ilev = 2, size(p%baseprecv) 
          if (.not.allocated(p%baseprecv(ilev)%iprcparm)) then 
            write(iout,*) 'Inconsistent MLPREC part!'
            return
          endif
          
          write(iout,*) 'Multilevel: Level No', ilev
          write(iout,*) 'Multilevel type: ',&
               &   ml_names(p%baseprecv(ilev)%iprcparm(ml_type_))
          if (p%baseprecv(ilev)%iprcparm(ml_type_)>no_ml_) then 
            write(iout,*) 'Multilevel aggregation: ', &
                 &   aggr_names(p%baseprecv(ilev)%iprcparm(aggr_alg_))
            write(iout,*) 'Smoother:               ', &
                 &  smooth_kinds(p%baseprecv(ilev)%iprcparm(smth_kind_))
            if (p%baseprecv(ilev)%iprcparm(smth_kind_) /= no_smth_) then 
              write(iout,*) 'Smoothing omega: ', &
                   & p%baseprecv(ilev)%dprcparm(smooth_omega_)
              write(iout,*) 'Smoothing position: ',&
                   & smooth_names(p%baseprecv(ilev)%iprcparm(smth_pos_))
            end if
            write(iout,*) 'Coarse matrix: ',&
                 & matrix_names(p%baseprecv(ilev)%iprcparm(coarse_mat_))
            if (allocated(p%baseprecv(ilev)%nlaggr)) then 
              write(iout,*) 'Aggregation sizes: ', &
                   &  sum( p%baseprecv(ilev)%nlaggr(:)),' : ',p%baseprecv(ilev)%nlaggr(:)
            end if
            write(iout,*) 'Factorization type: ',&
                 & fact_names(p%baseprecv(ilev)%iprcparm(f_type_))
            select case(p%baseprecv(ilev)%iprcparm(f_type_))
            case(f_ilu_n_)      
              write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(ilu_fill_in_)
            case(f_ilu_e_)         
              write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(fact_eps_)
            case(f_slu_,f_umf_,f_slud_) 
            case default
              write(iout,*) 'Should never get here!'
            end select
            write(iout,*) 'Number of Jacobi sweeps: ', &
                 &   (p%baseprecv(ilev)%iprcparm(jac_sweeps_))
          end if
        end do
      end if

    else
      write(iout,*) 'No Base preconditioner available, something is wrong!'
      return
    endif

  end subroutine psb_file_prec_descr

  function  psb_prec_short_descr(p)
    use psb_base_mod
    type(psb_dprec_type), intent(in) :: p
    character(len=20) :: psb_prec_short_descr
    psb_prec_short_descr = ' '
!!$    write(iout,*) 'Preconditioner description'
!!$    if (associated(p%baseprecv)) then 
!!$      if (size(p%baseprecv)>=1) then 
!!$        write(iout,*) 'Base preconditioner'
!!$        select case(p%baseprecv(1)%iprcparm(p_type_))
!!$        case(noprec_)
!!$          write(iout,*) 'No preconditioning'
!!$        case(diag_)
!!$          write(iout,*) 'Diagonal scaling'
!!$        case(bjac_)
!!$          write(iout,*) 'Block Jacobi with: ',&
!!$               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
!!$        case(asm_,ras_,ash_,rash_)
!!$          write(iout,*) 'Additive Schwarz with: ',&
!!$               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
!!$          write(iout,*) 'Overlap:',&
!!$               &  p%baseprecv(1)%iprcparm(n_ovr_)
!!$          write(iout,*) 'Restriction: ',&
!!$               &  restrict_names(p%baseprecv(1)%iprcparm(restr_))
!!$          write(iout,*) 'Prolongation: ',&
!!$               &  prolong_names(p%baseprecv(1)%iprcparm(prol_))
!!$        end select
!!$      end if
!!$      if (size(p%baseprecv)>=2) then 
!!$        if (.not.associated(p%baseprecv(2)%iprcparm)) then 
!!$          write(iout,*) 'Inconsistent MLPREC part!'
!!$          return
!!$        endif
!!$        write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_))
!!$        if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then 
!!$          write(iout,*) 'Multilevel aggregation: ', &
!!$               &   aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
!!$          write(iout,*) 'Smoother:               ', &
!!$               &  smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_))
!!$          write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_)
!!$          write(iout,*) 'Smoothing position: ',&
!!$               & smooth_names(p%baseprecv(2)%iprcparm(smth_pos_))
!!$          write(iout,*) 'Coarse matrix: ',&
!!$               & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
!!$          write(iout,*) 'Factorization type: ',&
!!$               & fact_names(p%baseprecv(2)%iprcparm(f_type_))
!!$          select case(p%baseprecv(2)%iprcparm(f_type_))
!!$          case(f_ilu_n_)      
!!$            write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(ilu_fill_in_)
!!$          case(f_ilu_e_)         
!!$            write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_)
!!$          case(f_slu_,f_umf_,f_slud_)         
!!$          case default
!!$            write(iout,*) 'Should never get here!'
!!$          end select
!!$          write(iout,*) 'Number of Jacobi sweeps: ', &
!!$               &   (p%baseprecv(2)%iprcparm(jac_sweeps_))
!!$
!!$        end if
!!$      end if
!!$
!!$    else
!!$      write(iout,*) 'No Base preconditioner available, something is wrong!'
!!$      return
!!$    endif

  end function psb_prec_short_descr


  subroutine psb_zfile_prec_descr(iout,p)
    use psb_base_mod
    integer, intent(in)              :: iout
    type(psb_zprec_type), intent(in) :: p

    write(iout,*) 'Preconditioner description'
    if (allocated(p%baseprecv)) then 
      if (size(p%baseprecv)>=1) then 
        write(iout,*) 'Base preconditioner'
        select case(p%baseprecv(1)%iprcparm(p_type_))
        case(noprec_)
          write(iout,*) 'No preconditioning'
        case(diag_)
          write(iout,*) 'Diagonal scaling'
        case(bjac_)
          write(iout,*) 'Block Jacobi with: ',&
               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
        case(asm_,ras_,ash_,rash_)
          write(iout,*) 'Additive Schwarz with: ',&
               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
          write(iout,*) 'Overlap:',&
               &  p%baseprecv(1)%iprcparm(n_ovr_)
          write(iout,*) 'Restriction: ',&
               &  restrict_names(p%baseprecv(1)%iprcparm(restr_))
          write(iout,*) 'Prolongation: ',&
               &  prolong_names(p%baseprecv(1)%iprcparm(prol_))
        end select
      end if
      if (size(p%baseprecv)>=2) then 
        if (.not.allocated(p%baseprecv(2)%iprcparm)) then 
          write(iout,*) 'Inconsistent MLPREC part!'
          return
        endif
        write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_))
        if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then 
          write(iout,*) 'Multilevel aggregation: ', &
               &   aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
          write(iout,*) 'Smoother:               ', &
               &  smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_))
          if (p%baseprecv(2)%iprcparm(smth_kind_) /= no_smth_) then 
            write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_)
            write(iout,*) 'Smoothing position: ',&
                 & smooth_names(p%baseprecv(2)%iprcparm(smth_pos_))
          end if
          
          write(iout,*) 'Coarse matrix: ',&
               & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
          if (allocated(p%baseprecv(ilev)%nlaggr)) then 
            write(iout,*) 'Aggregation sizes: ', &
                 &  sum( p%baseprecv(2)%nlaggr(:)),' : ',p%baseprecv(2)%nlaggr(:)
          endif
          write(iout,*) 'Factorization type: ',&
               & fact_names(p%baseprecv(2)%iprcparm(f_type_))
          select case(p%baseprecv(2)%iprcparm(f_type_))
          case(f_ilu_n_)      
            write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(ilu_fill_in_)
          case(f_ilu_e_)         
            write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_)
          case(f_slu_,f_umf_,f_slud_)         
          case default
            write(iout,*) 'Should never get here!'
          end select
          write(iout,*) 'Number of Jacobi sweeps: ', &
               &   (p%baseprecv(2)%iprcparm(jac_sweeps_))

        end if
      end if

    else
      write(iout,*) 'No Base preconditioner available, something is wrong!'
      return
    endif

  end subroutine psb_zfile_prec_descr

  function  psb_zprec_short_descr(p)
    use psb_base_mod
    type(psb_zprec_type), intent(in) :: p
    character(len=20) :: psb_zprec_short_descr
    psb_zprec_short_descr = ' '
!!$    write(iout,*) 'Preconditioner description'
!!$    if (associated(p%baseprecv)) then 
!!$      if (size(p%baseprecv)>=1) then 
!!$        write(iout,*) 'Base preconditioner'
!!$        select case(p%baseprecv(1)%iprcparm(p_type_))
!!$        case(noprec_)
!!$          write(iout,*) 'No preconditioning'
!!$        case(diag_)
!!$          write(iout,*) 'Diagonal scaling'
!!$        case(bjac_)
!!$          write(iout,*) 'Block Jacobi with: ',&
!!$               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
!!$        case(asm_,ras_,ash_,rash_)
!!$          write(iout,*) 'Additive Schwarz with: ',&
!!$               &  fact_names(p%baseprecv(1)%iprcparm(f_type_))
!!$          write(iout,*) 'Overlap:',&
!!$               &  p%baseprecv(1)%iprcparm(n_ovr_)
!!$          write(iout,*) 'Restriction: ',&
!!$               &  restrict_names(p%baseprecv(1)%iprcparm(restr_))
!!$          write(iout,*) 'Prolongation: ',&
!!$               &  prolong_names(p%baseprecv(1)%iprcparm(prol_))
!!$        end select
!!$      end if
!!$      if (size(p%baseprecv)>=2) then 
!!$        if (.not.associated(p%baseprecv(2)%iprcparm)) then 
!!$          write(iout,*) 'Inconsistent MLPREC part!'
!!$          return
!!$        endif
!!$        write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_))
!!$        if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then 
!!$          write(iout,*) 'Multilevel aggregation: ', &
!!$               &   aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
!!$          write(iout,*) 'Smoother:               ', &
!!$               &  smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_))
!!$          write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_)
!!$          write(iout,*) 'Smoothing position: ',&
!!$               & smooth_names(p%baseprecv(2)%iprcparm(smth_pos_))
!!$          write(iout,*) 'Coarse matrix: ',&
!!$               & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
!!$          write(iout,*) 'Factorization type: ',&
!!$               & fact_names(p%baseprecv(2)%iprcparm(f_type_))
!!$          select case(p%baseprecv(2)%iprcparm(f_type_))
!!$          case(f_ilu_n_)      
!!$            write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(ilu_fill_in_)
!!$          case(f_ilu_e_)         
!!$            write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_)
!!$          case(f_slu_,f_umf_,f_slud_)         
!!$          case default
!!$            write(iout,*) 'Should never get here!'
!!$          end select
!!$          write(iout,*) 'Number of Jacobi sweeps: ', &
!!$               &   (p%baseprecv(2)%iprcparm(jac_sweeps_))
!!$
!!$        end if
!!$      end if
!!$
!!$    else
!!$      write(iout,*) 'No Base preconditioner available, something is wrong!'
!!$      return
!!$    endif

  end function psb_zprec_short_descr




  function is_legal_base_prec(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_base_prec

    is_legal_base_prec = ((ip>=noprec_).and.(ip<=rash_))
    return
  end function is_legal_base_prec
  function is_legal_n_ovr(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_n_ovr

    is_legal_n_ovr = (ip >=0) 
    return
  end function is_legal_n_ovr
  function is_legal_renum(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_renum
    ! For the time being we are disabling renumbering options. 
    is_legal_renum = (ip ==0) 
    return
  end function is_legal_renum
  function is_legal_jac_sweeps(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_jac_sweeps

    is_legal_jac_sweeps = (ip >= 1) 
    return
  end function is_legal_jac_sweeps
  function is_legal_prolong(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_prolong
    is_legal_prolong = ((ip>=psb_none_).and.(ip<=psb_square_root_))
    return
  end function is_legal_prolong
  function is_legal_restrict(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_restrict
    is_legal_restrict = ((ip==psb_nohalo_).or.(ip==psb_halo_))
    return
  end function is_legal_restrict
  function is_legal_ml_type(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_type

    is_legal_ml_type = ((ip>=no_ml_).and.(ip<=max_ml_))
    return
  end function is_legal_ml_type
  function is_legal_ml_aggr_kind(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_aggr_kind

    is_legal_ml_aggr_kind = ((ip>=loc_aggr_).and.(ip<=max_aggr_))
    return
  end function is_legal_ml_aggr_kind
  function is_legal_ml_smooth_pos(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_smooth_pos

    is_legal_ml_smooth_pos = ((ip>=pre_smooth_).and.(ip<=max_smooth_))
    return
  end function is_legal_ml_smooth_pos
  function is_legal_ml_smth_kind(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_smth_kind

    is_legal_ml_smth_kind = ((ip>=no_smth_).and.(ip<=smth_biz_))
    return
  end function is_legal_ml_smth_kind
  function is_legal_ml_coarse_mat(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_coarse_mat

    is_legal_ml_coarse_mat = ((ip>=mat_distr_).and.(ip<=mat_repl_))
    return
  end function is_legal_ml_coarse_mat
  function is_legal_ml_fact(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_fact

    is_legal_ml_fact = ((ip>=f_ilu_n_).and.(ip<=f_slud_))
    return
  end function is_legal_ml_fact
  function is_legal_ml_lev(ip)
    use psb_base_mod
    integer, intent(in) :: ip
    logical             :: is_legal_ml_lev

    is_legal_ml_lev = (ip>=0)
    return
  end function is_legal_ml_lev
  function is_legal_omega(ip)
    use psb_base_mod
    real(kind(1.d0)), intent(in) :: ip
    logical             :: is_legal_omega

    is_legal_omega = ((ip>=0.0d0).and.(ip<=2.0d0))
    return
  end function is_legal_omega
  function is_legal_ml_eps(ip)
    use psb_base_mod
    real(kind(1.d0)), intent(in) :: ip
    logical             :: is_legal_ml_eps

    is_legal_ml_eps = (ip>=0.0d0)
    return
  end function is_legal_ml_eps


  subroutine psb_icheck_def(ip,name,id,is_legal)
    use psb_base_mod
    integer, intent(inout) :: ip
    integer, intent(in)    :: id
    character(len=*), intent(in) :: name
    interface 
      function is_legal(i)
        integer, intent(in) :: i
        logical             :: is_legal
      end function is_legal
    end interface

    if (.not.is_legal(ip)) then     
      write(0,*) 'Illegal value for ',name,' :',ip, '. defaulting to ',id
      ip = id
    end if
  end subroutine psb_icheck_def

  subroutine psb_dcheck_def(ip,name,id,is_legal)
    use psb_base_mod
    real(kind(1.d0)), intent(inout) :: ip
    real(kind(1.d0)), intent(in)    :: id
    character(len=*), intent(in) :: name
    interface 
      function is_legal(i)
        real(kind(1.d0)), intent(in) :: i
        logical             :: is_legal
      end function is_legal
    end interface

    if (.not.is_legal(ip)) then     
      write(0,*) 'Illegal value for ',name,' :',ip, '. defaulting to ',id
      ip = id
    end if
  end subroutine psb_dcheck_def

  subroutine psb_dbase_precfree(p,info)
    use psb_base_mod

    type(psb_dbaseprc_type), intent(inout) :: p
    integer, intent(out)                :: info
    integer :: i

    info = 0

    ! Actually we migh just deallocate the top level array, except 
    ! for the inner UMFPACK or SLU stuff

    if (allocated(p%d)) then 
      deallocate(p%d,stat=info)
    end if

    if (allocated(p%av))  then 
      do i=1,size(p%av) 
        call psb_sp_free(p%av(i),info)
        if (info /= 0) then 
          ! Actually, we don't care here about this.
          ! Just let it go.
          ! return
        end if
      enddo
      deallocate(p%av,stat=info)
    end if

    if (allocated(p%desc_data%matrix_data)) &
         & call psb_cdfree(p%desc_data,info)
    if (allocated(p%desc_ac%matrix_data)) &
         & call psb_cdfree(p%desc_ac,info)
    
    if (allocated(p%dprcparm)) then 
      deallocate(p%dprcparm,stat=info)
    end if
    ! This is a pointer to something else, must not free it here. 
    nullify(p%base_a) 
    ! This is a pointer to something else, must not free it here. 
    nullify(p%base_desc) 

    if (allocated(p%dorig)) then 
      deallocate(p%dorig,stat=info)
    endif

    if (allocated(p%mlia)) then 
      deallocate(p%mlia,stat=info)
    endif

    if (allocated(p%nlaggr)) then 
      deallocate(p%nlaggr,stat=info)
    endif

    if (allocated(p%perm)) then 
      deallocate(p%perm,stat=info)
    endif

    if (allocated(p%invperm)) then 
      deallocate(p%invperm,stat=info)
    endif

    if (allocated(p%iprcparm)) then 
      if (p%iprcparm(f_type_)==f_slu_) then 
        call psb_dslu_free(p%iprcparm(slu_ptr_),info)
      end if
      if (p%iprcparm(f_type_)==f_slud_) then 
        call psb_dsludist_free(p%iprcparm(slud_ptr_),info)
      end if
      if (p%iprcparm(f_type_)==f_umf_) then 
        call psb_dumf_free(p%iprcparm(umf_symptr_),&
             & p%iprcparm(umf_numptr_),info)
      end if
      deallocate(p%iprcparm,stat=info)
    end if
    call psb_nullify_baseprec(p)
  end subroutine psb_dbase_precfree

  subroutine psb_nullify_dbaseprec(p)
    use psb_base_mod

    type(psb_dbaseprc_type), intent(inout) :: p

    nullify(p%base_a) 
    nullify(p%base_desc) 
!!$    nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,&
!!$         & p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac)

  end subroutine psb_nullify_dbaseprec

  subroutine psb_zbase_precfree(p,info)
    use psb_base_mod
    type(psb_zbaseprc_type), intent(inout) :: p
    integer, intent(out)                :: info
    integer :: i

    info = 0

    if (allocated(p%d)) then 
      deallocate(p%d,stat=info)
    end if

    if (allocated(p%av))  then 
      do i=1,size(p%av) 
        call psb_sp_free(p%av(i),info)
        if (info /= 0) then 
          ! Actually, we don't care here about this.
          ! Just let it go.
          ! return
        end if
      enddo
      deallocate(p%av,stat=info)

    end if
    ! call psb_cdfree(p%desc_data,info)
    ! call psb_cdfree(p%desc_ac,info)

    if (allocated(p%dprcparm)) then 
      deallocate(p%dprcparm,stat=info)
    end if
    ! This is a pointer to something else, must not free it here. 
    nullify(p%base_a) 
    ! This is a pointer to something else, must not free it here. 
    nullify(p%base_desc) 

    if (allocated(p%dorig)) then 
      deallocate(p%dorig,stat=info)
    endif

    if (allocated(p%mlia)) then 
      deallocate(p%mlia,stat=info)
    endif

    if (allocated(p%nlaggr)) then 
      deallocate(p%nlaggr,stat=info)
    endif

    if (allocated(p%perm)) then 
      deallocate(p%perm,stat=info)
    endif

    if (allocated(p%invperm)) then 
      deallocate(p%invperm,stat=info)
    endif

    if (allocated(p%iprcparm)) then 
      if (p%iprcparm(f_type_)==f_slu_) then 
        call psb_zslu_free(p%iprcparm(slu_ptr_),info)
      end if
      if (p%iprcparm(f_type_)==f_umf_) then 
        call psb_zumf_free(p%iprcparm(umf_symptr_),&
             & p%iprcparm(umf_numptr_),info)
      end if
      deallocate(p%iprcparm,stat=info)
    end if
    call psb_nullify_baseprec(p)
  end subroutine psb_zbase_precfree

  subroutine psb_nullify_zbaseprec(p)
    use psb_base_mod

    type(psb_zbaseprc_type), intent(inout) :: p


    nullify(p%base_a) 
    nullify(p%base_desc) 

  end subroutine psb_nullify_zbaseprec


  function pr_to_str(iprec)
    use psb_base_mod

    integer, intent(in)  :: iprec
    character(len=10)     :: pr_to_str

    select case(iprec)
    case(noprec_)
      pr_to_str='NOPREC'
    case(diag_)         
      pr_to_str='DIAG'
    case(bjac_)         
      pr_to_str='BJAC'
    case(asm_)      
      pr_to_str='ASM'
    case(ash_)      
      pr_to_str='ASM'
    case(ras_)      
      pr_to_str='ASM'
    case(rash_)      
      pr_to_str='ASM'
    end select

  end function pr_to_str

end module psb_prec_type