Fix cmp_avg_cr

Fix missing library in C interface linking

.VERSION

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

CMakeLists.txt

@@ -13,15 +13,10 @@ endif()
 
 
 # Check for the installation path for psblas
-#if(NOT DEFINED PSBLAS_INSTALL_DIR)
-#    message(FATAL_ERROR "Please specify the path to the psblas installation directory using -DPSBLAS_INSTALL_DIR=<path>")
-#endif()
-message(STATUS "psblas directory is ${PSBLAS_INSTALL_DIR};;")
-
+message(STATUS "psblas directory is ${PSBLAS_INSTALL_DIR};")
 
 message(STATUS "PSBLAS DIRECTORY INC ${INCDIR}; MOD ${MODDIR}; LIB ${LIBDIR};")
 
-
 #set(CMAKE_CXX_STANDARD 17) # Set cxx standard for the c++ part of the library
 
 
@@ -31,7 +26,7 @@ find_package(psblas REQUIRED PATHS ${PSBLAS_INSTALL_DIR})
 if(NOT psblas_FOUND)
     message(FATAL_ERROR "PSBLAS not found!")
 else()
-    message(STATUS "Found PSBLAS: ${psblas_LIBRARIES}")
+    message(STATUS "Found PSBLAS: ${PSBLAS_LIBRARIES}")
 endif()
 
 if(CMAKE_BUILD_TYPE STREQUAL "Debug")
@@ -42,9 +37,9 @@ if(CMAKE_BUILD_TYPE STREQUAL "Debug")
     message(STATUS "Fortran and CXX debug flags added: -g")
 endif()
 
-    string(APPEND CMAKE_Fortran_FLAGS " -O2")
-    string(APPEND CMAKE_CXX_FLAGS " -O2")
-    message(STATUS "Fortran and CXX optimization flags added: -O2")
+string(APPEND CMAKE_Fortran_FLAGS " -O2")
+string(APPEND CMAKE_CXX_FLAGS " -O2")
+message(STATUS "Fortran and CXX optimization flags added: -O2")
 
 
 
@@ -430,7 +425,9 @@ target_include_directories(amgcbind PUBLIC ${INCDIR} ${MODDIR})
 target_link_libraries(amgcbind
   #PUBLIC ${LAPACK_LINKER_FLAGS} ${LAPACK_LIBRARIES} ${LAPACK95_LIBRARIES}
   #PUBLIC ${BLAS_LINKER_FLAGS} ${BLAS_LIBRARIES} ${BLAS95_LIBRARIES}
-  PUBLIC amgprec psblas::util psblas::linsolve psblas::prec psblas::ext psblas::cbind psblas::base)  #TODO  check actual libraries needed
+  PUBLIC amgprec psblas::util psblas::linsolve psblas::prec
+  psblas::ext psblas::cbind psblas::base)
+  #TODO  check actual libraries needed
 
 
 

LICENSE

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

Make.inc.in

@@ -71,7 +71,7 @@ EXTRALIBS=@EXTRA_LIBS@
 
 #
 AMGCDEFINES=$(SLUFLAGS) $(UMFFLAGS) $(SLUDISTFLAGS)
-#$(PSBCDEFINES) $(MUMPSFLAGS)
+#$(PSBCDEFINES) $(MUMPSFLAGS) 
 CDEFINES=$(AMGCDEFINES)
 AMGFDEFINES=@AMGFDEFINES@ $(PSBFDEFINES)
 FDEFINES=$(AMGFDEFINES)

Makefile

@@ -19,7 +19,7 @@ libdir:
 amgobjs: mods
 	cd amgprec && $(MAKE) objs
 
-cbnd: mods
+cbnd: amgobjs
 	cd cbind && $(MAKE)  objs
 
 install: all

amgprec/Makefile

@@ -58,7 +58,8 @@ MODOBJS=amg_base_prec_type.o amg_prec_type.o  amg_prec_mod.o \
 	$(SMODOBJS) $(DMODOBJS) $(CMODOBJS) $(ZMODOBJS)
 
 
-LOCAL_MODS=$(MODOBJS:.o=$(.mod))
+LOCAL_MODS=$(MODOBJS:.o=$(.mod)) amg_c_l1_diag_solver$(.mod) amg_s_l1_diag_solver$(.mod) \
+	amg_d_l1_diag_solver$(.mod) amg_z_l1_diag_solver$(.mod)
 LIBNAME=libamg_prec.a
 
 all: mods  objs impld
@@ -70,7 +71,7 @@ objs: mods impld
 impld: mods
 	cd impl && $(MAKE) 
 
-lib: mods impld
+lib: objs
 	cd impl && $(MAKE) lib
 	$(AR) $(HERE)/$(LIBNAME) $(MODOBJS)
 	$(RANLIB) $(HERE)/$(LIBNAME)
@@ -79,7 +80,7 @@ lib: mods impld
 
 $(MODOBJS): $(PSBLAS_MODDIR)/$(PSBBASEMODNAME)$(.mod)
 
-#amg_base_prec_type.o: amg_const.h
+amg_base_prec_type.o: amg_config.h
 amg_s_prec_type.o amg_d_prec_type.o amg_c_prec_type.o amg_z_prec_type.o	: amg_base_prec_type.o
 amg_prec_type.o: amg_s_prec_type.o amg_d_prec_type.o amg_c_prec_type.o amg_z_prec_type.o
 amg_prec_mod.o: amg_prec_type.o amg_s_prec_mod.o amg_d_prec_mod.o amg_c_prec_mod.o amg_z_prec_mod.o

amgprec/amg_base_prec_type.F90

@@ -84,7 +84,7 @@ module amg_base_prec_type
   character(len=*), parameter   :: amg_version_string_ = "1.2.0"
   integer(psb_ipk_), parameter  :: amg_version_major_  = 1
   integer(psb_ipk_), parameter  :: amg_version_minor_  = 2
-  integer(psb_ipk_), parameter  :: amg_patchlevel_     = 0
+  integer(psb_ipk_), parameter  :: amg_patchlevel_     = 1
 
   type amg_ml_parms
     integer(psb_ipk_) :: sweeps_pre, sweeps_post

amgprec/amg_c_ainv_solver.F90

@@ -103,7 +103,7 @@ module amg_c_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_c_ainv_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -225,7 +225,7 @@ module amg_c_ainv_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_d_vect_type, psb_c_base_vect_type, psb_spk_, psb_ipk_
       implicit none
-      type(psb_cspmat_type), intent(in), target   :: a
+      type(psb_cspmat_type), intent(inout), target  :: a
       integer(psb_ipk_), intent(in)               :: fillin,alg
       real(psb_spk_), intent(in)                  :: thresh
       type(psb_cspmat_type), intent(inout)        :: wmat, zmat

amgprec/amg_c_as_smoother.f90

@@ -230,7 +230,7 @@ module amg_c_as_smoother
            & psb_desc_type, psb_c_base_sparse_mat, psb_ipk_,&
            & psb_i_base_vect_type
       implicit none 
-      type(psb_cspmat_type), intent(in), target        :: a
+      type(psb_cspmat_type), intent(inout), target     :: a
       Type(psb_desc_type), Intent(inout)                 :: desc_a 
       class(amg_c_as_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                     :: info

amgprec/amg_c_base_aggregator_mod.f90

@@ -89,7 +89,7 @@ module amg_c_base_aggregator_mod
     procedure, pass(ag) :: bld_tprol   => amg_c_base_aggregator_build_tprol
     procedure, pass(ag) :: mat_bld     => amg_c_base_aggregator_mat_bld
     procedure, pass(ag) :: mat_asb     => amg_c_base_aggregator_mat_asb
-    procedure, pass(ag) :: bld_map     => amg_c_base_aggregator_bld_map
+    procedure, pass(ag) :: bld_linmap  => amg_c_base_aggregator_bld_linmap
     procedure, pass(ag) :: update_next => amg_c_base_aggregator_update_next
     procedure, pass(ag) :: clone       => amg_c_base_aggregator_clone
     procedure, pass(ag) :: free        => amg_c_base_aggregator_free
@@ -458,7 +458,7 @@ contains
   end subroutine amg_c_base_aggregator_mat_asb
 
   !
-  !> Function   bld_map
+  !> Function   bld_linmap
   !! \memberof  amg_c_base_aggregator_type
   !! \brief     Build linear map between hierarchy levels 
   !!
@@ -473,7 +473,7 @@ contains
   !!  \param map       The output map
   !!  \param info    Return code
   !!  
-  subroutine  amg_c_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
+  subroutine  amg_c_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
        & op_restr,op_prol,map,info)
     use psb_base_mod
     implicit none
@@ -484,7 +484,7 @@ contains
     type(psb_clinmap_type), intent(out)    :: map
     integer(psb_ipk_), intent(out)       :: info
     integer(psb_ipk_) :: err_act
-    character(len=20) :: name='c_base_aggregator_bld_map'
+    character(len=20) :: name='c_base_aggregator_bld_linmap'
 
     info = psb_success_
     call psb_erractionsave(err_act)
@@ -508,7 +508,6 @@ contains
 9999 call psb_error_handler(err_act)
 
     return
-  end subroutine amg_c_base_aggregator_bld_map
-
+  end subroutine amg_c_base_aggregator_bld_linmap
 
 end module amg_c_base_aggregator_mod

amgprec/amg_c_base_smoother_mod.f90

@@ -237,7 +237,7 @@ module amg_c_base_smoother_mod
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & amg_c_base_smoother_type, psb_ipk_, psb_i_base_vect_type
       ! Arguments
-      type(psb_cspmat_type), intent(in), target     :: a
+      type(psb_cspmat_type), intent(inout), target  :: a
       Type(psb_desc_type), Intent(inout)              :: desc_a 
       class(amg_c_base_smoother_type), intent(inout) :: sm 
       integer(psb_ipk_), intent(out)                  :: info

amgprec/amg_c_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_c_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_cspmat_type), intent(in), target             :: a
+      type(psb_cspmat_type), intent(inout), target          :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_c_base_solver_type), intent(inout)          :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_c_diag_solver.f90

@@ -119,7 +119,7 @@ module amg_c_diag_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & amg_c_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_c_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_c_l1_diag_solver
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & amg_c_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_c_l1_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_c_gs_solver.f90

@@ -181,7 +181,7 @@ module amg_c_gs_solver
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_gs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_c_gs_solver
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_bwgs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_cspmat_type), intent(in), target           :: a
+    type(psb_cspmat_type), intent(inout), target        :: a
     Type(psb_desc_type), Intent(inout)                  :: desc_a 
     class(amg_c_id_solver_type), intent(inout)          :: sv
     integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_ilu_solver.f90

@@ -144,7 +144,7 @@ module amg_c_ilu_solver
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_ilu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_inner_mod.f90

@@ -56,7 +56,7 @@ module amg_c_inner_mod
            & psb_spk_, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       import :: amg_cprec_type
       implicit none
-      type(psb_cspmat_type), intent(in), target          :: a
+      type(psb_cspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       type(amg_cprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -67,7 +67,7 @@ module amg_c_inner_mod
   end interface amg_mlprec_bld
 
   interface amg_mlprec_aply
-    subroutine amg_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
+    subroutine amg_cmlprec_aply_a(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_
       import :: amg_cprec_type
       implicit none 
@@ -79,7 +79,7 @@ module amg_c_inner_mod
       character,intent(in)               :: trans
       complex(psb_spk_),target             :: work(:)
       integer(psb_ipk_), intent(out)     :: info
-    end subroutine amg_cmlprec_aply
+    end subroutine amg_cmlprec_aply_a
     subroutine amg_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_cspmat_type, psb_desc_type, &
            & psb_spk_, psb_c_vect_type, psb_ipk_

amgprec/amg_c_invk_solver.f90

@@ -94,7 +94,7 @@ module amg_c_invk_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_c_invk_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_invt_solver.f90

@@ -94,7 +94,7 @@ module amg_c_invt_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_c_invt_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_jac_smoother.f90

@@ -151,7 +151,7 @@ module amg_c_jac_smoother
       import :: psb_desc_type, amg_c_jac_smoother_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_c_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -274,7 +274,7 @@ module amg_c_jac_smoother
       import :: psb_desc_type, amg_c_l1_jac_smoother_type, psb_c_vect_type, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_c_l1_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -345,6 +345,7 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
+      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_c_jac_solver.f90

@@ -143,7 +143,7 @@ module amg_c_jac_solver
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_c_jac_solver
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_l1_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_krm_solver.f90

@@ -174,7 +174,7 @@ module amg_c_krm_solver
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_c_krm_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_c_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_mumps_solver_type), intent(inout)       :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_c_onelev_mod.f90

@@ -154,13 +154,35 @@ module amg_c_onelev_mod
   private :: c_wrk_alloc, c_wrk_free, &
        & c_wrk_clone, c_wrk_move_alloc, c_wrk_cnv, c_wrk_sizeof
 
+  !
+  ! Remap.
+  !  This keeps track of remapping.
+  !  Logic here is as follows:
+  !  1. AC_PRE_REMAP  need to figure out if we
+  !                really need it.
+  !  2. DESC_AC_PRE_REMAP  contains the descriptor before
+  !                remapping. Meaning that it is possible
+  !                to implement the RESTRICTOR operator by
+  !                 a. Doing LINMAP_U2V onto this one
+  !                 b. For each process, send the data to
+  !                    IDEST.
+  !                This assumes that remapping goes by
+  !                a factor of 2.
+  !                For the PROLONGATOR operators, we first
+  !                use DESC_AC, then split and send onto
+  !                the processes in DESC_AC_PRE_REMAP.
+  !
+  ! To be fixed: what happens if NP the starting processes
+  ! is not an even number? Coordinate with _X_remap in PSBLAS
+  ! 
   type amg_c_remap_data_type
     type(psb_cspmat_type)          :: ac_pre_remap
     type(psb_desc_type)            :: desc_ac_pre_remap
     integer(psb_ipk_)              :: idest
     integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
   contains
-    procedure, pass(rmp) :: clone   => c_remap_data_clone
+    procedure, pass(rmp) :: clone      => c_remap_data_clone
+    procedure, pass(rmp) :: move_alloc => c_remap_move_alloc
   end type amg_c_remap_data_type
 
   type amg_c_onelev_type
@@ -206,7 +228,7 @@ module amg_c_onelev_mod
     procedure, pass(lv) :: get_wrksz => c_base_onelev_get_wrksize
     procedure, pass(lv) :: allocate_wrk   => c_base_onelev_allocate_wrk
     procedure, pass(lv) :: free_wrk       => c_base_onelev_free_wrk
-    procedure, nopass   :: stringval => amg_stringval
+    procedure, nopass   :: stringval  => amg_stringval
     procedure, pass(lv) :: move_alloc => c_base_onelev_move_alloc
 
 
@@ -618,7 +640,7 @@ contains
     ! Arguments
     class(amg_c_onelev_type), target, intent(inout) :: lv
     class(amg_c_onelev_type), target, intent(inout) :: lvout
-    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_ipk_), intent(out)                  :: info
 
     info = psb_success_
     if (allocated(lv%sm)) then
@@ -684,6 +706,7 @@ contains
     if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info)
     if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info)
     if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info)
+    if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info)
     b%base_a    => lv%base_a
     b%base_desc => lv%base_desc
 
@@ -738,6 +761,7 @@ contains
     info = psb_success_
     nwv = lv%get_wrksz()
     if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
+!!$    write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb()
     if (info == 0) then
       if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
         !
@@ -785,47 +809,32 @@ contains
 
     info = psb_success_
     call wk%free(info)
-    if  (present(desc2)) then
-!!$      write(0,*) 'Check on wrk_alloc 2',&
-!!$           & desc2%get_local_rows(), desc%get_local_rows(),&
-!!$           & desc2%get_local_cols(),desc%get_local_cols()
-!!$      flush(0)
+!!$    write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',&
+!!$         & present(desc2),desc%is_valid()
+
+    allocate(wk%wv(nwv),stat=info)    
+    if  (present(desc2).and.(desc%is_valid())) then
+!!$      write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
       if (desc2%get_local_cols()>desc%get_local_cols()) then
-        call psb_geasb(wk%vx2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc2,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
       else
-!!$        write(0,*) 'Check on wrk_alloc 1.5 ',&
-!!$             & desc%get_local_rows(),&
-!!$             & desc%get_local_cols()
-        call psb_geasb(wk%vx2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)
       end if
-    else
-!!$      write(0,*) 'Check on wrk_alloc 1 ',&
-!!$           & desc%get_local_rows(),&
-!!$           & desc%get_local_cols()
+    else if (present(desc2)) then
+      call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
+    else if (desc%is_valid()) then
+      call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)      
+    end if
+
+  contains
+    subroutine inner_do_wrk_alloc(wk,nwv,desc,vmold)
+      class(amg_cmlprec_wrk_type), target, intent(inout) :: wk
+      integer(psb_ipk_), intent(in)                     :: nwv
+      type(psb_desc_type), intent(in)                   :: desc
+      class(psb_c_base_vect_type), intent(in), optional  :: vmold
+      
+      integer(psb_ipk_) :: i
+      
       call psb_geasb(wk%vx2l,desc,info,&
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vy2l,desc,info,&
@@ -834,12 +843,11 @@ contains
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vty,desc,info,&
            & scratch=.true.,mold=vmold)
-      allocate(wk%wv(nwv),stat=info)
       do i=1,nwv
         call psb_geasb(wk%wv(i),desc,info,&
              & scratch=.true.,mold=vmold)
       end do
-    end if
+    end subroutine inner_do_wrk_alloc
   end subroutine c_wrk_alloc
 
   subroutine c_wrk_free(wk,info)
@@ -992,4 +1000,25 @@ contains
     call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
   end subroutine c_remap_data_clone
 
+  subroutine c_remap_move_alloc(rmp, remap_out, info)
+    use psb_base_mod
+    implicit none
+    ! Arguments
+    class(amg_c_remap_data_type), target, intent(inout) :: rmp
+    class(amg_c_remap_data_type), target, intent(inout) :: remap_out
+    integer(psb_ipk_), intent(out)                    :: info
+    !
+    integer(psb_ipk_) :: i
+
+    info = psb_success_
+
+    call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info)
+    if (info == psb_success_) &
+         & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info)
+    remap_out%idest = rmp%idest
+    call move_alloc(rmp%isrc,remap_out%isrc)
+    call move_alloc(rmp%nrsrc,remap_out%nrsrc)
+    call move_alloc(rmp%naggr,remap_out%naggr)
+  end subroutine c_remap_move_alloc
+  
 end module amg_c_onelev_mod

amgprec/amg_c_prec_type.f90

@@ -102,6 +102,7 @@ module amg_c_prec_type
     ! The multilevel hierarchy
     !
     type(amg_c_onelev_type), allocatable :: precv(:)
+    integer(psb_ipk_)                    :: nlevs
   contains
     procedure, pass(prec)               :: psb_c_apply2_vect => amg_c_apply2_vect
     procedure, pass(prec)               :: psb_c_apply1_vect => amg_c_apply1_vect
@@ -119,6 +120,7 @@ module amg_c_prec_type
     procedure, pass(prec)               :: cmp_complexity => amg_c_cmp_compl
     procedure, pass(prec)               :: get_avg_cr => amg_c_get_avg_cr
     procedure, pass(prec)               :: cmp_avg_cr => amg_c_cmp_avg_cr
+    procedure, pass(prec)               :: set_nlevs  => amg_c_set_nlevs
     procedure, pass(prec)               :: get_nlevs  => amg_c_get_nlevs
     procedure, pass(prec)               :: get_nzeros => amg_c_get_nzeros
     procedure, pass(prec)               :: sizeof => amg_cprec_sizeof
@@ -311,7 +313,7 @@ module amg_c_prec_type
            & psb_c_base_sparse_mat, psb_c_base_vect_type, &
            & psb_i_base_vect_type, amg_cprec_type, psb_ipk_
       implicit none
-      type(psb_cspmat_type), intent(in), target          :: a
+      type(psb_cspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_cprec_type), intent(inout), target       :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -323,14 +325,15 @@ module amg_c_prec_type
   end interface amg_precbld
 
   interface amg_hierarchy_bld
-    subroutine amg_c_hierarchy_bld(a,desc_a,prec,info)
+    subroutine amg_c_hierarchy_bld(a,desc_a,prec,info,cpymat)
       import :: psb_cspmat_type, psb_desc_type, psb_spk_, &
            & amg_cprec_type, psb_ipk_
       implicit none
-      type(psb_cspmat_type), intent(in), target          :: a
+      type(psb_cspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_cprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
+      logical, intent(in), optional                        :: cpymat
       !      character, intent(in),optional             :: upd
     end subroutine amg_c_hierarchy_bld
   end interface amg_hierarchy_bld
@@ -434,10 +437,19 @@ contains
     class(amg_cprec_type), intent(in) :: prec
     integer(psb_ipk_) :: val
     val = 0
-    if (allocated(prec%precv)) then
-      val = size(prec%precv)
-    end if
+!!$    if (allocated(prec%precv)) then
+!!$      val = size(prec%precv)
+!!$    end if
+    val = prec%nlevs
+!!$    write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
   end function amg_c_get_nlevs
+
+  subroutine amg_c_set_nlevs(prec,nl) 
+    implicit none
+    class(amg_cprec_type), intent(inout) :: prec
+    integer(psb_ipk_) :: nl
+    prec%nlevs = nl
+  end subroutine amg_c_set_nlevs
   !
   ! Function returning the size of the amg_prec_type data structure
   ! in bytes or in number of nonzeros of the operator(s) involved.
@@ -508,7 +520,7 @@ contains
 
     real(psb_spk_) :: num, den, nmin
     type(psb_ctxt_type) :: ctxt
-    integer(psb_ipk_)   :: il
+    integer(psb_ipk_)   :: il,nl
 
     num = -sone
     den = sone
@@ -518,7 +530,10 @@ contains
       num = prec%precv(il)%base_a%get_nzeros()
       if (num >= szero) then
         den = num
-        do il=2,size(prec%precv)
+        nl = prec%get_nlevs()
+!!$        write(0,*) 'Inside cmp_compl ',nl,size(prec%precv)
+        do il=2, nl
+!!$          write(0,*) '                 ',il,associated(prec%precv(il)%base_a)
           num = num + max(0,prec%precv(il)%base_a%get_nzeros())
         end do
       end if
@@ -549,7 +564,6 @@ contains
   end function amg_c_get_avg_cr
 
   subroutine amg_c_cmp_avg_cr(prec)
-
     implicit none
     class(amg_cprec_type), intent(inout) :: prec
 
@@ -557,17 +571,18 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_)   :: il, nl, iam, np
 
-
     avgcr = szero
+    nl = prec%get_nlevs()
+    do il=2,nl
+      if (prec%precv(il)%base_desc%is_ok()) then 
+        ctxt = prec%precv(il)%base_desc%get_ctxt()
+        call psb_info(ctxt,iam,np)
+        if (iam >=0)  avgcr = avgcr + max(szero,prec%precv(il)%szratio)
+      end if
+    end do
+    avgcr = avgcr / (nl-1)
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    if (allocated(prec%precv)) then
-      nl = size(prec%precv)
-      do il=2,nl
-        avgcr = avgcr + max(szero,prec%precv(il)%szratio)
-      end do
-      avgcr = avgcr / (nl-1)
-    end if
     call psb_sum(ctxt,avgcr)
     prec%ag_data%avg_cr = avgcr/np
   end subroutine amg_c_cmp_avg_cr
@@ -585,9 +600,7 @@ contains
   !             error code.
   !
   subroutine amg_cprecfree(p,info)
-
     implicit none
-
     ! Arguments
     type(amg_cprec_type), intent(inout) :: p
     integer(psb_ipk_), intent(out)        :: info
@@ -613,9 +626,7 @@ contains
   end subroutine amg_cprecfree
 
   subroutine amg_c_prec_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_cprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -647,9 +658,7 @@ contains
   end subroutine amg_c_prec_free
 
   subroutine amg_c_smoothers_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_cprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -665,12 +674,8 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
-    !
-    ! In the internals, do FREE on components,
-    ! but do not deallocate them
-    !
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         call prec%precv(i)%free_smoothers(info)
       end do
     end if
@@ -684,9 +689,7 @@ contains
   end subroutine amg_c_smoothers_free
 
   subroutine amg_c_hierarchy_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_cprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -713,7 +716,6 @@ contains
 
   end subroutine amg_c_hierarchy_free
 
-
   !
   ! Top level methods.
   !
@@ -778,7 +780,6 @@ contains
 
   end subroutine amg_c_apply1_vect
 
-
   subroutine amg_c_apply2v(prec,x,y,desc_data,info,trans,work)
     implicit none
     type(psb_desc_type),intent(in)    :: desc_data
@@ -843,7 +844,6 @@ contains
   subroutine amg_c_dump(prec,info,istart,iend,iproc,prefix,head,&
        & ac,rp,smoother,solver,tprol,&
        & global_num)
-
     implicit none
     class(amg_cprec_type), intent(in)     :: prec
     integer(psb_ipk_), intent(out)          :: info
@@ -860,7 +860,7 @@ contains
     info = 0
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    iln = size(prec%precv)
+    iln = prec%get_nlevs()
     if (present(istart)) then
       il1 = max(1,istart)
     else
@@ -884,7 +884,6 @@ contains
   end subroutine amg_c_dump
 
   subroutine amg_c_cnv(prec,info,amold,vmold,imold)
-
     implicit none
     class(amg_cprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)       :: info
@@ -896,7 +895,7 @@ contains
 
     info = psb_success_
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         if (info == psb_success_ ) &
              & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
       end do
@@ -905,7 +904,6 @@ contains
   end subroutine amg_c_cnv
 
   subroutine amg_c_clone(prec,precout,info)
-
     implicit none
     class(amg_cprec_type), intent(inout) :: prec
     class(psb_cprec_type), intent(inout) :: precout
@@ -917,7 +915,6 @@ contains
   end subroutine amg_c_clone
 
   subroutine amg_c_inner_clone(prec,precout,info)
-
     implicit none
     class(amg_cprec_type), intent(inout)         :: prec
     class(psb_cprec_type), target, intent(inout) :: precout
@@ -933,8 +930,9 @@ contains
       pout%ctxt          = prec%ctxt
       pout%ag_data       = prec%ag_data
       pout%outer_sweeps  = prec%outer_sweeps
+      pout%nlevs         = prec%nlevs
       if (allocated(prec%precv)) then
-        ln = size(prec%precv)
+        ln = prec%get_nlevs()
         allocate(pout%precv(ln),stat=info)
         if (info /= psb_success_) goto 9999
         if (ln >= 1) then
@@ -942,6 +940,7 @@ contains
         end if
         do lev=2, ln
           if (info /= psb_success_) exit
+!!$          write(0,*) 'Inner_clone must be checked and reimplemented! '
           call prec%precv(lev)%clone(pout%precv(lev),info)
           if (info == psb_success_) then
             pout%precv(lev)%base_a       => pout%precv(lev)%ac
@@ -1006,7 +1005,7 @@ contains
     integer(psb_ipk_), intent(out)        :: info
     class(psb_c_base_vect_type), intent(in), optional  :: vmold
     !
-    ! In MLD the DESC optional argument is ignored, since
+    ! In AMG the DESC optional argument is ignored, since
     ! the necessary info is contained in the various entries of the
     ! PRECV component.
     type(psb_desc_type), intent(in), optional  :: desc
@@ -1021,7 +1020,7 @@ contains
     if (psb_errstatus_fatal()) then
       info = psb_err_internal_error_; goto 9999
     end if
-    nlev   = size(prec%precv)
+    nlev   = prec%get_nlevs()
     level = 1
     do level = 1, nlev
       call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@@ -1044,7 +1043,6 @@ contains
   subroutine amg_c_free_wrk(prec,info)
     use psb_base_mod
     implicit none
-
     ! Arguments
     class(amg_cprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -1061,7 +1059,7 @@ contains
     end if
 
     if (allocated(prec%precv)) then
-      nlev   = size(prec%precv)
+      nlev   = prec%get_nlevs()
       do level = 1, nlev
         call prec%precv(level)%free_wrk(info)
       end do

amgprec/amg_c_slu_solver.F90

@@ -124,7 +124,7 @@ module amg_c_slu_solver
       Implicit None
 
       ! Arguments
-      type(psb_cspmat_type), intent(in), target           :: a
+      type(psb_cspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_c_slu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_ainv_solver.F90

@@ -103,7 +103,7 @@ module amg_d_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_ainv_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -225,7 +225,7 @@ module amg_d_ainv_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
            & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, psb_ipk_
       implicit none
-      type(psb_dspmat_type), intent(in), target   :: a
+      type(psb_dspmat_type), intent(inout), target  :: a
       integer(psb_ipk_), intent(in)               :: fillin,alg
       real(psb_dpk_), intent(in)                  :: thresh
       type(psb_dspmat_type), intent(inout)        :: wmat, zmat

amgprec/amg_d_as_smoother.f90

@@ -230,7 +230,7 @@ module amg_d_as_smoother
            & psb_desc_type, psb_d_base_sparse_mat, psb_ipk_,&
            & psb_i_base_vect_type
       implicit none 
-      type(psb_dspmat_type), intent(in), target        :: a
+      type(psb_dspmat_type), intent(inout), target     :: a
       Type(psb_desc_type), Intent(inout)                 :: desc_a 
       class(amg_d_as_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                     :: info

amgprec/amg_d_base_aggregator_mod.f90

@@ -89,7 +89,7 @@ module amg_d_base_aggregator_mod
     procedure, pass(ag) :: bld_tprol   => amg_d_base_aggregator_build_tprol
     procedure, pass(ag) :: mat_bld     => amg_d_base_aggregator_mat_bld
     procedure, pass(ag) :: mat_asb     => amg_d_base_aggregator_mat_asb
-    procedure, pass(ag) :: bld_map     => amg_d_base_aggregator_bld_map
+    procedure, pass(ag) :: bld_linmap  => amg_d_base_aggregator_bld_linmap
     procedure, pass(ag) :: update_next => amg_d_base_aggregator_update_next
     procedure, pass(ag) :: clone       => amg_d_base_aggregator_clone
     procedure, pass(ag) :: free        => amg_d_base_aggregator_free
@@ -458,7 +458,7 @@ contains
   end subroutine amg_d_base_aggregator_mat_asb
 
   !
-  !> Function   bld_map
+  !> Function   bld_linmap
   !! \memberof  amg_d_base_aggregator_type
   !! \brief     Build linear map between hierarchy levels 
   !!
@@ -473,7 +473,7 @@ contains
   !!  \param map       The output map
   !!  \param info    Return code
   !!  
-  subroutine  amg_d_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
+  subroutine  amg_d_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
        & op_restr,op_prol,map,info)
     use psb_base_mod
     implicit none
@@ -484,7 +484,7 @@ contains
     type(psb_dlinmap_type), intent(out)    :: map
     integer(psb_ipk_), intent(out)       :: info
     integer(psb_ipk_) :: err_act
-    character(len=20) :: name='d_base_aggregator_bld_map'
+    character(len=20) :: name='d_base_aggregator_bld_linmap'
 
     info = psb_success_
     call psb_erractionsave(err_act)
@@ -508,7 +508,6 @@ contains
 9999 call psb_error_handler(err_act)
 
     return
-  end subroutine amg_d_base_aggregator_bld_map
-
+  end subroutine amg_d_base_aggregator_bld_linmap
 
 end module amg_d_base_aggregator_mod

amgprec/amg_d_base_smoother_mod.f90

@@ -237,7 +237,7 @@ module amg_d_base_smoother_mod
            & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
            & amg_d_base_smoother_type, psb_ipk_, psb_i_base_vect_type
       ! Arguments
-      type(psb_dspmat_type), intent(in), target     :: a
+      type(psb_dspmat_type), intent(inout), target  :: a
       Type(psb_desc_type), Intent(inout)              :: desc_a 
       class(amg_d_base_smoother_type), intent(inout) :: sm 
       integer(psb_ipk_), intent(out)                  :: info

amgprec/amg_d_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_d_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_dspmat_type), intent(in), target             :: a
+      type(psb_dspmat_type), intent(inout), target          :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_d_base_solver_type), intent(inout)          :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_d_diag_solver.f90

@@ -119,7 +119,7 @@ module amg_d_diag_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
            & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
            & amg_d_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_d_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_d_l1_diag_solver
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
            & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
            & amg_d_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_d_l1_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_d_gs_solver.f90

@@ -181,7 +181,7 @@ module amg_d_gs_solver
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_gs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_d_gs_solver
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_bwgs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_dspmat_type), intent(in), target           :: a
+    type(psb_dspmat_type), intent(inout), target        :: a
     Type(psb_desc_type), Intent(inout)                  :: desc_a 
     class(amg_d_id_solver_type), intent(inout)          :: sv
     integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_ilu_solver.f90

@@ -144,7 +144,7 @@ module amg_d_ilu_solver
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_ilu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_inner_mod.f90

@@ -56,7 +56,7 @@ module amg_d_inner_mod
            & psb_dpk_, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_
       import :: amg_dprec_type
       implicit none
-      type(psb_dspmat_type), intent(in), target          :: a
+      type(psb_dspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       type(amg_dprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -67,7 +67,7 @@ module amg_d_inner_mod
   end interface amg_mlprec_bld
 
   interface amg_mlprec_aply
-    subroutine amg_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
+    subroutine amg_dmlprec_aply_a(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
       import :: amg_dprec_type
       implicit none 
@@ -79,7 +79,7 @@ module amg_d_inner_mod
       character,intent(in)               :: trans
       real(psb_dpk_),target             :: work(:)
       integer(psb_ipk_), intent(out)     :: info
-    end subroutine amg_dmlprec_aply
+    end subroutine amg_dmlprec_aply_a
     subroutine amg_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_dspmat_type, psb_desc_type, &
            & psb_dpk_, psb_d_vect_type, psb_ipk_

amgprec/amg_d_invk_solver.f90

@@ -94,7 +94,7 @@ module amg_d_invk_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_invk_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_invt_solver.f90

@@ -94,7 +94,7 @@ module amg_d_invt_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_invt_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_jac_smoother.f90

@@ -151,7 +151,7 @@ module amg_d_jac_smoother
       import :: psb_desc_type, amg_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, &
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -274,7 +274,7 @@ module amg_d_jac_smoother
       import :: psb_desc_type, amg_d_l1_jac_smoother_type, psb_d_vect_type, &
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_l1_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -345,6 +345,7 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
+      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_d_jac_solver.f90

@@ -143,7 +143,7 @@ module amg_d_jac_solver
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_d_jac_solver
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_l1_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_krm_solver.f90

@@ -174,7 +174,7 @@ module amg_d_krm_solver
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_krm_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_d_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_mumps_solver_type), intent(inout)       :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_onelev_mod.f90

@@ -155,13 +155,35 @@ module amg_d_onelev_mod
   private :: d_wrk_alloc, d_wrk_free, &
        & d_wrk_clone, d_wrk_move_alloc, d_wrk_cnv, d_wrk_sizeof
 
+  !
+  ! Remap.
+  !  This keeps track of remapping.
+  !  Logic here is as follows:
+  !  1. AC_PRE_REMAP  need to figure out if we
+  !                really need it.
+  !  2. DESC_AC_PRE_REMAP  contains the descriptor before
+  !                remapping. Meaning that it is possible
+  !                to implement the RESTRICTOR operator by
+  !                 a. Doing LINMAP_U2V onto this one
+  !                 b. For each process, send the data to
+  !                    IDEST.
+  !                This assumes that remapping goes by
+  !                a factor of 2.
+  !                For the PROLONGATOR operators, we first
+  !                use DESC_AC, then split and send onto
+  !                the processes in DESC_AC_PRE_REMAP.
+  !
+  ! To be fixed: what happens if NP the starting processes
+  ! is not an even number? Coordinate with _X_remap in PSBLAS
+  ! 
   type amg_d_remap_data_type
     type(psb_dspmat_type)          :: ac_pre_remap
     type(psb_desc_type)            :: desc_ac_pre_remap
     integer(psb_ipk_)              :: idest
     integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
   contains
-    procedure, pass(rmp) :: clone   => d_remap_data_clone
+    procedure, pass(rmp) :: clone      => d_remap_data_clone
+    procedure, pass(rmp) :: move_alloc => d_remap_move_alloc
   end type amg_d_remap_data_type
 
   type amg_d_onelev_type
@@ -207,7 +229,7 @@ module amg_d_onelev_mod
     procedure, pass(lv) :: get_wrksz => d_base_onelev_get_wrksize
     procedure, pass(lv) :: allocate_wrk   => d_base_onelev_allocate_wrk
     procedure, pass(lv) :: free_wrk       => d_base_onelev_free_wrk
-    procedure, nopass   :: stringval => amg_stringval
+    procedure, nopass   :: stringval  => amg_stringval
     procedure, pass(lv) :: move_alloc => d_base_onelev_move_alloc
 
 
@@ -619,7 +641,7 @@ contains
     ! Arguments
     class(amg_d_onelev_type), target, intent(inout) :: lv
     class(amg_d_onelev_type), target, intent(inout) :: lvout
-    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_ipk_), intent(out)                  :: info
 
     info = psb_success_
     if (allocated(lv%sm)) then
@@ -685,6 +707,7 @@ contains
     if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info)
     if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info)
     if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info)
+    if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info)
     b%base_a    => lv%base_a
     b%base_desc => lv%base_desc
 
@@ -739,6 +762,7 @@ contains
     info = psb_success_
     nwv = lv%get_wrksz()
     if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
+!!$    write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb()
     if (info == 0) then
       if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
         !
@@ -786,47 +810,32 @@ contains
 
     info = psb_success_
     call wk%free(info)
-    if  (present(desc2)) then
-!!$      write(0,*) 'Check on wrk_alloc 2',&
-!!$           & desc2%get_local_rows(), desc%get_local_rows(),&
-!!$           & desc2%get_local_cols(),desc%get_local_cols()
-!!$      flush(0)
+!!$    write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',&
+!!$         & present(desc2),desc%is_valid()
+
+    allocate(wk%wv(nwv),stat=info)    
+    if  (present(desc2).and.(desc%is_valid())) then
+!!$      write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
       if (desc2%get_local_cols()>desc%get_local_cols()) then
-        call psb_geasb(wk%vx2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc2,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
       else
-!!$        write(0,*) 'Check on wrk_alloc 1.5 ',&
-!!$             & desc%get_local_rows(),&
-!!$             & desc%get_local_cols()
-        call psb_geasb(wk%vx2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)
       end if
-    else
-!!$      write(0,*) 'Check on wrk_alloc 1 ',&
-!!$           & desc%get_local_rows(),&
-!!$           & desc%get_local_cols()
+    else if (present(desc2)) then
+      call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
+    else if (desc%is_valid()) then
+      call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)      
+    end if
+
+  contains
+    subroutine inner_do_wrk_alloc(wk,nwv,desc,vmold)
+      class(amg_dmlprec_wrk_type), target, intent(inout) :: wk
+      integer(psb_ipk_), intent(in)                     :: nwv
+      type(psb_desc_type), intent(in)                   :: desc
+      class(psb_d_base_vect_type), intent(in), optional  :: vmold
+      
+      integer(psb_ipk_) :: i
+      
       call psb_geasb(wk%vx2l,desc,info,&
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vy2l,desc,info,&
@@ -835,12 +844,11 @@ contains
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vty,desc,info,&
            & scratch=.true.,mold=vmold)
-      allocate(wk%wv(nwv),stat=info)
       do i=1,nwv
         call psb_geasb(wk%wv(i),desc,info,&
              & scratch=.true.,mold=vmold)
       end do
-    end if
+    end subroutine inner_do_wrk_alloc
   end subroutine d_wrk_alloc
 
   subroutine d_wrk_free(wk,info)
@@ -993,4 +1001,25 @@ contains
     call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
   end subroutine d_remap_data_clone
 
+  subroutine d_remap_move_alloc(rmp, remap_out, info)
+    use psb_base_mod
+    implicit none
+    ! Arguments
+    class(amg_d_remap_data_type), target, intent(inout) :: rmp
+    class(amg_d_remap_data_type), target, intent(inout) :: remap_out
+    integer(psb_ipk_), intent(out)                    :: info
+    !
+    integer(psb_ipk_) :: i
+
+    info = psb_success_
+
+    call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info)
+    if (info == psb_success_) &
+         & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info)
+    remap_out%idest = rmp%idest
+    call move_alloc(rmp%isrc,remap_out%isrc)
+    call move_alloc(rmp%nrsrc,remap_out%nrsrc)
+    call move_alloc(rmp%naggr,remap_out%naggr)
+  end subroutine d_remap_move_alloc
+  
 end module amg_d_onelev_mod

amgprec/amg_d_parmatch_aggregator_mod.F90

@@ -136,7 +136,7 @@ module amg_d_parmatch_aggregator_mod
     procedure, pass(ag) :: mat_bld      => amg_d_parmatch_aggregator_mat_bld
     procedure, pass(ag) :: mat_asb      => amg_d_parmatch_aggregator_mat_asb
     procedure, pass(ag) :: inner_mat_asb      => amg_d_parmatch_aggregator_inner_mat_asb
-    procedure, pass(ag) :: bld_map      => amg_d_parmatch_aggregator_bld_map
+    procedure, pass(ag) :: bld_linmap   => amg_d_parmatch_aggregator_bld_linmap
     procedure, pass(ag) :: csetc        => amg_d_parmatch_aggr_csetc
     procedure, pass(ag) :: cseti        => amg_d_parmatch_aggr_cseti
     procedure, pass(ag) :: default      => amg_d_parmatch_aggr_set_default
@@ -643,7 +643,7 @@ contains
     end select
   end subroutine amg_d_parmatch_aggregator_clone
 
-  subroutine  amg_d_parmatch_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
+  subroutine  amg_d_parmatch_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
        & op_restr,op_prol,map,info)
     use psb_base_mod
     implicit none
@@ -654,7 +654,7 @@ contains
     type(psb_dlinmap_type), intent(out)     :: map
     integer(psb_ipk_), intent(out)          :: info
     integer(psb_ipk_) :: err_act
-    character(len=20) :: name='d_parmatch_aggregator_bld_map'
+    character(len=20) :: name='d_parmatch_aggregator_bld_linmap'
 
     info = psb_success_
     call psb_erractionsave(err_act)
@@ -680,5 +680,5 @@ contains
 9999 call psb_error_handler(err_act)
 
     return
-  end subroutine amg_d_parmatch_aggregator_bld_map
+  end subroutine amg_d_parmatch_aggregator_bld_linmap
 end module amg_d_parmatch_aggregator_mod

amgprec/amg_d_poly_smoother.f90

@@ -140,7 +140,7 @@ module amg_d_poly_smoother
       import :: psb_desc_type, amg_d_poly_smoother_type, psb_d_vect_type, psb_dpk_, &
            & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_d_poly_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -279,6 +279,7 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
+      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_d_prec_type.f90

@@ -102,6 +102,7 @@ module amg_d_prec_type
     ! The multilevel hierarchy
     !
     type(amg_d_onelev_type), allocatable :: precv(:)
+    integer(psb_ipk_)                    :: nlevs
   contains
     procedure, pass(prec)               :: psb_d_apply2_vect => amg_d_apply2_vect
     procedure, pass(prec)               :: psb_d_apply1_vect => amg_d_apply1_vect
@@ -119,6 +120,7 @@ module amg_d_prec_type
     procedure, pass(prec)               :: cmp_complexity => amg_d_cmp_compl
     procedure, pass(prec)               :: get_avg_cr => amg_d_get_avg_cr
     procedure, pass(prec)               :: cmp_avg_cr => amg_d_cmp_avg_cr
+    procedure, pass(prec)               :: set_nlevs  => amg_d_set_nlevs
     procedure, pass(prec)               :: get_nlevs  => amg_d_get_nlevs
     procedure, pass(prec)               :: get_nzeros => amg_d_get_nzeros
     procedure, pass(prec)               :: sizeof => amg_dprec_sizeof
@@ -311,7 +313,7 @@ module amg_d_prec_type
            & psb_d_base_sparse_mat, psb_d_base_vect_type, &
            & psb_i_base_vect_type, amg_dprec_type, psb_ipk_
       implicit none
-      type(psb_dspmat_type), intent(in), target          :: a
+      type(psb_dspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_dprec_type), intent(inout), target       :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -323,14 +325,15 @@ module amg_d_prec_type
   end interface amg_precbld
 
   interface amg_hierarchy_bld
-    subroutine amg_d_hierarchy_bld(a,desc_a,prec,info)
+    subroutine amg_d_hierarchy_bld(a,desc_a,prec,info,cpymat)
       import :: psb_dspmat_type, psb_desc_type, psb_dpk_, &
            & amg_dprec_type, psb_ipk_
       implicit none
-      type(psb_dspmat_type), intent(in), target          :: a
+      type(psb_dspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_dprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
+      logical, intent(in), optional                        :: cpymat
       !      character, intent(in),optional             :: upd
     end subroutine amg_d_hierarchy_bld
   end interface amg_hierarchy_bld
@@ -434,10 +437,19 @@ contains
     class(amg_dprec_type), intent(in) :: prec
     integer(psb_ipk_) :: val
     val = 0
-    if (allocated(prec%precv)) then
-      val = size(prec%precv)
-    end if
+!!$    if (allocated(prec%precv)) then
+!!$      val = size(prec%precv)
+!!$    end if
+    val = prec%nlevs
+!!$    write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
   end function amg_d_get_nlevs
+
+  subroutine amg_d_set_nlevs(prec,nl) 
+    implicit none
+    class(amg_dprec_type), intent(inout) :: prec
+    integer(psb_ipk_) :: nl
+    prec%nlevs = nl
+  end subroutine amg_d_set_nlevs
   !
   ! Function returning the size of the amg_prec_type data structure
   ! in bytes or in number of nonzeros of the operator(s) involved.
@@ -508,7 +520,7 @@ contains
 
     real(psb_dpk_) :: num, den, nmin
     type(psb_ctxt_type) :: ctxt
-    integer(psb_ipk_)   :: il
+    integer(psb_ipk_)   :: il,nl
 
     num = -done
     den = done
@@ -518,7 +530,10 @@ contains
       num = prec%precv(il)%base_a%get_nzeros()
       if (num >= dzero) then
         den = num
-        do il=2,size(prec%precv)
+        nl = prec%get_nlevs()
+!!$        write(0,*) 'Inside cmp_compl ',nl,size(prec%precv)
+        do il=2, nl
+!!$          write(0,*) '                 ',il,associated(prec%precv(il)%base_a)
           num = num + max(0,prec%precv(il)%base_a%get_nzeros())
         end do
       end if
@@ -549,7 +564,6 @@ contains
   end function amg_d_get_avg_cr
 
   subroutine amg_d_cmp_avg_cr(prec)
-
     implicit none
     class(amg_dprec_type), intent(inout) :: prec
 
@@ -557,17 +571,18 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_)   :: il, nl, iam, np
 
-
     avgcr = dzero
+    nl = prec%get_nlevs()
+    do il=2,nl
+      if (prec%precv(il)%base_desc%is_ok()) then 
+        ctxt = prec%precv(il)%base_desc%get_ctxt()
+        call psb_info(ctxt,iam,np)
+        if (iam >=0)  avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
+      end if
+    end do
+    avgcr = avgcr / (nl-1)
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    if (allocated(prec%precv)) then
-      nl = size(prec%precv)
-      do il=2,nl
-        avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
-      end do
-      avgcr = avgcr / (nl-1)
-    end if
     call psb_sum(ctxt,avgcr)
     prec%ag_data%avg_cr = avgcr/np
   end subroutine amg_d_cmp_avg_cr
@@ -585,9 +600,7 @@ contains
   !             error code.
   !
   subroutine amg_dprecfree(p,info)
-
     implicit none
-
     ! Arguments
     type(amg_dprec_type), intent(inout) :: p
     integer(psb_ipk_), intent(out)        :: info
@@ -613,9 +626,7 @@ contains
   end subroutine amg_dprecfree
 
   subroutine amg_d_prec_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_dprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -647,9 +658,7 @@ contains
   end subroutine amg_d_prec_free
 
   subroutine amg_d_smoothers_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_dprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -665,12 +674,8 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
-    !
-    ! In the internals, do FREE on components,
-    ! but do not deallocate them
-    !
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         call prec%precv(i)%free_smoothers(info)
       end do
     end if
@@ -684,9 +689,7 @@ contains
   end subroutine amg_d_smoothers_free
 
   subroutine amg_d_hierarchy_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_dprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -713,7 +716,6 @@ contains
 
   end subroutine amg_d_hierarchy_free
 
-
   !
   ! Top level methods.
   !
@@ -778,7 +780,6 @@ contains
 
   end subroutine amg_d_apply1_vect
 
-
   subroutine amg_d_apply2v(prec,x,y,desc_data,info,trans,work)
     implicit none
     type(psb_desc_type),intent(in)    :: desc_data
@@ -843,7 +844,6 @@ contains
   subroutine amg_d_dump(prec,info,istart,iend,iproc,prefix,head,&
        & ac,rp,smoother,solver,tprol,&
        & global_num)
-
     implicit none
     class(amg_dprec_type), intent(in)     :: prec
     integer(psb_ipk_), intent(out)          :: info
@@ -860,7 +860,7 @@ contains
     info = 0
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    iln = size(prec%precv)
+    iln = prec%get_nlevs()
     if (present(istart)) then
       il1 = max(1,istart)
     else
@@ -884,7 +884,6 @@ contains
   end subroutine amg_d_dump
 
   subroutine amg_d_cnv(prec,info,amold,vmold,imold)
-
     implicit none
     class(amg_dprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)       :: info
@@ -896,7 +895,7 @@ contains
 
     info = psb_success_
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         if (info == psb_success_ ) &
              & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
       end do
@@ -905,7 +904,6 @@ contains
   end subroutine amg_d_cnv
 
   subroutine amg_d_clone(prec,precout,info)
-
     implicit none
     class(amg_dprec_type), intent(inout) :: prec
     class(psb_dprec_type), intent(inout) :: precout
@@ -917,7 +915,6 @@ contains
   end subroutine amg_d_clone
 
   subroutine amg_d_inner_clone(prec,precout,info)
-
     implicit none
     class(amg_dprec_type), intent(inout)         :: prec
     class(psb_dprec_type), target, intent(inout) :: precout
@@ -933,8 +930,9 @@ contains
       pout%ctxt          = prec%ctxt
       pout%ag_data       = prec%ag_data
       pout%outer_sweeps  = prec%outer_sweeps
+      pout%nlevs         = prec%nlevs
       if (allocated(prec%precv)) then
-        ln = size(prec%precv)
+        ln = prec%get_nlevs()
         allocate(pout%precv(ln),stat=info)
         if (info /= psb_success_) goto 9999
         if (ln >= 1) then
@@ -942,6 +940,7 @@ contains
         end if
         do lev=2, ln
           if (info /= psb_success_) exit
+!!$          write(0,*) 'Inner_clone must be checked and reimplemented! '
           call prec%precv(lev)%clone(pout%precv(lev),info)
           if (info == psb_success_) then
             pout%precv(lev)%base_a       => pout%precv(lev)%ac
@@ -1006,7 +1005,7 @@ contains
     integer(psb_ipk_), intent(out)        :: info
     class(psb_d_base_vect_type), intent(in), optional  :: vmold
     !
-    ! In MLD the DESC optional argument is ignored, since
+    ! In AMG the DESC optional argument is ignored, since
     ! the necessary info is contained in the various entries of the
     ! PRECV component.
     type(psb_desc_type), intent(in), optional  :: desc
@@ -1021,7 +1020,7 @@ contains
     if (psb_errstatus_fatal()) then
       info = psb_err_internal_error_; goto 9999
     end if
-    nlev   = size(prec%precv)
+    nlev   = prec%get_nlevs()
     level = 1
     do level = 1, nlev
       call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@@ -1044,7 +1043,6 @@ contains
   subroutine amg_d_free_wrk(prec,info)
     use psb_base_mod
     implicit none
-
     ! Arguments
     class(amg_dprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -1061,7 +1059,7 @@ contains
     end if
 
     if (allocated(prec%precv)) then
-      nlev   = size(prec%precv)
+      nlev   = prec%get_nlevs()
       do level = 1, nlev
         call prec%precv(level)%free_wrk(info)
       end do

amgprec/amg_d_slu_solver.F90

@@ -124,7 +124,7 @@ module amg_d_slu_solver
       Implicit None
 
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_slu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_d_sludist_solver.F90

@@ -259,7 +259,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_dspmat_type), intent(in), target           :: a
+    type(psb_dspmat_type), intent(inout), target        :: a
     Type(psb_desc_type), Intent(inout)                  :: desc_a 
     class(amg_d_sludist_solver_type), intent(inout)     :: sv
     integer, intent(out)                                :: info

amgprec/amg_d_umf_solver.F90

@@ -163,7 +163,7 @@ module amg_d_umf_solver
       Implicit None
       
       ! Arguments
-      type(psb_dspmat_type), intent(in), target           :: a
+      type(psb_dspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_d_umf_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_s_ainv_solver.F90

@@ -103,7 +103,7 @@ module amg_s_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_ainv_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -225,7 +225,7 @@ module amg_s_ainv_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
            & psb_d_vect_type, psb_s_base_vect_type, psb_spk_, psb_ipk_
       implicit none
-      type(psb_sspmat_type), intent(in), target   :: a
+      type(psb_sspmat_type), intent(inout), target  :: a
       integer(psb_ipk_), intent(in)               :: fillin,alg
       real(psb_spk_), intent(in)                  :: thresh
       type(psb_sspmat_type), intent(inout)        :: wmat, zmat

amgprec/amg_s_as_smoother.f90

@@ -230,7 +230,7 @@ module amg_s_as_smoother
            & psb_desc_type, psb_s_base_sparse_mat, psb_ipk_,&
            & psb_i_base_vect_type
       implicit none 
-      type(psb_sspmat_type), intent(in), target        :: a
+      type(psb_sspmat_type), intent(inout), target     :: a
       Type(psb_desc_type), Intent(inout)                 :: desc_a 
       class(amg_s_as_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                     :: info

amgprec/amg_s_base_aggregator_mod.f90

@@ -89,7 +89,7 @@ module amg_s_base_aggregator_mod
     procedure, pass(ag) :: bld_tprol   => amg_s_base_aggregator_build_tprol
     procedure, pass(ag) :: mat_bld     => amg_s_base_aggregator_mat_bld
     procedure, pass(ag) :: mat_asb     => amg_s_base_aggregator_mat_asb
-    procedure, pass(ag) :: bld_map     => amg_s_base_aggregator_bld_map
+    procedure, pass(ag) :: bld_linmap  => amg_s_base_aggregator_bld_linmap
     procedure, pass(ag) :: update_next => amg_s_base_aggregator_update_next
     procedure, pass(ag) :: clone       => amg_s_base_aggregator_clone
     procedure, pass(ag) :: free        => amg_s_base_aggregator_free
@@ -458,7 +458,7 @@ contains
   end subroutine amg_s_base_aggregator_mat_asb
 
   !
-  !> Function   bld_map
+  !> Function   bld_linmap
   !! \memberof  amg_s_base_aggregator_type
   !! \brief     Build linear map between hierarchy levels 
   !!
@@ -473,7 +473,7 @@ contains
   !!  \param map       The output map
   !!  \param info    Return code
   !!  
-  subroutine  amg_s_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
+  subroutine  amg_s_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
        & op_restr,op_prol,map,info)
     use psb_base_mod
     implicit none
@@ -484,7 +484,7 @@ contains
     type(psb_slinmap_type), intent(out)    :: map
     integer(psb_ipk_), intent(out)       :: info
     integer(psb_ipk_) :: err_act
-    character(len=20) :: name='s_base_aggregator_bld_map'
+    character(len=20) :: name='s_base_aggregator_bld_linmap'
 
     info = psb_success_
     call psb_erractionsave(err_act)
@@ -508,7 +508,6 @@ contains
 9999 call psb_error_handler(err_act)
 
     return
-  end subroutine amg_s_base_aggregator_bld_map
-
+  end subroutine amg_s_base_aggregator_bld_linmap
 
 end module amg_s_base_aggregator_mod

amgprec/amg_s_base_smoother_mod.f90

@@ -237,7 +237,7 @@ module amg_s_base_smoother_mod
            & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
            & amg_s_base_smoother_type, psb_ipk_, psb_i_base_vect_type
       ! Arguments
-      type(psb_sspmat_type), intent(in), target     :: a
+      type(psb_sspmat_type), intent(inout), target  :: a
       Type(psb_desc_type), Intent(inout)              :: desc_a 
       class(amg_s_base_smoother_type), intent(inout) :: sm 
       integer(psb_ipk_), intent(out)                  :: info

amgprec/amg_s_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_s_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_sspmat_type), intent(in), target             :: a
+      type(psb_sspmat_type), intent(inout), target          :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_s_base_solver_type), intent(inout)          :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_s_diag_solver.f90

@@ -119,7 +119,7 @@ module amg_s_diag_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
            & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
            & amg_s_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_s_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_s_l1_diag_solver
       import :: psb_desc_type, psb_sspmat_type,  psb_s_base_sparse_mat, &
            & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
            & amg_s_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_s_l1_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_s_gs_solver.f90

@@ -181,7 +181,7 @@ module amg_s_gs_solver
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_gs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_s_gs_solver
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_bwgs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_sspmat_type), intent(in), target           :: a
+    type(psb_sspmat_type), intent(inout), target        :: a
     Type(psb_desc_type), Intent(inout)                  :: desc_a 
     class(amg_s_id_solver_type), intent(inout)          :: sv
     integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_ilu_solver.f90

@@ -144,7 +144,7 @@ module amg_s_ilu_solver
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_ilu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_inner_mod.f90

@@ -56,7 +56,7 @@ module amg_s_inner_mod
            & psb_spk_, psb_s_base_sparse_mat, psb_s_base_vect_type, psb_ipk_
       import :: amg_sprec_type
       implicit none
-      type(psb_sspmat_type), intent(in), target          :: a
+      type(psb_sspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       type(amg_sprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -67,7 +67,7 @@ module amg_s_inner_mod
   end interface amg_mlprec_bld
 
   interface amg_mlprec_aply
-    subroutine amg_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
+    subroutine amg_smlprec_aply_a(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_
       import :: amg_sprec_type
       implicit none 
@@ -79,7 +79,7 @@ module amg_s_inner_mod
       character,intent(in)               :: trans
       real(psb_spk_),target             :: work(:)
       integer(psb_ipk_), intent(out)     :: info
-    end subroutine amg_smlprec_aply
+    end subroutine amg_smlprec_aply_a
     subroutine amg_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_sspmat_type, psb_desc_type, &
            & psb_spk_, psb_s_vect_type, psb_ipk_

amgprec/amg_s_invk_solver.f90

@@ -94,7 +94,7 @@ module amg_s_invk_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_invk_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_invt_solver.f90

@@ -94,7 +94,7 @@ module amg_s_invt_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_invt_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_jac_smoother.f90

@@ -151,7 +151,7 @@ module amg_s_jac_smoother
       import :: psb_desc_type, amg_s_jac_smoother_type, psb_s_vect_type, psb_spk_, &
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -274,7 +274,7 @@ module amg_s_jac_smoother
       import :: psb_desc_type, amg_s_l1_jac_smoother_type, psb_s_vect_type, &
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_l1_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -345,6 +345,7 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
+      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_s_jac_solver.f90

@@ -143,7 +143,7 @@ module amg_s_jac_solver
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_s_jac_solver
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_l1_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_krm_solver.f90

@@ -174,7 +174,7 @@ module amg_s_krm_solver
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_krm_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_s_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_mumps_solver_type), intent(inout)       :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_s_onelev_mod.f90

@@ -155,13 +155,35 @@ module amg_s_onelev_mod
   private :: s_wrk_alloc, s_wrk_free, &
        & s_wrk_clone, s_wrk_move_alloc, s_wrk_cnv, s_wrk_sizeof
 
+  !
+  ! Remap.
+  !  This keeps track of remapping.
+  !  Logic here is as follows:
+  !  1. AC_PRE_REMAP  need to figure out if we
+  !                really need it.
+  !  2. DESC_AC_PRE_REMAP  contains the descriptor before
+  !                remapping. Meaning that it is possible
+  !                to implement the RESTRICTOR operator by
+  !                 a. Doing LINMAP_U2V onto this one
+  !                 b. For each process, send the data to
+  !                    IDEST.
+  !                This assumes that remapping goes by
+  !                a factor of 2.
+  !                For the PROLONGATOR operators, we first
+  !                use DESC_AC, then split and send onto
+  !                the processes in DESC_AC_PRE_REMAP.
+  !
+  ! To be fixed: what happens if NP the starting processes
+  ! is not an even number? Coordinate with _X_remap in PSBLAS
+  ! 
   type amg_s_remap_data_type
     type(psb_sspmat_type)          :: ac_pre_remap
     type(psb_desc_type)            :: desc_ac_pre_remap
     integer(psb_ipk_)              :: idest
     integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
   contains
-    procedure, pass(rmp) :: clone   => s_remap_data_clone
+    procedure, pass(rmp) :: clone      => s_remap_data_clone
+    procedure, pass(rmp) :: move_alloc => s_remap_move_alloc
   end type amg_s_remap_data_type
 
   type amg_s_onelev_type
@@ -207,7 +229,7 @@ module amg_s_onelev_mod
     procedure, pass(lv) :: get_wrksz => s_base_onelev_get_wrksize
     procedure, pass(lv) :: allocate_wrk   => s_base_onelev_allocate_wrk
     procedure, pass(lv) :: free_wrk       => s_base_onelev_free_wrk
-    procedure, nopass   :: stringval => amg_stringval
+    procedure, nopass   :: stringval  => amg_stringval
     procedure, pass(lv) :: move_alloc => s_base_onelev_move_alloc
 
 
@@ -619,7 +641,7 @@ contains
     ! Arguments
     class(amg_s_onelev_type), target, intent(inout) :: lv
     class(amg_s_onelev_type), target, intent(inout) :: lvout
-    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_ipk_), intent(out)                  :: info
 
     info = psb_success_
     if (allocated(lv%sm)) then
@@ -685,6 +707,7 @@ contains
     if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info)
     if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info)
     if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info)
+    if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info)
     b%base_a    => lv%base_a
     b%base_desc => lv%base_desc
 
@@ -739,6 +762,7 @@ contains
     info = psb_success_
     nwv = lv%get_wrksz()
     if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
+!!$    write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb()
     if (info == 0) then
       if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
         !
@@ -786,47 +810,32 @@ contains
 
     info = psb_success_
     call wk%free(info)
-    if  (present(desc2)) then
-!!$      write(0,*) 'Check on wrk_alloc 2',&
-!!$           & desc2%get_local_rows(), desc%get_local_rows(),&
-!!$           & desc2%get_local_cols(),desc%get_local_cols()
-!!$      flush(0)
+!!$    write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',&
+!!$         & present(desc2),desc%is_valid()
+
+    allocate(wk%wv(nwv),stat=info)    
+    if  (present(desc2).and.(desc%is_valid())) then
+!!$      write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
       if (desc2%get_local_cols()>desc%get_local_cols()) then
-        call psb_geasb(wk%vx2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc2,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
       else
-!!$        write(0,*) 'Check on wrk_alloc 1.5 ',&
-!!$             & desc%get_local_rows(),&
-!!$             & desc%get_local_cols()
-        call psb_geasb(wk%vx2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)
       end if
-    else
-!!$      write(0,*) 'Check on wrk_alloc 1 ',&
-!!$           & desc%get_local_rows(),&
-!!$           & desc%get_local_cols()
+    else if (present(desc2)) then
+      call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
+    else if (desc%is_valid()) then
+      call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)      
+    end if
+
+  contains
+    subroutine inner_do_wrk_alloc(wk,nwv,desc,vmold)
+      class(amg_smlprec_wrk_type), target, intent(inout) :: wk
+      integer(psb_ipk_), intent(in)                     :: nwv
+      type(psb_desc_type), intent(in)                   :: desc
+      class(psb_s_base_vect_type), intent(in), optional  :: vmold
+      
+      integer(psb_ipk_) :: i
+      
       call psb_geasb(wk%vx2l,desc,info,&
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vy2l,desc,info,&
@@ -835,12 +844,11 @@ contains
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vty,desc,info,&
            & scratch=.true.,mold=vmold)
-      allocate(wk%wv(nwv),stat=info)
       do i=1,nwv
         call psb_geasb(wk%wv(i),desc,info,&
              & scratch=.true.,mold=vmold)
       end do
-    end if
+    end subroutine inner_do_wrk_alloc
   end subroutine s_wrk_alloc
 
   subroutine s_wrk_free(wk,info)
@@ -993,4 +1001,25 @@ contains
     call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
   end subroutine s_remap_data_clone
 
+  subroutine s_remap_move_alloc(rmp, remap_out, info)
+    use psb_base_mod
+    implicit none
+    ! Arguments
+    class(amg_s_remap_data_type), target, intent(inout) :: rmp
+    class(amg_s_remap_data_type), target, intent(inout) :: remap_out
+    integer(psb_ipk_), intent(out)                    :: info
+    !
+    integer(psb_ipk_) :: i
+
+    info = psb_success_
+
+    call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info)
+    if (info == psb_success_) &
+         & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info)
+    remap_out%idest = rmp%idest
+    call move_alloc(rmp%isrc,remap_out%isrc)
+    call move_alloc(rmp%nrsrc,remap_out%nrsrc)
+    call move_alloc(rmp%naggr,remap_out%naggr)
+  end subroutine s_remap_move_alloc
+  
 end module amg_s_onelev_mod

amgprec/amg_s_parmatch_aggregator_mod.F90

@@ -136,7 +136,7 @@ module amg_s_parmatch_aggregator_mod
     procedure, pass(ag) :: mat_bld      => amg_s_parmatch_aggregator_mat_bld
     procedure, pass(ag) :: mat_asb      => amg_s_parmatch_aggregator_mat_asb
     procedure, pass(ag) :: inner_mat_asb      => amg_s_parmatch_aggregator_inner_mat_asb
-    procedure, pass(ag) :: bld_map      => amg_s_parmatch_aggregator_bld_map
+    procedure, pass(ag) :: bld_linmap   => amg_s_parmatch_aggregator_bld_linmap
     procedure, pass(ag) :: csetc        => amg_s_parmatch_aggr_csetc
     procedure, pass(ag) :: cseti        => amg_s_parmatch_aggr_cseti
     procedure, pass(ag) :: default      => amg_s_parmatch_aggr_set_default
@@ -643,7 +643,7 @@ contains
     end select
   end subroutine amg_s_parmatch_aggregator_clone
 
-  subroutine  amg_s_parmatch_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
+  subroutine  amg_s_parmatch_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
        & op_restr,op_prol,map,info)
     use psb_base_mod
     implicit none
@@ -654,7 +654,7 @@ contains
     type(psb_slinmap_type), intent(out)     :: map
     integer(psb_ipk_), intent(out)          :: info
     integer(psb_ipk_) :: err_act
-    character(len=20) :: name='s_parmatch_aggregator_bld_map'
+    character(len=20) :: name='s_parmatch_aggregator_bld_linmap'
 
     info = psb_success_
     call psb_erractionsave(err_act)
@@ -680,5 +680,5 @@ contains
 9999 call psb_error_handler(err_act)
 
     return
-  end subroutine amg_s_parmatch_aggregator_bld_map
+  end subroutine amg_s_parmatch_aggregator_bld_linmap
 end module amg_s_parmatch_aggregator_mod

amgprec/amg_s_poly_smoother.f90

@@ -140,7 +140,7 @@ module amg_s_poly_smoother
       import :: psb_desc_type, amg_s_poly_smoother_type, psb_s_vect_type, psb_spk_, &
            & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_s_poly_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -279,6 +279,7 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
+      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_s_prec_type.f90

@@ -102,6 +102,7 @@ module amg_s_prec_type
     ! The multilevel hierarchy
     !
     type(amg_s_onelev_type), allocatable :: precv(:)
+    integer(psb_ipk_)                    :: nlevs
   contains
     procedure, pass(prec)               :: psb_s_apply2_vect => amg_s_apply2_vect
     procedure, pass(prec)               :: psb_s_apply1_vect => amg_s_apply1_vect
@@ -119,6 +120,7 @@ module amg_s_prec_type
     procedure, pass(prec)               :: cmp_complexity => amg_s_cmp_compl
     procedure, pass(prec)               :: get_avg_cr => amg_s_get_avg_cr
     procedure, pass(prec)               :: cmp_avg_cr => amg_s_cmp_avg_cr
+    procedure, pass(prec)               :: set_nlevs  => amg_s_set_nlevs
     procedure, pass(prec)               :: get_nlevs  => amg_s_get_nlevs
     procedure, pass(prec)               :: get_nzeros => amg_s_get_nzeros
     procedure, pass(prec)               :: sizeof => amg_sprec_sizeof
@@ -311,7 +313,7 @@ module amg_s_prec_type
            & psb_s_base_sparse_mat, psb_s_base_vect_type, &
            & psb_i_base_vect_type, amg_sprec_type, psb_ipk_
       implicit none
-      type(psb_sspmat_type), intent(in), target          :: a
+      type(psb_sspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_sprec_type), intent(inout), target       :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -323,14 +325,15 @@ module amg_s_prec_type
   end interface amg_precbld
 
   interface amg_hierarchy_bld
-    subroutine amg_s_hierarchy_bld(a,desc_a,prec,info)
+    subroutine amg_s_hierarchy_bld(a,desc_a,prec,info,cpymat)
       import :: psb_sspmat_type, psb_desc_type, psb_spk_, &
            & amg_sprec_type, psb_ipk_
       implicit none
-      type(psb_sspmat_type), intent(in), target          :: a
+      type(psb_sspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_sprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
+      logical, intent(in), optional                        :: cpymat
       !      character, intent(in),optional             :: upd
     end subroutine amg_s_hierarchy_bld
   end interface amg_hierarchy_bld
@@ -434,10 +437,19 @@ contains
     class(amg_sprec_type), intent(in) :: prec
     integer(psb_ipk_) :: val
     val = 0
-    if (allocated(prec%precv)) then
-      val = size(prec%precv)
-    end if
+!!$    if (allocated(prec%precv)) then
+!!$      val = size(prec%precv)
+!!$    end if
+    val = prec%nlevs
+!!$    write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
   end function amg_s_get_nlevs
+
+  subroutine amg_s_set_nlevs(prec,nl) 
+    implicit none
+    class(amg_sprec_type), intent(inout) :: prec
+    integer(psb_ipk_) :: nl
+    prec%nlevs = nl
+  end subroutine amg_s_set_nlevs
   !
   ! Function returning the size of the amg_prec_type data structure
   ! in bytes or in number of nonzeros of the operator(s) involved.
@@ -508,7 +520,7 @@ contains
 
     real(psb_spk_) :: num, den, nmin
     type(psb_ctxt_type) :: ctxt
-    integer(psb_ipk_)   :: il
+    integer(psb_ipk_)   :: il,nl
 
     num = -sone
     den = sone
@@ -518,7 +530,10 @@ contains
       num = prec%precv(il)%base_a%get_nzeros()
       if (num >= szero) then
         den = num
-        do il=2,size(prec%precv)
+        nl = prec%get_nlevs()
+!!$        write(0,*) 'Inside cmp_compl ',nl,size(prec%precv)
+        do il=2, nl
+!!$          write(0,*) '                 ',il,associated(prec%precv(il)%base_a)
           num = num + max(0,prec%precv(il)%base_a%get_nzeros())
         end do
       end if
@@ -549,7 +564,6 @@ contains
   end function amg_s_get_avg_cr
 
   subroutine amg_s_cmp_avg_cr(prec)
-
     implicit none
     class(amg_sprec_type), intent(inout) :: prec
 
@@ -557,17 +571,18 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_)   :: il, nl, iam, np
 
-
     avgcr = szero
+    nl = prec%get_nlevs()
+    do il=2,nl
+      if (prec%precv(il)%base_desc%is_ok()) then 
+        ctxt = prec%precv(il)%base_desc%get_ctxt()
+        call psb_info(ctxt,iam,np)
+        if (iam >=0)  avgcr = avgcr + max(szero,prec%precv(il)%szratio)
+      end if
+    end do
+    avgcr = avgcr / (nl-1)
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    if (allocated(prec%precv)) then
-      nl = size(prec%precv)
-      do il=2,nl
-        avgcr = avgcr + max(szero,prec%precv(il)%szratio)
-      end do
-      avgcr = avgcr / (nl-1)
-    end if
     call psb_sum(ctxt,avgcr)
     prec%ag_data%avg_cr = avgcr/np
   end subroutine amg_s_cmp_avg_cr
@@ -585,9 +600,7 @@ contains
   !             error code.
   !
   subroutine amg_sprecfree(p,info)
-
     implicit none
-
     ! Arguments
     type(amg_sprec_type), intent(inout) :: p
     integer(psb_ipk_), intent(out)        :: info
@@ -613,9 +626,7 @@ contains
   end subroutine amg_sprecfree
 
   subroutine amg_s_prec_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_sprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -647,9 +658,7 @@ contains
   end subroutine amg_s_prec_free
 
   subroutine amg_s_smoothers_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_sprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -665,12 +674,8 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
-    !
-    ! In the internals, do FREE on components,
-    ! but do not deallocate them
-    !
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         call prec%precv(i)%free_smoothers(info)
       end do
     end if
@@ -684,9 +689,7 @@ contains
   end subroutine amg_s_smoothers_free
 
   subroutine amg_s_hierarchy_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_sprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -713,7 +716,6 @@ contains
 
   end subroutine amg_s_hierarchy_free
 
-
   !
   ! Top level methods.
   !
@@ -778,7 +780,6 @@ contains
 
   end subroutine amg_s_apply1_vect
 
-
   subroutine amg_s_apply2v(prec,x,y,desc_data,info,trans,work)
     implicit none
     type(psb_desc_type),intent(in)    :: desc_data
@@ -843,7 +844,6 @@ contains
   subroutine amg_s_dump(prec,info,istart,iend,iproc,prefix,head,&
        & ac,rp,smoother,solver,tprol,&
        & global_num)
-
     implicit none
     class(amg_sprec_type), intent(in)     :: prec
     integer(psb_ipk_), intent(out)          :: info
@@ -860,7 +860,7 @@ contains
     info = 0
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    iln = size(prec%precv)
+    iln = prec%get_nlevs()
     if (present(istart)) then
       il1 = max(1,istart)
     else
@@ -884,7 +884,6 @@ contains
   end subroutine amg_s_dump
 
   subroutine amg_s_cnv(prec,info,amold,vmold,imold)
-
     implicit none
     class(amg_sprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)       :: info
@@ -896,7 +895,7 @@ contains
 
     info = psb_success_
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         if (info == psb_success_ ) &
              & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
       end do
@@ -905,7 +904,6 @@ contains
   end subroutine amg_s_cnv
 
   subroutine amg_s_clone(prec,precout,info)
-
     implicit none
     class(amg_sprec_type), intent(inout) :: prec
     class(psb_sprec_type), intent(inout) :: precout
@@ -917,7 +915,6 @@ contains
   end subroutine amg_s_clone
 
   subroutine amg_s_inner_clone(prec,precout,info)
-
     implicit none
     class(amg_sprec_type), intent(inout)         :: prec
     class(psb_sprec_type), target, intent(inout) :: precout
@@ -933,8 +930,9 @@ contains
       pout%ctxt          = prec%ctxt
       pout%ag_data       = prec%ag_data
       pout%outer_sweeps  = prec%outer_sweeps
+      pout%nlevs         = prec%nlevs
       if (allocated(prec%precv)) then
-        ln = size(prec%precv)
+        ln = prec%get_nlevs()
         allocate(pout%precv(ln),stat=info)
         if (info /= psb_success_) goto 9999
         if (ln >= 1) then
@@ -942,6 +940,7 @@ contains
         end if
         do lev=2, ln
           if (info /= psb_success_) exit
+!!$          write(0,*) 'Inner_clone must be checked and reimplemented! '
           call prec%precv(lev)%clone(pout%precv(lev),info)
           if (info == psb_success_) then
             pout%precv(lev)%base_a       => pout%precv(lev)%ac
@@ -1006,7 +1005,7 @@ contains
     integer(psb_ipk_), intent(out)        :: info
     class(psb_s_base_vect_type), intent(in), optional  :: vmold
     !
-    ! In MLD the DESC optional argument is ignored, since
+    ! In AMG the DESC optional argument is ignored, since
     ! the necessary info is contained in the various entries of the
     ! PRECV component.
     type(psb_desc_type), intent(in), optional  :: desc
@@ -1021,7 +1020,7 @@ contains
     if (psb_errstatus_fatal()) then
       info = psb_err_internal_error_; goto 9999
     end if
-    nlev   = size(prec%precv)
+    nlev   = prec%get_nlevs()
     level = 1
     do level = 1, nlev
       call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@@ -1044,7 +1043,6 @@ contains
   subroutine amg_s_free_wrk(prec,info)
     use psb_base_mod
     implicit none
-
     ! Arguments
     class(amg_sprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -1061,7 +1059,7 @@ contains
     end if
 
     if (allocated(prec%precv)) then
-      nlev   = size(prec%precv)
+      nlev   = prec%get_nlevs()
       do level = 1, nlev
         call prec%precv(level)%free_wrk(info)
       end do

amgprec/amg_s_slu_solver.F90

@@ -124,7 +124,7 @@ module amg_s_slu_solver
       Implicit None
 
       ! Arguments
-      type(psb_sspmat_type), intent(in), target           :: a
+      type(psb_sspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_s_slu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_ainv_solver.F90

@@ -103,7 +103,7 @@ module amg_z_ainv_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_z_ainv_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -225,7 +225,7 @@ module amg_z_ainv_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
            & psb_d_vect_type, psb_z_base_vect_type, psb_dpk_, psb_ipk_
       implicit none
-      type(psb_zspmat_type), intent(in), target   :: a
+      type(psb_zspmat_type), intent(inout), target  :: a
       integer(psb_ipk_), intent(in)               :: fillin,alg
       real(psb_dpk_), intent(in)                  :: thresh
       type(psb_zspmat_type), intent(inout)        :: wmat, zmat

amgprec/amg_z_as_smoother.f90

@@ -230,7 +230,7 @@ module amg_z_as_smoother
            & psb_desc_type, psb_z_base_sparse_mat, psb_ipk_,&
            & psb_i_base_vect_type
       implicit none 
-      type(psb_zspmat_type), intent(in), target        :: a
+      type(psb_zspmat_type), intent(inout), target     :: a
       Type(psb_desc_type), Intent(inout)                 :: desc_a 
       class(amg_z_as_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                     :: info

amgprec/amg_z_base_aggregator_mod.f90

@@ -89,7 +89,7 @@ module amg_z_base_aggregator_mod
     procedure, pass(ag) :: bld_tprol   => amg_z_base_aggregator_build_tprol
     procedure, pass(ag) :: mat_bld     => amg_z_base_aggregator_mat_bld
     procedure, pass(ag) :: mat_asb     => amg_z_base_aggregator_mat_asb
-    procedure, pass(ag) :: bld_map     => amg_z_base_aggregator_bld_map
+    procedure, pass(ag) :: bld_linmap  => amg_z_base_aggregator_bld_linmap
     procedure, pass(ag) :: update_next => amg_z_base_aggregator_update_next
     procedure, pass(ag) :: clone       => amg_z_base_aggregator_clone
     procedure, pass(ag) :: free        => amg_z_base_aggregator_free
@@ -458,7 +458,7 @@ contains
   end subroutine amg_z_base_aggregator_mat_asb
 
   !
-  !> Function   bld_map
+  !> Function   bld_linmap
   !! \memberof  amg_z_base_aggregator_type
   !! \brief     Build linear map between hierarchy levels 
   !!
@@ -473,7 +473,7 @@ contains
   !!  \param map       The output map
   !!  \param info    Return code
   !!  
-  subroutine  amg_z_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
+  subroutine  amg_z_base_aggregator_bld_linmap(ag,desc_a,desc_ac,ilaggr,nlaggr,&
        & op_restr,op_prol,map,info)
     use psb_base_mod
     implicit none
@@ -484,7 +484,7 @@ contains
     type(psb_zlinmap_type), intent(out)    :: map
     integer(psb_ipk_), intent(out)       :: info
     integer(psb_ipk_) :: err_act
-    character(len=20) :: name='z_base_aggregator_bld_map'
+    character(len=20) :: name='z_base_aggregator_bld_linmap'
 
     info = psb_success_
     call psb_erractionsave(err_act)
@@ -508,7 +508,6 @@ contains
 9999 call psb_error_handler(err_act)
 
     return
-  end subroutine amg_z_base_aggregator_bld_map
-
+  end subroutine amg_z_base_aggregator_bld_linmap
 
 end module amg_z_base_aggregator_mod

amgprec/amg_z_base_smoother_mod.f90

@@ -237,7 +237,7 @@ module amg_z_base_smoother_mod
            & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
            & amg_z_base_smoother_type, psb_ipk_, psb_i_base_vect_type
       ! Arguments
-      type(psb_zspmat_type), intent(in), target     :: a
+      type(psb_zspmat_type), intent(inout), target  :: a
       Type(psb_desc_type), Intent(inout)              :: desc_a 
       class(amg_z_base_smoother_type), intent(inout) :: sm 
       integer(psb_ipk_), intent(out)                  :: info

amgprec/amg_z_base_solver_mod.f90

@@ -170,7 +170,7 @@ module amg_z_base_solver_mod
       Implicit None
       
       ! Arguments
-      type(psb_zspmat_type), intent(in), target             :: a
+      type(psb_zspmat_type), intent(inout), target          :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_z_base_solver_type), intent(inout)          :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_z_diag_solver.f90

@@ -119,7 +119,7 @@ module amg_z_diag_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
            & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
            & amg_z_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_z_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info
@@ -331,7 +331,7 @@ module amg_z_l1_diag_solver
       import :: psb_desc_type, psb_zspmat_type,  psb_z_base_sparse_mat, &
            & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
            & amg_z_l1_diag_solver_type, psb_ipk_, psb_i_base_vect_type      
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                    :: desc_a 
       class(amg_z_l1_diag_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/amg_z_gs_solver.f90

@@ -181,7 +181,7 @@ module amg_z_gs_solver
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_gs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -195,7 +195,7 @@ module amg_z_gs_solver
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_bwgs_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_id_solver.f90

@@ -123,7 +123,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_zspmat_type), intent(in), target           :: a
+    type(psb_zspmat_type), intent(inout), target        :: a
     Type(psb_desc_type), Intent(inout)                  :: desc_a 
     class(amg_z_id_solver_type), intent(inout)          :: sv
     integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_ilu_solver.f90

@@ -144,7 +144,7 @@ module amg_z_ilu_solver
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_ilu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_inner_mod.f90

@@ -56,7 +56,7 @@ module amg_z_inner_mod
            & psb_dpk_, psb_z_base_sparse_mat, psb_z_base_vect_type, psb_ipk_
       import :: amg_zprec_type
       implicit none
-      type(psb_zspmat_type), intent(in), target          :: a
+      type(psb_zspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       type(amg_zprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -67,7 +67,7 @@ module amg_z_inner_mod
   end interface amg_mlprec_bld
 
   interface amg_mlprec_aply
-    subroutine amg_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
+    subroutine amg_zmlprec_aply_a(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_
       import :: amg_zprec_type
       implicit none 
@@ -79,7 +79,7 @@ module amg_z_inner_mod
       character,intent(in)               :: trans
       complex(psb_dpk_),target             :: work(:)
       integer(psb_ipk_), intent(out)     :: info
-    end subroutine amg_zmlprec_aply
+    end subroutine amg_zmlprec_aply_a
     subroutine amg_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
       import :: psb_zspmat_type, psb_desc_type, &
            & psb_dpk_, psb_z_vect_type, psb_ipk_

amgprec/amg_z_invk_solver.f90

@@ -94,7 +94,7 @@ module amg_z_invk_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_z_invk_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_invt_solver.f90

@@ -94,7 +94,7 @@ module amg_z_invt_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_z_invt_solver_type), intent(inout)        :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_jac_smoother.f90

@@ -151,7 +151,7 @@ module amg_z_jac_smoother
       import :: psb_desc_type, amg_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, &
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_z_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -274,7 +274,7 @@ module amg_z_jac_smoother
       import :: psb_desc_type, amg_z_l1_jac_smoother_type, psb_z_vect_type, &
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_z_l1_jac_smoother_type), intent(inout)       :: sm
       integer(psb_ipk_), intent(out)                      :: info
@@ -345,6 +345,7 @@ contains
 
     if (allocated(sm%sv)) then
       call sm%sv%free(info)
+      if (info == psb_success_) deallocate(sm%sv,stat=info)
       if (info /= psb_success_) then
         info = psb_err_alloc_dealloc_
         call psb_errpush(info,name)

amgprec/amg_z_jac_solver.f90

@@ -143,7 +143,7 @@ module amg_z_jac_solver
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info
@@ -160,7 +160,7 @@ module amg_z_jac_solver
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none 
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_l1_jac_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_krm_solver.f90

@@ -174,7 +174,7 @@ module amg_z_krm_solver
            & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
            & psb_ipk_, psb_i_base_vect_type
       implicit none
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a
       class(amg_z_krm_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_mumps_solver.F90

@@ -163,7 +163,7 @@ module amg_z_mumps_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_mumps_solver_type), intent(inout)       :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_onelev_mod.f90

@@ -154,13 +154,35 @@ module amg_z_onelev_mod
   private :: z_wrk_alloc, z_wrk_free, &
        & z_wrk_clone, z_wrk_move_alloc, z_wrk_cnv, z_wrk_sizeof
 
+  !
+  ! Remap.
+  !  This keeps track of remapping.
+  !  Logic here is as follows:
+  !  1. AC_PRE_REMAP  need to figure out if we
+  !                really need it.
+  !  2. DESC_AC_PRE_REMAP  contains the descriptor before
+  !                remapping. Meaning that it is possible
+  !                to implement the RESTRICTOR operator by
+  !                 a. Doing LINMAP_U2V onto this one
+  !                 b. For each process, send the data to
+  !                    IDEST.
+  !                This assumes that remapping goes by
+  !                a factor of 2.
+  !                For the PROLONGATOR operators, we first
+  !                use DESC_AC, then split and send onto
+  !                the processes in DESC_AC_PRE_REMAP.
+  !
+  ! To be fixed: what happens if NP the starting processes
+  ! is not an even number? Coordinate with _X_remap in PSBLAS
+  ! 
   type amg_z_remap_data_type
     type(psb_zspmat_type)          :: ac_pre_remap
     type(psb_desc_type)            :: desc_ac_pre_remap
     integer(psb_ipk_)              :: idest
     integer(psb_ipk_), allocatable :: isrc(:), nrsrc(:), naggr(:)
   contains
-    procedure, pass(rmp) :: clone   => z_remap_data_clone
+    procedure, pass(rmp) :: clone      => z_remap_data_clone
+    procedure, pass(rmp) :: move_alloc => z_remap_move_alloc
   end type amg_z_remap_data_type
 
   type amg_z_onelev_type
@@ -206,7 +228,7 @@ module amg_z_onelev_mod
     procedure, pass(lv) :: get_wrksz => z_base_onelev_get_wrksize
     procedure, pass(lv) :: allocate_wrk   => z_base_onelev_allocate_wrk
     procedure, pass(lv) :: free_wrk       => z_base_onelev_free_wrk
-    procedure, nopass   :: stringval => amg_stringval
+    procedure, nopass   :: stringval  => amg_stringval
     procedure, pass(lv) :: move_alloc => z_base_onelev_move_alloc
 
 
@@ -618,7 +640,7 @@ contains
     ! Arguments
     class(amg_z_onelev_type), target, intent(inout) :: lv
     class(amg_z_onelev_type), target, intent(inout) :: lvout
-    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_ipk_), intent(out)                  :: info
 
     info = psb_success_
     if (allocated(lv%sm)) then
@@ -684,6 +706,7 @@ contains
     if (info == psb_success_) call psb_move_alloc(lv%tprol,b%tprol,info)
     if (info == psb_success_) call psb_move_alloc(lv%desc_ac,b%desc_ac,info)
     if (info == psb_success_) call psb_move_alloc(lv%linmap,b%linmap,info)
+    if (info == psb_success_) call lv%remap_data%move_alloc(b%remap_data,info)
     b%base_a    => lv%base_a
     b%base_desc => lv%base_desc
 
@@ -738,6 +761,7 @@ contains
     info = psb_success_
     nwv = lv%get_wrksz()
     if (.not.allocated(lv%wrk)) allocate(lv%wrk,stat=info)
+!!$    write(0,*) 'From allocate_wrk :',lv%remap_data%desc_ac_pre_remap%is_asb()
     if (info == 0) then
       if (lv%remap_data%desc_ac_pre_remap%is_asb()) then
         !
@@ -785,47 +809,32 @@ contains
 
     info = psb_success_
     call wk%free(info)
-    if  (present(desc2)) then
-!!$      write(0,*) 'Check on wrk_alloc 2',&
-!!$           & desc2%get_local_rows(), desc%get_local_rows(),&
-!!$           & desc2%get_local_cols(),desc%get_local_cols()
-!!$      flush(0)
+!!$    write(0,*) 'wrk_alloc D: "',trim(desc%get_fmt()),'"',&
+!!$         & present(desc2),desc%is_valid()
+
+    allocate(wk%wv(nwv),stat=info)    
+    if  (present(desc2).and.(desc%is_valid())) then
+!!$      write(0,*) 'wrk_alloc D2:',desc2%get_fmt(),desc2%is_asb()
       if (desc2%get_local_cols()>desc%get_local_cols()) then
-        call psb_geasb(wk%vx2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc2,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc2,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
       else
-!!$        write(0,*) 'Check on wrk_alloc 1.5 ',&
-!!$             & desc%get_local_rows(),&
-!!$             & desc%get_local_cols()
-        call psb_geasb(wk%vx2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vy2l,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vtx,desc,info,&
-             & scratch=.true.,mold=vmold)
-        call psb_geasb(wk%vty,desc,info,&
-             & scratch=.true.,mold=vmold)
-        allocate(wk%wv(nwv),stat=info)
-        do i=1,nwv
-          call psb_geasb(wk%wv(i),desc,info,&
-               & scratch=.true.,mold=vmold)
-        end do
+        call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)
       end if
-    else
-!!$      write(0,*) 'Check on wrk_alloc 1 ',&
-!!$           & desc%get_local_rows(),&
-!!$           & desc%get_local_cols()
+    else if (present(desc2)) then
+      call inner_do_wrk_alloc(wk,nwv,desc2,vmold=vmold)
+    else if (desc%is_valid()) then
+      call inner_do_wrk_alloc(wk,nwv,desc,vmold=vmold)      
+    end if
+
+  contains
+    subroutine inner_do_wrk_alloc(wk,nwv,desc,vmold)
+      class(amg_zmlprec_wrk_type), target, intent(inout) :: wk
+      integer(psb_ipk_), intent(in)                     :: nwv
+      type(psb_desc_type), intent(in)                   :: desc
+      class(psb_z_base_vect_type), intent(in), optional  :: vmold
+      
+      integer(psb_ipk_) :: i
+      
       call psb_geasb(wk%vx2l,desc,info,&
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vy2l,desc,info,&
@@ -834,12 +843,11 @@ contains
            & scratch=.true.,mold=vmold)
       call psb_geasb(wk%vty,desc,info,&
            & scratch=.true.,mold=vmold)
-      allocate(wk%wv(nwv),stat=info)
       do i=1,nwv
         call psb_geasb(wk%wv(i),desc,info,&
              & scratch=.true.,mold=vmold)
       end do
-    end if
+    end subroutine inner_do_wrk_alloc
   end subroutine z_wrk_alloc
 
   subroutine z_wrk_free(wk,info)
@@ -992,4 +1000,25 @@ contains
     call psb_safe_ab_cpy(rmp%nrsrc,remap_out%nrsrc,info)
   end subroutine z_remap_data_clone
 
+  subroutine z_remap_move_alloc(rmp, remap_out, info)
+    use psb_base_mod
+    implicit none
+    ! Arguments
+    class(amg_z_remap_data_type), target, intent(inout) :: rmp
+    class(amg_z_remap_data_type), target, intent(inout) :: remap_out
+    integer(psb_ipk_), intent(out)                    :: info
+    !
+    integer(psb_ipk_) :: i
+
+    info = psb_success_
+
+    call psb_move_alloc(rmp%ac_pre_remap,remap_out%ac_pre_remap,info)
+    if (info == psb_success_) &
+         & call psb_move_alloc(rmp%desc_ac_pre_remap,remap_out%desc_ac_pre_remap,info)
+    remap_out%idest = rmp%idest
+    call move_alloc(rmp%isrc,remap_out%isrc)
+    call move_alloc(rmp%nrsrc,remap_out%nrsrc)
+    call move_alloc(rmp%naggr,remap_out%naggr)
+  end subroutine z_remap_move_alloc
+  
 end module amg_z_onelev_mod

amgprec/amg_z_prec_type.f90

@@ -102,6 +102,7 @@ module amg_z_prec_type
     ! The multilevel hierarchy
     !
     type(amg_z_onelev_type), allocatable :: precv(:)
+    integer(psb_ipk_)                    :: nlevs
   contains
     procedure, pass(prec)               :: psb_z_apply2_vect => amg_z_apply2_vect
     procedure, pass(prec)               :: psb_z_apply1_vect => amg_z_apply1_vect
@@ -119,6 +120,7 @@ module amg_z_prec_type
     procedure, pass(prec)               :: cmp_complexity => amg_z_cmp_compl
     procedure, pass(prec)               :: get_avg_cr => amg_z_get_avg_cr
     procedure, pass(prec)               :: cmp_avg_cr => amg_z_cmp_avg_cr
+    procedure, pass(prec)               :: set_nlevs  => amg_z_set_nlevs
     procedure, pass(prec)               :: get_nlevs  => amg_z_get_nlevs
     procedure, pass(prec)               :: get_nzeros => amg_z_get_nzeros
     procedure, pass(prec)               :: sizeof => amg_zprec_sizeof
@@ -311,7 +313,7 @@ module amg_z_prec_type
            & psb_z_base_sparse_mat, psb_z_base_vect_type, &
            & psb_i_base_vect_type, amg_zprec_type, psb_ipk_
       implicit none
-      type(psb_zspmat_type), intent(in), target          :: a
+      type(psb_zspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_zprec_type), intent(inout), target       :: prec
       integer(psb_ipk_), intent(out)                       :: info
@@ -323,14 +325,15 @@ module amg_z_prec_type
   end interface amg_precbld
 
   interface amg_hierarchy_bld
-    subroutine amg_z_hierarchy_bld(a,desc_a,prec,info)
+    subroutine amg_z_hierarchy_bld(a,desc_a,prec,info,cpymat)
       import :: psb_zspmat_type, psb_desc_type, psb_dpk_, &
            & amg_zprec_type, psb_ipk_
       implicit none
-      type(psb_zspmat_type), intent(in), target          :: a
+      type(psb_zspmat_type), intent(inout), target       :: a
       type(psb_desc_type), intent(inout), target           :: desc_a
       class(amg_zprec_type), intent(inout), target        :: prec
       integer(psb_ipk_), intent(out)                       :: info
+      logical, intent(in), optional                        :: cpymat
       !      character, intent(in),optional             :: upd
     end subroutine amg_z_hierarchy_bld
   end interface amg_hierarchy_bld
@@ -434,10 +437,19 @@ contains
     class(amg_zprec_type), intent(in) :: prec
     integer(psb_ipk_) :: val
     val = 0
-    if (allocated(prec%precv)) then
-      val = size(prec%precv)
-    end if
+!!$    if (allocated(prec%precv)) then
+!!$      val = size(prec%precv)
+!!$    end if
+    val = prec%nlevs
+!!$    write(0,*) ' NLEVS: ',prec%nlevs, val,size(prec%precv)
   end function amg_z_get_nlevs
+
+  subroutine amg_z_set_nlevs(prec,nl) 
+    implicit none
+    class(amg_zprec_type), intent(inout) :: prec
+    integer(psb_ipk_) :: nl
+    prec%nlevs = nl
+  end subroutine amg_z_set_nlevs
   !
   ! Function returning the size of the amg_prec_type data structure
   ! in bytes or in number of nonzeros of the operator(s) involved.
@@ -508,7 +520,7 @@ contains
 
     real(psb_dpk_) :: num, den, nmin
     type(psb_ctxt_type) :: ctxt
-    integer(psb_ipk_)   :: il
+    integer(psb_ipk_)   :: il,nl
 
     num = -done
     den = done
@@ -518,7 +530,10 @@ contains
       num = prec%precv(il)%base_a%get_nzeros()
       if (num >= dzero) then
         den = num
-        do il=2,size(prec%precv)
+        nl = prec%get_nlevs()
+!!$        write(0,*) 'Inside cmp_compl ',nl,size(prec%precv)
+        do il=2, nl
+!!$          write(0,*) '                 ',il,associated(prec%precv(il)%base_a)
           num = num + max(0,prec%precv(il)%base_a%get_nzeros())
         end do
       end if
@@ -549,7 +564,6 @@ contains
   end function amg_z_get_avg_cr
 
   subroutine amg_z_cmp_avg_cr(prec)
-
     implicit none
     class(amg_zprec_type), intent(inout) :: prec
 
@@ -557,17 +571,18 @@ contains
     type(psb_ctxt_type) :: ctxt
     integer(psb_ipk_)   :: il, nl, iam, np
 
-
     avgcr = dzero
+    nl = prec%get_nlevs()
+    do il=2,nl
+      if (prec%precv(il)%base_desc%is_ok()) then 
+        ctxt = prec%precv(il)%base_desc%get_ctxt()
+        call psb_info(ctxt,iam,np)
+        if (iam >=0)  avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
+      end if
+    end do
+    avgcr = avgcr / (nl-1)
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    if (allocated(prec%precv)) then
-      nl = size(prec%precv)
-      do il=2,nl
-        avgcr = avgcr + max(dzero,prec%precv(il)%szratio)
-      end do
-      avgcr = avgcr / (nl-1)
-    end if
     call psb_sum(ctxt,avgcr)
     prec%ag_data%avg_cr = avgcr/np
   end subroutine amg_z_cmp_avg_cr
@@ -585,9 +600,7 @@ contains
   !             error code.
   !
   subroutine amg_zprecfree(p,info)
-
     implicit none
-
     ! Arguments
     type(amg_zprec_type), intent(inout) :: p
     integer(psb_ipk_), intent(out)        :: info
@@ -613,9 +626,7 @@ contains
   end subroutine amg_zprecfree
 
   subroutine amg_z_prec_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_zprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -647,9 +658,7 @@ contains
   end subroutine amg_z_prec_free
 
   subroutine amg_z_smoothers_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_zprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -665,12 +674,8 @@ contains
       info = psb_err_internal_error_; goto 9999
     end if
 
-    !
-    ! In the internals, do FREE on components,
-    ! but do not deallocate them
-    !
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         call prec%precv(i)%free_smoothers(info)
       end do
     end if
@@ -684,9 +689,7 @@ contains
   end subroutine amg_z_smoothers_free
 
   subroutine amg_z_hierarchy_free(prec,info)
-
     implicit none
-
     ! Arguments
     class(amg_zprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -713,7 +716,6 @@ contains
 
   end subroutine amg_z_hierarchy_free
 
-
   !
   ! Top level methods.
   !
@@ -778,7 +780,6 @@ contains
 
   end subroutine amg_z_apply1_vect
 
-
   subroutine amg_z_apply2v(prec,x,y,desc_data,info,trans,work)
     implicit none
     type(psb_desc_type),intent(in)    :: desc_data
@@ -843,7 +844,6 @@ contains
   subroutine amg_z_dump(prec,info,istart,iend,iproc,prefix,head,&
        & ac,rp,smoother,solver,tprol,&
        & global_num)
-
     implicit none
     class(amg_zprec_type), intent(in)     :: prec
     integer(psb_ipk_), intent(out)          :: info
@@ -860,7 +860,7 @@ contains
     info = 0
     ctxt = prec%ctxt
     call psb_info(ctxt,iam,np)
-    iln = size(prec%precv)
+    iln = prec%get_nlevs()
     if (present(istart)) then
       il1 = max(1,istart)
     else
@@ -884,7 +884,6 @@ contains
   end subroutine amg_z_dump
 
   subroutine amg_z_cnv(prec,info,amold,vmold,imold)
-
     implicit none
     class(amg_zprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)       :: info
@@ -896,7 +895,7 @@ contains
 
     info = psb_success_
     if (allocated(prec%precv)) then
-      do i=1,size(prec%precv)
+      do i=1,prec%get_nlevs()
         if (info == psb_success_ ) &
              & call prec%precv(i)%cnv(info,amold=amold,vmold=vmold,imold=imold)
       end do
@@ -905,7 +904,6 @@ contains
   end subroutine amg_z_cnv
 
   subroutine amg_z_clone(prec,precout,info)
-
     implicit none
     class(amg_zprec_type), intent(inout) :: prec
     class(psb_zprec_type), intent(inout) :: precout
@@ -917,7 +915,6 @@ contains
   end subroutine amg_z_clone
 
   subroutine amg_z_inner_clone(prec,precout,info)
-
     implicit none
     class(amg_zprec_type), intent(inout)         :: prec
     class(psb_zprec_type), target, intent(inout) :: precout
@@ -933,8 +930,9 @@ contains
       pout%ctxt          = prec%ctxt
       pout%ag_data       = prec%ag_data
       pout%outer_sweeps  = prec%outer_sweeps
+      pout%nlevs         = prec%nlevs
       if (allocated(prec%precv)) then
-        ln = size(prec%precv)
+        ln = prec%get_nlevs()
         allocate(pout%precv(ln),stat=info)
         if (info /= psb_success_) goto 9999
         if (ln >= 1) then
@@ -942,6 +940,7 @@ contains
         end if
         do lev=2, ln
           if (info /= psb_success_) exit
+!!$          write(0,*) 'Inner_clone must be checked and reimplemented! '
           call prec%precv(lev)%clone(pout%precv(lev),info)
           if (info == psb_success_) then
             pout%precv(lev)%base_a       => pout%precv(lev)%ac
@@ -1006,7 +1005,7 @@ contains
     integer(psb_ipk_), intent(out)        :: info
     class(psb_z_base_vect_type), intent(in), optional  :: vmold
     !
-    ! In MLD the DESC optional argument is ignored, since
+    ! In AMG the DESC optional argument is ignored, since
     ! the necessary info is contained in the various entries of the
     ! PRECV component.
     type(psb_desc_type), intent(in), optional  :: desc
@@ -1021,7 +1020,7 @@ contains
     if (psb_errstatus_fatal()) then
       info = psb_err_internal_error_; goto 9999
     end if
-    nlev   = size(prec%precv)
+    nlev   = prec%get_nlevs()
     level = 1
     do level = 1, nlev
       call prec%precv(level)%allocate_wrk(info,vmold=vmold)
@@ -1044,7 +1043,6 @@ contains
   subroutine amg_z_free_wrk(prec,info)
     use psb_base_mod
     implicit none
-
     ! Arguments
     class(amg_zprec_type), intent(inout) :: prec
     integer(psb_ipk_), intent(out)        :: info
@@ -1061,7 +1059,7 @@ contains
     end if
 
     if (allocated(prec%precv)) then
-      nlev   = size(prec%precv)
+      nlev   = prec%get_nlevs()
       do level = 1, nlev
         call prec%precv(level)%free_wrk(info)
       end do

amgprec/amg_z_slu_solver.F90

@@ -124,7 +124,7 @@ module amg_z_slu_solver
       Implicit None
 
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_slu_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                      :: info

amgprec/amg_z_sludist_solver.F90

@@ -259,7 +259,7 @@ contains
     Implicit None
 
     ! Arguments
-    type(psb_zspmat_type), intent(in), target           :: a
+    type(psb_zspmat_type), intent(inout), target        :: a
     Type(psb_desc_type), Intent(inout)                  :: desc_a 
     class(amg_z_sludist_solver_type), intent(inout)     :: sv
     integer, intent(out)                                :: info

amgprec/amg_z_umf_solver.F90

@@ -163,7 +163,7 @@ module amg_z_umf_solver
       Implicit None
       
       ! Arguments
-      type(psb_zspmat_type), intent(in), target           :: a
+      type(psb_zspmat_type), intent(inout), target        :: a
       Type(psb_desc_type), Intent(inout)                  :: desc_a 
       class(amg_z_umf_solver_type), intent(inout)         :: sv
       integer(psb_ipk_), intent(out)                        :: info

amgprec/impl/Makefile

@@ -25,22 +25,22 @@ MPCOBJS=amg_dslud_interface.o amg_zslud_interface.o
 
 DINNEROBJS= amg_dmlprec_bld.o amg_dfile_prec_descr.o amg_dfile_prec_memory_use.o \
 	amg_d_smoothers_bld.o amg_d_hierarchy_bld.o amg_d_hierarchy_rebld.o \
-	amg_dmlprec_aply.o   \
+	amg_dmlprec_aply.o amg_dmlprec_aply_a.o \
 	$(DMPFOBJS) amg_d_extprol_bld.o
 
 SINNEROBJS=   amg_smlprec_bld.o amg_sfile_prec_descr.o amg_sfile_prec_memory_use.o \
 	amg_s_smoothers_bld.o amg_s_hierarchy_bld.o amg_s_hierarchy_rebld.o \
-	amg_smlprec_aply.o   \
+	amg_smlprec_aply.o amg_smlprec_aply_a.o \
 	$(SMPFOBJS) amg_s_extprol_bld.o 
 
 ZINNEROBJS=   amg_zmlprec_bld.o amg_zfile_prec_descr.o amg_zfile_prec_memory_use.o \
 	amg_z_smoothers_bld.o amg_z_hierarchy_bld.o amg_z_hierarchy_rebld.o \
-	amg_zmlprec_aply.o   \
+	amg_zmlprec_aply.o amg_zmlprec_aply_a.o \
 	$(ZMPFOBJS)  amg_z_extprol_bld.o
 
 CINNEROBJS=   amg_cmlprec_bld.o amg_cfile_prec_descr.o amg_cfile_prec_memory_use.o \
 	amg_c_smoothers_bld.o amg_c_hierarchy_bld.o amg_c_hierarchy_rebld.o \
-	amg_cmlprec_aply.o   \
+	amg_cmlprec_aply.o amg_cmlprec_aply_a.o \
 	$(CMPFOBJS)  amg_c_extprol_bld.o
 
 INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS)

amgprec/impl/aggregator/amg_c_dec_aggregator_mat_asb.f90

@@ -114,63 +114,67 @@ subroutine  amg_c_dec_aggregator_mat_asb(ag,parms,a,desc_a,&
   info  = psb_success_
   ctxt = desc_a%get_context()
   call psb_info(ctxt,me,np)
+  
+  if (me >=0) then 
+    select case(parms%coarse_mat)
 
-  select case(parms%coarse_mat)
+    case(amg_distr_mat_) 
 
-  case(amg_distr_mat_) 
+      call ac%cscnv(info,type='csr')
+      call op_prol%cscnv(info,type='csr')
+      call op_restr%cscnv(info,type='csr')
 
-    call ac%cscnv(info,type='csr')
-    call op_prol%cscnv(info,type='csr')
-    call op_restr%cscnv(info,type='csr')
+      if (debug_level >= psb_debug_outer_) &
+           & write(debug_unit,*) me,' ',trim(name),&
+           & 'Done ac '
 
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done ac '
+    case(amg_repl_mat_) 
+      !
+      ! We are assuming here that an c matrix
+      ! can hold all entries
+      !
+      if (desc_ac%get_global_rows() < huge(1_psb_ipk_) ) then 
+        ntaggr = desc_ac%get_global_rows()
+        i_nr   = ntaggr
+      else
+        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='invalid amg_coarse_mat_')
+        goto 9999
+      end if
 
-  case(amg_repl_mat_) 
-    !
-    ! We are assuming here that an c matrix
-    ! can hold all entries
-    !
-    if (desc_ac%get_global_rows() < huge(1_psb_ipk_) ) then 
-      ntaggr = desc_ac%get_global_rows()
-      i_nr   = ntaggr
-    else
+      call op_prol%mv_to(tmpcoo)
+      nzl = tmpcoo%get_nzeros()
+      call psb_loc_to_glob(tmpcoo%ja(1:nzl),desc_ac,info,'I')
+      call tmpcoo%set_ncols(i_nr)
+      call op_prol%mv_from(tmpcoo)
+
+      call op_restr%mv_to(tmpcoo)
+      nzl = tmpcoo%get_nzeros()
+      call psb_loc_to_glob(tmpcoo%ia(1:nzl),desc_ac,info,'I')
+      call tmpcoo%set_nrows(i_nr)
+      call op_restr%mv_from(tmpcoo)
+
+
+      call psb_gather(tmp_ac,ac,desc_ac,info,root=-ione,&
+           & dupl=psb_dupl_add_,keeploc=.false.)
+      call tmp_ac%mv_to(tmpcoo)
+      call ac%mv_from(tmpcoo)
+
+      call psb_cdall(ctxt,desc_ac,info,mg=ntaggr,repl=.true.)
+      if (info == psb_success_) call psb_cdasb(desc_ac,info)
+
+      if (info /= psb_success_) goto 9999
+
+    case default 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_coarse_mat_')
       goto 9999
-    end if
-    
-    call op_prol%mv_to(tmpcoo)
-    nzl = tmpcoo%get_nzeros()
-    call psb_loc_to_glob(tmpcoo%ja(1:nzl),desc_ac,info,'I')
-    call tmpcoo%set_ncols(i_nr)
-    call op_prol%mv_from(tmpcoo)
-
-    call op_restr%mv_to(tmpcoo)
-    nzl = tmpcoo%get_nzeros()
-    call psb_loc_to_glob(tmpcoo%ia(1:nzl),desc_ac,info,'I')
-    call tmpcoo%set_nrows(i_nr)
-    call op_restr%mv_from(tmpcoo)
-
-
-    call psb_gather(tmp_ac,ac,desc_ac,info,root=-ione,&
-         & dupl=psb_dupl_add_,keeploc=.false.)
-    call tmp_ac%mv_to(tmpcoo)
-    call ac%mv_from(tmpcoo)
-    
-    call psb_cdall(ctxt,desc_ac,info,mg=ntaggr,repl=.true.)
-    if (info == psb_success_) call psb_cdasb(desc_ac,info)
-
-    if (info /= psb_success_) goto 9999
-
-  case default 
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='invalid amg_coarse_mat_')
-    goto 9999
-  end select
-
-
+    end select
+  else
+    call op_prol%allocate(izero,izero,info)
+    call op_restr%allocate(izero,izero,info)
+    call ac%allocate(izero,izero,info)
+  end if
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/amg_c_dec_aggregator_mat_bld.f90

@@ -139,7 +139,7 @@ subroutine  amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   use amg_c_prec_type, amg_protect_name => amg_c_dec_aggregator_mat_bld
   use amg_c_inner_mod 
   implicit none
-  
+
   class(amg_c_dec_aggregator_type), target, intent(inout) :: ag
   type(amg_sml_parms), intent(inout)      :: parms 
   type(psb_cspmat_type), intent(in)    :: a
@@ -169,39 +169,46 @@ subroutine  amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   ctxt = desc_a%get_context()
   call psb_info(ctxt,me,np)
 
-  !
-  ! Build the coarse-level matrix from the fine-level one, starting from 
-  ! the mapping defined by amg_aggrmap_bld and applying the aggregation
-  ! algorithm specified by 
-  !
-  select case (parms%aggr_prol)
-  case (amg_no_smooth_) 
+  if (me >=0) then 
+    !
+    ! Build the coarse-level matrix from the fine-level one, starting from 
+    ! the mapping defined by amg_aggrmap_bld and applying the aggregation
+    ! algorithm specified by 
+    !
+    select case (parms%aggr_prol)
+    case (amg_no_smooth_) 
 
-    call amg_caggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+      call amg_caggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+           nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
+    case(amg_smooth_prol_,amg_l1_smooth_prol_) 
 
-    call amg_caggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
-            op_restr,t_prol,info)
+      call amg_caggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+           ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+           op_restr,t_prol,info)
 
 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_caggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-  
-  case(amg_min_energy_) 
 
-    call amg_caggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    case(amg_min_energy_) 
 
-  case default
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='Invalid aggr kind')
-    goto 9999
+      call amg_caggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+           nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='Invalid aggr kind')
+      goto 9999
+
+    end select
+  else
+    call op_prol%allocate(izero,izero,info)
+    call op_restr%allocate(izero,izero,info)
+    call ac%allocate(izero,izero,info)
+  end if
 
-  end select
   if (info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat bld')
     goto 9999
@@ -214,5 +221,5 @@ subroutine  amg_c_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
 9999 call psb_error_handler(err_act)
   return
 
-  
+
 end subroutine amg_c_dec_aggregator_mat_bld

amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90

@@ -122,19 +122,23 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
   call amg_check_def(parms%aggr_ord,'Ordering',&
        &   amg_aggr_ord_nat_,is_legal_ml_aggr_ord)
   call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',szero,is_legal_s_aggr_thrs)
+  if (me >=0) then 
+    !
+    ! The decoupled aggregator based on SOC measures ignores
+    ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
+    !
+    if (do_timings) call psb_tic(idx_map_bld)
+    clean_zeros = ag%do_clean_zeros
+    call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+    if (do_timings) call psb_toc(idx_map_bld)
+    if (do_timings) call psb_tic(idx_map_tprol)
 
-  !
-  ! The decoupled aggregator based on SOC measures ignores
-  ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
-  !
-  if (do_timings) call psb_tic(idx_map_bld)
-  clean_zeros = ag%do_clean_zeros
-  call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
-  if (do_timings) call psb_toc(idx_map_bld)
-  if (do_timings) call psb_tic(idx_map_tprol)
-
-  if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
-  if (do_timings) call psb_toc(idx_map_tprol)
+    if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+    if (do_timings) call psb_toc(idx_map_tprol)
+  else
+    allocate(nlaggr(0),ilaggr(0))
+    call t_prol%allocate(lzero,lzero,info)
+  end if
   if (info /= psb_success_) then
     info=psb_err_from_subroutine_
     call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')

amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90

@@ -155,7 +155,7 @@ subroutine amg_caggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   debug_level = psb_get_debug_level()
 
   ctxt = desc_a%get_context()
-  icomm = desc_a%get_mpic()
+  icomm = ctxt%get_mpic()
 
   call psb_info(ctxt, me, np)
 

amgprec/impl/aggregator/amg_d_dec_aggregator_mat_asb.f90

@@ -114,63 +114,67 @@ subroutine  amg_d_dec_aggregator_mat_asb(ag,parms,a,desc_a,&
   info  = psb_success_
   ctxt = desc_a%get_context()
   call psb_info(ctxt,me,np)
+  
+  if (me >=0) then 
+    select case(parms%coarse_mat)
 
-  select case(parms%coarse_mat)
+    case(amg_distr_mat_) 
 
-  case(amg_distr_mat_) 
+      call ac%cscnv(info,type='csr')
+      call op_prol%cscnv(info,type='csr')
+      call op_restr%cscnv(info,type='csr')
 
-    call ac%cscnv(info,type='csr')
-    call op_prol%cscnv(info,type='csr')
-    call op_restr%cscnv(info,type='csr')
+      if (debug_level >= psb_debug_outer_) &
+           & write(debug_unit,*) me,' ',trim(name),&
+           & 'Done ac '
 
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done ac '
+    case(amg_repl_mat_) 
+      !
+      ! We are assuming here that an d matrix
+      ! can hold all entries
+      !
+      if (desc_ac%get_global_rows() < huge(1_psb_ipk_) ) then 
+        ntaggr = desc_ac%get_global_rows()
+        i_nr   = ntaggr
+      else
+        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='invalid amg_coarse_mat_')
+        goto 9999
+      end if
 
-  case(amg_repl_mat_) 
-    !
-    ! We are assuming here that an d matrix
-    ! can hold all entries
-    !
-    if (desc_ac%get_global_rows() < huge(1_psb_ipk_) ) then 
-      ntaggr = desc_ac%get_global_rows()
-      i_nr   = ntaggr
-    else
+      call op_prol%mv_to(tmpcoo)
+      nzl = tmpcoo%get_nzeros()
+      call psb_loc_to_glob(tmpcoo%ja(1:nzl),desc_ac,info,'I')
+      call tmpcoo%set_ncols(i_nr)
+      call op_prol%mv_from(tmpcoo)
+
+      call op_restr%mv_to(tmpcoo)
+      nzl = tmpcoo%get_nzeros()
+      call psb_loc_to_glob(tmpcoo%ia(1:nzl),desc_ac,info,'I')
+      call tmpcoo%set_nrows(i_nr)
+      call op_restr%mv_from(tmpcoo)
+
+
+      call psb_gather(tmp_ac,ac,desc_ac,info,root=-ione,&
+           & dupl=psb_dupl_add_,keeploc=.false.)
+      call tmp_ac%mv_to(tmpcoo)
+      call ac%mv_from(tmpcoo)
+
+      call psb_cdall(ctxt,desc_ac,info,mg=ntaggr,repl=.true.)
+      if (info == psb_success_) call psb_cdasb(desc_ac,info)
+
+      if (info /= psb_success_) goto 9999
+
+    case default 
       info = psb_err_internal_error_
       call psb_errpush(info,name,a_err='invalid amg_coarse_mat_')
       goto 9999
-    end if
-    
-    call op_prol%mv_to(tmpcoo)
-    nzl = tmpcoo%get_nzeros()
-    call psb_loc_to_glob(tmpcoo%ja(1:nzl),desc_ac,info,'I')
-    call tmpcoo%set_ncols(i_nr)
-    call op_prol%mv_from(tmpcoo)
-
-    call op_restr%mv_to(tmpcoo)
-    nzl = tmpcoo%get_nzeros()
-    call psb_loc_to_glob(tmpcoo%ia(1:nzl),desc_ac,info,'I')
-    call tmpcoo%set_nrows(i_nr)
-    call op_restr%mv_from(tmpcoo)
-
-
-    call psb_gather(tmp_ac,ac,desc_ac,info,root=-ione,&
-         & dupl=psb_dupl_add_,keeploc=.false.)
-    call tmp_ac%mv_to(tmpcoo)
-    call ac%mv_from(tmpcoo)
-    
-    call psb_cdall(ctxt,desc_ac,info,mg=ntaggr,repl=.true.)
-    if (info == psb_success_) call psb_cdasb(desc_ac,info)
-
-    if (info /= psb_success_) goto 9999
-
-  case default 
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='invalid amg_coarse_mat_')
-    goto 9999
-  end select
-
-
+    end select
+  else
+    call op_prol%allocate(izero,izero,info)
+    call op_restr%allocate(izero,izero,info)
+    call ac%allocate(izero,izero,info)
+  end if
   call psb_erractionrestore(err_act)
   return
 

amgprec/impl/aggregator/amg_d_dec_aggregator_mat_bld.f90

@@ -139,7 +139,7 @@ subroutine  amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   use amg_d_prec_type, amg_protect_name => amg_d_dec_aggregator_mat_bld
   use amg_d_inner_mod 
   implicit none
-  
+
   class(amg_d_dec_aggregator_type), target, intent(inout) :: ag
   type(amg_dml_parms), intent(inout)      :: parms 
   type(psb_dspmat_type), intent(in)    :: a
@@ -169,39 +169,46 @@ subroutine  amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
   ctxt = desc_a%get_context()
   call psb_info(ctxt,me,np)
 
-  !
-  ! Build the coarse-level matrix from the fine-level one, starting from 
-  ! the mapping defined by amg_aggrmap_bld and applying the aggregation
-  ! algorithm specified by 
-  !
-  select case (parms%aggr_prol)
-  case (amg_no_smooth_) 
+  if (me >=0) then 
+    !
+    ! Build the coarse-level matrix from the fine-level one, starting from 
+    ! the mapping defined by amg_aggrmap_bld and applying the aggregation
+    ! algorithm specified by 
+    !
+    select case (parms%aggr_prol)
+    case (amg_no_smooth_) 
 
-    call amg_daggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+      call amg_daggrmat_nosmth_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+           nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
 
-  case(amg_smooth_prol_,amg_l1_smooth_prol_) 
+    case(amg_smooth_prol_,amg_l1_smooth_prol_) 
 
-    call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
-            ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
-            op_restr,t_prol,info)
+      call amg_daggrmat_smth_bld(parms%aggr_prol,a,desc_a,&
+           ilaggr,nlaggr,parms,ac,desc_ac,op_prol,&
+           op_restr,t_prol,info)
 
 !!$  case(amg_biz_prol_) 
 !!$
 !!$    call amg_daggrmat_biz_bld(a,desc_a,ilaggr,nlaggr, &
 !!$         & parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
-  
-  case(amg_min_energy_) 
 
-    call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
-            nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+    case(amg_min_energy_) 
 
-  case default
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='Invalid aggr kind')
-    goto 9999
+      call amg_daggrmat_minnrg_bld(parms%aggr_prol,a,desc_a,ilaggr,&
+           nlaggr,parms,ac,desc_ac,op_prol,op_restr,t_prol,info)
+
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='Invalid aggr kind')
+      goto 9999
+
+    end select
+  else
+    call op_prol%allocate(izero,izero,info)
+    call op_restr%allocate(izero,izero,info)
+    call ac%allocate(izero,izero,info)
+  end if
 
-  end select
   if (info /= psb_success_) then
     call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner aggrmat bld')
     goto 9999
@@ -214,5 +221,5 @@ subroutine  amg_d_dec_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
 9999 call psb_error_handler(err_act)
   return
 
-  
+
 end subroutine amg_d_dec_aggregator_mat_bld

amgprec/impl/aggregator/amg_d_dec_aggregator_tprol.f90

@@ -122,19 +122,23 @@ subroutine  amg_d_dec_aggregator_build_tprol(ag,parms,ag_data,&
   call amg_check_def(parms%aggr_ord,'Ordering',&
        &   amg_aggr_ord_nat_,is_legal_ml_aggr_ord)
   call amg_check_def(parms%aggr_thresh,'Aggr_Thresh',dzero,is_legal_d_aggr_thrs)
+  if (me >=0) then 
+    !
+    ! The decoupled aggregator based on SOC measures ignores
+    ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
+    !
+    if (do_timings) call psb_tic(idx_map_bld)
+    clean_zeros = ag%do_clean_zeros
+    call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
+    if (do_timings) call psb_toc(idx_map_bld)
+    if (do_timings) call psb_tic(idx_map_tprol)
 
-  !
-  ! The decoupled aggregator based on SOC measures ignores
-  ! ag_data except for clean_zeros; soc_map_bld is a procedure pointer.
-  !
-  if (do_timings) call psb_tic(idx_map_bld)
-  clean_zeros = ag%do_clean_zeros
-  call ag%soc_map_bld(parms%aggr_ord,parms%aggr_thresh,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
-  if (do_timings) call psb_toc(idx_map_bld)
-  if (do_timings) call psb_tic(idx_map_tprol)
-
-  if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
-  if (do_timings) call psb_toc(idx_map_tprol)
+    if (info==psb_success_) call amg_map_to_tprol(desc_a,ilaggr,nlaggr,t_prol,info)
+    if (do_timings) call psb_toc(idx_map_tprol)
+  else
+    allocate(nlaggr(0),ilaggr(0))
+    call t_prol%allocate(lzero,lzero,info)
+  end if
   if (info /= psb_success_) then
     info=psb_err_from_subroutine_
     call psb_errpush(info,name,a_err='soc_map_bld/map_to_tprol')

amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90

@@ -155,7 +155,7 @@ subroutine amg_daggrmat_minnrg_bld(dol1smoothing,a,desc_a,ilaggr,nlaggr,&
   debug_level = psb_get_debug_level()
 
   ctxt = desc_a%get_context()
-  icomm = desc_a%get_mpic()
+  icomm = ctxt%get_mpic()
 
   call psb_info(ctxt, me, np)