Do not print matching statistics

Makefile

@@ -3,7 +3,7 @@ include Make.inc
 
 all:  library 
 
-library: libdir amgp
+library: libdir amgp cbnd
 #cbnd
 
 libdir:
@@ -33,8 +33,8 @@ install: all
 	mkdir -p  $(INSTALL_SAMPLESDIR) && \
 		 mkdir -p  $(INSTALL_SAMPLESDIR)/simple &&\
 	 	 mkdir -p  $(INSTALL_SAMPLESDIR)/advanced && \
-		(cd examples; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/simple ) && \
-		(cd tests; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/advanced )
+		(cd samples/simple; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/simple ) && \
+		(cd samples/advanced; /bin/cp -fr pdegen fileread $(INSTALL_SAMPLESDIR)/advanced )
 cleanlib:
 	(cd lib; /bin/rm -f *.a *$(.mod) *$(.fh))
 	(cd include; /bin/rm -f *.a *$(.mod) *$(.fh))
@@ -42,13 +42,13 @@ cleanlib:
 
 veryclean: cleanlib
 	(cd amgprec; make veryclean)
-	(cd examples/fileread; make clean)
-	(cd examples/pdegen; make clean)
-	(cd tests/fileread; make clean)
-	(cd tests/pdegen; make clean)
+	(cd samples/simple/fileread; make clean)
+	(cd samples/simple/pdegen; make clean)
+	(cd samples/advanced/fileread; make clean)
+	(cd samples/advanced/pdegen; make clean)
 
 check: all
-	make check -C tests/pdegen
+	make check -C samples/advanced/pdegen
 
 clean:
 	(cd amgprec; make clean)

amgprec/amg_base_prec_type.F90

@@ -136,6 +136,8 @@ module amg_base_prec_type
     integer(psb_lpk_)                  :: target_coarse_size
     ! 2. maximum number of levels.   Defaults to  20
     integer(psb_ipk_)                  :: max_levs    = 20_psb_ipk_
+  contains
+    procedure, pass(ag) :: default => i_ag_default
   end type amg_iaggr_data
 
   type, extends(amg_iaggr_data) :: amg_saggr_data
@@ -143,6 +145,8 @@ module amg_base_prec_type
     real(psb_spk_)                     :: min_cr_ratio   = 1.5_psb_spk_
     real(psb_spk_)                     :: op_complexity  = szero
     real(psb_spk_)                     :: avg_cr         = szero
+  contains
+    procedure, pass(ag) ::   default => s_ag_default
   end type amg_saggr_data
 
   type, extends(amg_iaggr_data) :: amg_daggr_data
@@ -150,6 +154,8 @@ module amg_base_prec_type
     real(psb_dpk_)                     :: min_cr_ratio   = 1.5_psb_dpk_
     real(psb_dpk_)                     :: op_complexity  = dzero
     real(psb_dpk_)                     :: avg_cr         = dzero
+  contains
+    procedure, pass(ag) ::   default => d_ag_default
   end type amg_daggr_data
 
 
@@ -1240,4 +1246,31 @@ contains
          & (parms1%aggr_thresh      == parms2%aggr_thresh     )
   end function amg_d_equal_aggregation
 
+  subroutine i_ag_default(ag)
+    class(amg_iaggr_data), intent(inout) :: ag
+
+    ag%min_coarse_size = -ione
+    ag%min_coarse_size_per_process = -ione
+    ag%max_levs    = 20_psb_ipk_
+  end subroutine i_ag_default
+
+  subroutine s_ag_default(ag)
+    class(amg_saggr_data), intent(inout) :: ag
+
+    call ag%amg_iaggr_data%default()
+    ag%min_cr_ratio   = 1.5_psb_spk_
+    ag%op_complexity  = szero
+    ag%avg_cr         = szero
+  end subroutine s_ag_default
+
+  subroutine d_ag_default(ag)
+    class(amg_daggr_data), intent(inout) :: ag
+
+    call ag%amg_iaggr_data%default()
+    ag%min_cr_ratio   = 1.5_psb_dpk_
+    ag%op_complexity  = dzero
+    ag%avg_cr         = dzero
+  end subroutine d_ag_default
+
+
 end module amg_base_prec_type

amgprec/amg_c_base_aggregator_mod.f90

@@ -126,7 +126,7 @@ module amg_c_base_aggregator_mod
            & psb_c_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_c_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_sml_parms), intent(inout)         :: parms 
       type(psb_c_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
@@ -144,7 +144,7 @@ module amg_c_base_aggregator_mod
            & psb_c_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_c_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_sml_parms), intent(inout)         :: parms 
       type(psb_c_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr

amgprec/amg_c_prec_type.f90

@@ -155,11 +155,12 @@ module amg_c_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_cfile_prec_descr(prec,iout,root,verbosity)
+    subroutine amg_cfile_prec_descr(prec,info,iout,root,verbosity)
       import :: amg_cprec_type, psb_ipk_
       implicit none
       ! Arguments
       class(amg_cprec_type), intent(in)      :: prec
+      integer(psb_ipk_), intent(out) :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity

amgprec/amg_d_base_aggregator_mod.f90

@@ -126,7 +126,7 @@ module amg_d_base_aggregator_mod
            & psb_d_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_d_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_dml_parms), intent(inout)         :: parms 
       type(psb_d_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
@@ -144,7 +144,7 @@ module amg_d_base_aggregator_mod
            & psb_d_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_d_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_dml_parms), intent(inout)         :: parms 
       type(psb_d_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr

amgprec/amg_d_matchboxp_mod.f90

@@ -118,6 +118,7 @@ module dmatchboxp_mod
     module procedure dPMatchBox
   end interface PMatchBox
 
+  logical, parameter, private :: print_statistics=.false.
 contains
 
   subroutine dmatchboxp_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,&
@@ -421,11 +422,13 @@ contains
         nlpairs    = v(3)
 
       end block
-      if (iam == 0) then
-        write(0,*) 'Matching statistics: Unmatched nodes ',&
-             & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
+      if (print_statistics) then
+        if (iam == 0) then
+          write(0,*) 'Matching statistics: Unmatched nodes ',&
+               & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
+        end if
       end if
-
+      
       if (display_out_) then
         block
           integer(psb_ipk_) :: idx

amgprec/amg_d_parmatch_aggregator_mod.F90

@@ -168,7 +168,7 @@ module amg_d_parmatch_aggregator_mod
       type(amg_dml_parms), intent(inout)   :: parms
       type(amg_daggr_data), intent(in)     :: ag_data
       type(psb_dspmat_type), intent(inout) :: a
-      type(psb_desc_type), intent(inout)   :: desc_a
+      type(psb_desc_type), intent(inout)     :: desc_a
       integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
       type(psb_ldspmat_type), intent(out)  :: t_prol
       integer(psb_ipk_), intent(out)      :: info
@@ -235,7 +235,7 @@ module amg_d_parmatch_aggregator_mod
            & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data
       implicit none
       type(psb_dspmat_type), intent(in)       :: a
-      type(psb_desc_type), intent(in)         :: desc_a
+      type(psb_desc_type), intent(inout)       :: desc_a
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_dml_parms), intent(inout)      :: parms
       type(psb_ldspmat_type), intent(inout)   :: t_prol
@@ -257,7 +257,7 @@ module amg_d_parmatch_aggregator_mod
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_dml_parms), intent(inout)      :: parms
       type(psb_ldspmat_type), intent(inout)   :: t_prol
-      type(psb_dspmat_type), intent(out)      :: op_prol,ac, op_restr
+      type(psb_dspmat_type), intent(inout)      :: op_prol,ac, op_restr
       type(psb_desc_type), intent(inout)    :: desc_ac
       integer(psb_ipk_), intent(out)          :: info
     end subroutine amg_d_parmatch_unsmth_bld
@@ -275,7 +275,7 @@ module amg_d_parmatch_aggregator_mod
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_dml_parms), intent(inout)      :: parms
       type(psb_ldspmat_type), intent(inout)   :: t_prol
-      type(psb_dspmat_type), intent(out)      :: op_prol,ac, op_restr
+      type(psb_dspmat_type), intent(inout)    :: op_prol,ac, op_restr
       type(psb_desc_type), intent(inout)    :: desc_ac
       integer(psb_ipk_), intent(out)          :: info
     end subroutine amg_d_parmatch_smth_bld
@@ -288,11 +288,11 @@ module amg_d_parmatch_aggregator_mod
            & psb_ldspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_dml_parms, amg_daggr_data
       implicit none
       type(psb_dspmat_type), intent(inout)    :: a
-      type(psb_desc_type), intent(in)         :: desc_a
+      type(psb_desc_type), intent(inout)       :: desc_a
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_dml_parms), intent(inout)      :: parms
       type(psb_ldspmat_type), intent(inout)   :: t_prol
-      type(psb_dspmat_type), intent(out)      :: op_prol,ac, op_restr
+      type(psb_dspmat_type), intent(inout)    :: op_prol,ac, op_restr
       type(psb_desc_type), intent(out)          :: desc_ac
       integer(psb_ipk_), intent(out)          :: info
     end subroutine amg_d_parmatch_spmm_bld_ov
@@ -306,11 +306,11 @@ module amg_d_parmatch_aggregator_mod
            & psb_d_csr_sparse_mat, psb_ld_csr_sparse_mat
       implicit none
       type(psb_d_csr_sparse_mat), intent(inout) :: a
-      type(psb_desc_type), intent(in)           :: desc_a
+      type(psb_desc_type), intent(inout)          :: desc_a
       integer(psb_lpk_), intent(inout)          :: ilaggr(:), nlaggr(:)
       type(amg_dml_parms), intent(inout)        :: parms
       type(psb_ldspmat_type), intent(inout)     :: t_prol
-      type(psb_dspmat_type), intent(out)        :: op_prol,ac, op_restr
+      type(psb_dspmat_type), intent(inout)      :: op_prol,ac, op_restr
       type(psb_desc_type), intent(out)          :: desc_ac
       integer(psb_ipk_), intent(out)            :: info
     end subroutine amg_d_parmatch_spmm_bld_inner

amgprec/amg_d_prec_type.f90

@@ -155,11 +155,12 @@ module amg_d_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_dfile_prec_descr(prec,iout,root,verbosity)
+    subroutine amg_dfile_prec_descr(prec,info,iout,root,verbosity)
       import :: amg_dprec_type, psb_ipk_
       implicit none
       ! Arguments
       class(amg_dprec_type), intent(in)      :: prec
+      integer(psb_ipk_), intent(out) :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity

amgprec/amg_s_base_aggregator_mod.f90

@@ -126,7 +126,7 @@ module amg_s_base_aggregator_mod
            & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_s_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_sml_parms), intent(inout)         :: parms 
       type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
@@ -144,7 +144,7 @@ module amg_s_base_aggregator_mod
            & psb_s_coo_sparse_mat, amg_sml_parms, psb_spk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_s_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_sml_parms), intent(inout)         :: parms 
       type(psb_s_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr

amgprec/amg_s_matchboxp_mod.f90

@@ -118,6 +118,7 @@ module smatchboxp_mod
     module procedure sPMatchBox
   end interface PMatchBox
 
+  logical, parameter, private :: print_statistics=.false.
 contains
 
   subroutine smatchboxp_build_prol(w,a,desc_a,ilaggr,nlaggr,prol,info,&
@@ -421,11 +422,13 @@ contains
         nlpairs    = v(3)
 
       end block
-      if (iam == 0) then
-        write(0,*) 'Matching statistics: Unmatched nodes ',&
-             & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
+      if (print_statistics) then
+        if (iam == 0) then
+          write(0,*) 'Matching statistics: Unmatched nodes ',&
+               & nunmatched,' Singletons:',nlsingl,' Pairs:',nlpairs
+        end if
       end if
-
+      
       if (display_out_) then
         block
           integer(psb_ipk_) :: idx

amgprec/amg_s_parmatch_aggregator_mod.F90

@@ -168,7 +168,7 @@ module amg_s_parmatch_aggregator_mod
       type(amg_sml_parms), intent(inout)   :: parms
       type(amg_saggr_data), intent(in)     :: ag_data
       type(psb_sspmat_type), intent(inout) :: a
-      type(psb_desc_type), intent(inout)   :: desc_a
+      type(psb_desc_type), intent(inout)     :: desc_a
       integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
       type(psb_lsspmat_type), intent(out)  :: t_prol
       integer(psb_ipk_), intent(out)      :: info
@@ -235,7 +235,7 @@ module amg_s_parmatch_aggregator_mod
            & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data
       implicit none
       type(psb_sspmat_type), intent(in)       :: a
-      type(psb_desc_type), intent(in)         :: desc_a
+      type(psb_desc_type), intent(inout)       :: desc_a
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_sml_parms), intent(inout)      :: parms
       type(psb_lsspmat_type), intent(inout)   :: t_prol
@@ -257,7 +257,7 @@ module amg_s_parmatch_aggregator_mod
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_sml_parms), intent(inout)      :: parms
       type(psb_lsspmat_type), intent(inout)   :: t_prol
-      type(psb_sspmat_type), intent(out)      :: op_prol,ac, op_restr
+      type(psb_sspmat_type), intent(inout)      :: op_prol,ac, op_restr
       type(psb_desc_type), intent(inout)    :: desc_ac
       integer(psb_ipk_), intent(out)          :: info
     end subroutine amg_s_parmatch_unsmth_bld
@@ -275,7 +275,7 @@ module amg_s_parmatch_aggregator_mod
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_sml_parms), intent(inout)      :: parms
       type(psb_lsspmat_type), intent(inout)   :: t_prol
-      type(psb_sspmat_type), intent(out)      :: op_prol,ac, op_restr
+      type(psb_sspmat_type), intent(inout)    :: op_prol,ac, op_restr
       type(psb_desc_type), intent(inout)    :: desc_ac
       integer(psb_ipk_), intent(out)          :: info
     end subroutine amg_s_parmatch_smth_bld
@@ -288,11 +288,11 @@ module amg_s_parmatch_aggregator_mod
            & psb_lsspmat_type, psb_dpk_, psb_ipk_, psb_lpk_, amg_sml_parms, amg_saggr_data
       implicit none
       type(psb_sspmat_type), intent(inout)    :: a
-      type(psb_desc_type), intent(in)         :: desc_a
+      type(psb_desc_type), intent(inout)       :: desc_a
       integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
       type(amg_sml_parms), intent(inout)      :: parms
       type(psb_lsspmat_type), intent(inout)   :: t_prol
-      type(psb_sspmat_type), intent(out)      :: op_prol,ac, op_restr
+      type(psb_sspmat_type), intent(inout)    :: op_prol,ac, op_restr
       type(psb_desc_type), intent(out)          :: desc_ac
       integer(psb_ipk_), intent(out)          :: info
     end subroutine amg_s_parmatch_spmm_bld_ov
@@ -306,11 +306,11 @@ module amg_s_parmatch_aggregator_mod
            & psb_s_csr_sparse_mat, psb_ls_csr_sparse_mat
       implicit none
       type(psb_s_csr_sparse_mat), intent(inout) :: a
-      type(psb_desc_type), intent(in)           :: desc_a
+      type(psb_desc_type), intent(inout)          :: desc_a
       integer(psb_lpk_), intent(inout)          :: ilaggr(:), nlaggr(:)
       type(amg_sml_parms), intent(inout)        :: parms
       type(psb_lsspmat_type), intent(inout)     :: t_prol
-      type(psb_sspmat_type), intent(out)        :: op_prol,ac, op_restr
+      type(psb_sspmat_type), intent(inout)      :: op_prol,ac, op_restr
       type(psb_desc_type), intent(out)          :: desc_ac
       integer(psb_ipk_), intent(out)            :: info
     end subroutine amg_s_parmatch_spmm_bld_inner

amgprec/amg_s_prec_type.f90

@@ -155,11 +155,12 @@ module amg_s_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_sfile_prec_descr(prec,iout,root,verbosity)
+    subroutine amg_sfile_prec_descr(prec,info,iout,root,verbosity)
       import :: amg_sprec_type, psb_ipk_
       implicit none
       ! Arguments
       class(amg_sprec_type), intent(in)      :: prec
+      integer(psb_ipk_), intent(out) :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity

amgprec/amg_z_base_aggregator_mod.f90

@@ -126,7 +126,7 @@ module amg_z_base_aggregator_mod
            & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_z_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_dml_parms), intent(inout)         :: parms 
       type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
@@ -144,7 +144,7 @@ module amg_z_base_aggregator_mod
            & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
       implicit none
       type(psb_z_csr_sparse_mat), intent(inout) :: a_csr
-      type(psb_desc_type), intent(in)            :: desc_a
+      type(psb_desc_type), intent(inout)         :: desc_a
       integer(psb_lpk_), intent(inout)           :: nlaggr(:)
       type(amg_dml_parms), intent(inout)         :: parms 
       type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr

amgprec/amg_z_prec_type.f90

@@ -155,11 +155,12 @@ module amg_z_prec_type
 
 
   interface amg_precdescr
-    subroutine amg_zfile_prec_descr(prec,iout,root,verbosity)
+    subroutine amg_zfile_prec_descr(prec,info,iout,root,verbosity)
       import :: amg_zprec_type, psb_ipk_
       implicit none
       ! Arguments
       class(amg_zprec_type), intent(in)      :: prec
+      integer(psb_ipk_), intent(out) :: info
       integer(psb_ipk_), intent(in), optional :: iout
       integer(psb_ipk_), intent(in), optional :: root
       integer(psb_ipk_), intent(in), optional :: verbosity

amgprec/impl/aggregator/amg_c_dec_aggregator_tprol.f90

@@ -83,8 +83,8 @@ subroutine  amg_c_dec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_c_dec_aggregator_type), target, intent(inout) :: ag
   type(amg_sml_parms), intent(inout)  :: parms
   type(amg_saggr_data), intent(in)    :: ag_data
-  type(psb_cspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_cspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)    :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lcspmat_type), intent(out)  :: t_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_c_symdec_aggregator_tprol.f90

@@ -86,8 +86,8 @@ subroutine  amg_c_symdec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_c_symdec_aggregator_type), target, intent(inout) :: ag
   type(amg_sml_parms), intent(inout)  :: parms
   type(amg_saggr_data), intent(in)    :: ag_data
-  type(psb_cspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_cspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)   :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lcspmat_type), intent(out)  :: op_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_caggrmat_minnrg_bld.f90

@@ -105,7 +105,7 @@
 !
 !
 subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,info)
+     & ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_c_inner_mod, amg_protect_name => amg_caggrmat_minnrg_bld
@@ -117,8 +117,8 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)         :: parms 
-  type(psb_lcspmat_type), intent(inout)       :: op_prol
-  type(psb_lcspmat_type), intent(out)         :: ac,op_restr
+  type(psb_lcspmat_type), intent(inout)       :: t_prol
+  type(psb_cspmat_type), intent(inout)       :: op_prol, ac,op_restr
   type(psb_desc_type), intent(inout)            :: desc_ac
   integer(psb_ipk_), intent(out)                :: info
 
@@ -171,6 +171,8 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  !NEEDS TO BE REWORKED !!
+  
   ! naggr: number of local aggregates
   ! nrow: local rows. 
   ! 
@@ -183,361 +185,361 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
     goto 9999      
   end if
 
-  ! Get the diagonal D
-  adiag =  a%get_diag(info)
-  if (info == psb_success_) &
-       & call psb_realloc(ncol,adiag,info)    
-  if (info == psb_success_) &
-       & call psb_halo(adiag,desc_a,info)
-  if (info == psb_success_) call a%cp_to_l(la)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
-    goto 9999
-  end if
-
-  do i=1,size(adiag)
-    if (adiag(i) /= czero) then
-      adinv(i) = cone / adiag(i)
-    else
-      adinv(i) = cone
-    end if
-  end do
-
-
-
-  ! 1. Allocate Ptilde in sparse matrix form 
-  call op_prol%mv_to(tmpcoo)
-  call ptilde%mv_from(tmpcoo)
-  call ptilde%cscnv(info,type='csr')
-
-  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
-  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & ' Initial copies done.'
-
-  call da%scal(adinv,info)
-
-  call psb_spspmm(da,ptilde,dap,info)
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
-    goto 9999
-  end if
-
-  call dap%clone(atmp,info)
-
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  call psb_spspmm(da,atmp,dadap,info)
-  call atmp%free()
-
-  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
-  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
-  call dap%mv_to(csc_dap)
-  call dadap%mv_to(csc_dadap)
-
-  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
-  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-  ! !$  write(0,*) trim(name),' OMP :',omp
-  ! !$  write(0,*) trim(name),' ODEN:',oden
-
-  omp = omp/oden
-
-  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done NUMBMM 1'
-
-  call am3%mv_to(acsr3)
-  ! Compute omega_int
-  ommx = czero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = czero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  do i=1, nrow
-    omf(i) = ommx
-    do j=acsr3%irp(i),acsr3%irp(i+1)-1
-      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero
-    if(psb_minreal(omf(i)) < szero) omf(i) = czero
-  end do
-
-  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
-
-  if (filter_mat) then
-    !
-    ! Build the filtered matrix Af from A
-    ! 
-    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
-
-    do i=1,nrow
-      tmp = czero
-      jd  = -1 
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) jd = j 
-        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
-          tmp=tmp+acsrf%val(j)
-          acsrf%val(j)=czero
-        endif
-      enddo
-      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
-        acsrf%val(jd)=acsrf%val(jd)-tmp
-      end if
-    enddo
-    ! Take out zeroed terms 
-    call acsrf%clean_zeros(info)
-
-    !
-    ! Build the smoothed prolongator using the filtered matrix
-    ! 
-    do i=1,acsrf%get_nrows()
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) then 
-          acsrf%val(j) = cone - omf(i)*acsrf%val(j) 
-        else
-          acsrf%val(j) = - omf(i)*acsrf%val(j) 
-        end if
-      end do
-    end do
-
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-
-    call af%mv_from(acsrf)
-    !
-    ! op_prol = (I-w*D*Af)Ptilde
-    ! Doing it this way means to consider diag(Af_i)
-    ! 
-    !
-    call psb_spspmm(af,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done SPSPMM 1'
-  else
-    !
-    ! Build the smoothed prolongator using the original matrix
-    !
-    do i=1,acsr3%get_nrows()
-      do j=acsr3%irp(i),acsr3%irp(i+1)-1
-        if (acsr3%ja(j) == i) then 
-          acsr3%val(j) = cone - omf(i)*acsr3%val(j) 
-        else
-          acsr3%val(j) = - omf(i)*acsr3%val(j) 
-        end if
-      end do
-    end do
-
-    call am3%mv_from(acsr3)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-    !
-    ! 
-    ! op_prol = (I-w*D*A)Ptilde
-    ! 
-    !
-    call psb_spspmm(am3,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done NUMBMM 1'
-
-  end if
-
-
-  !
-  ! Ok, let's start over with the restrictor
-  ! 
-  call ptilde%transc(rtilde)
-  call la%cscnv(atmp,info,type='csr')
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.true.,rowscale=.true.)
-  nrt  = am4%get_nrows() 
-  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
-  call atmp2%cscnv(info,type='CSR')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
-  call am4%free()
-  call atmp2%free()
-
-  ! This is to compute the transpose. It ONLY works if the
-  ! original A has a symmetric pattern.
-  call atmp%transc(atmp2) 
-  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
-  call dat%cscnv(info,type='csr')
-  call dat%scal(adinv,info)
-
-  ! Now for the product. 
-  call psb_spspmm(dat,ptilde,datp,info)
-
-  call datp%clone(atmp2,info)
-  call psb_sphalo(atmp2,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-
-  call psb_symbmm(dat,atmp2,datdatp,info)
-  call psb_numbmm(dat,atmp2,datdatp)
-  call atmp2%free()
-
-  call datp%mv_to(csc_datp)    
-  call datdatp%mv_to(csc_datdatp)    
-
-  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
-  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-
-
-  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
-  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
-  omp = omp/oden
-  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
-  ! Compute omega_int
-  ommx = czero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = czero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  ! Going over the columns of atmp means going over the rows
-  ! of A^T. Hopefully ;-) 
-  call atmp%cp_to(acsc)
-
-  do i=1, nrow
-    omf(i) = ommx
-    do j= acsc%icp(i),acsc%icp(i+1)-1
-      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero
-    if(psb_minreal(omf(i)) < szero) omf(i) = czero
-  end do
-  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
-  call psb_halo(omf,desc_a,info)
-  call acsc%free() 
-
-
-  call atmp%mv_to(acsr1)
-
-  do i=1,acsr1%get_nrows()
-    do j=acsr1%irp(i),acsr1%irp(i+1)-1
-      if (acsr1%ja(j) == i) then 
-        acsr1%val(j) = cone - acsr1%val(j)*omf(acsr1%ja(j))
-      else
-        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
-      end if
-    end do
-  end do
-  call atmp%mv_from(acsr1)
-
-  call rtilde%mv_to(tmpcoo)
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call rtilde%mv_from(tmpcoo)
-  call rtilde%cscnv(info,type='csr')
-
-  call psb_spspmm(rtilde,atmp,op_restr,info)
-
-  !
-  ! Now we have to gather the halo of op_prol, and add it to itself
-  ! to multiply it by A,
-  !
-  call op_prol%clone(tmp_prol,info) 
-  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
-    goto 9999
-  end if
-
-  !
-  ! Now we have to fix this.  The only rows of B that are correct 
-  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
-  !
-  call op_restr%mv_to(tmpcoo)
-
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call op_restr%mv_from(tmpcoo)
-  call op_restr%cscnv(info,type='csr')
-
-
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'starting sphalo/ rwxtd'
-
-  call psb_spspmm(la,tmp_prol,am3,info)
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done SPSPMM 2'
-
-  call psb_sphalo(am3,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         & a_err='Extend am3')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done sphalo/ rwxtd'
-
-  call psb_spspmm(op_restr,am3,ac,info)
-  if (info == psb_success_) call am3%free()
-  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
-
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         &a_err='Build ac = op_restr x am3')
-    goto 9999
-  end if
+!!$  ! Get the diagonal D
+!!$  adiag =  a%get_diag(info)
+!!$  if (info == psb_success_) &
+!!$       & call psb_realloc(ncol,adiag,info)    
+!!$  if (info == psb_success_) &
+!!$       & call psb_halo(adiag,desc_a,info)
+!!$  if (info == psb_success_) call a%cp_to_l(la)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  do i=1,size(adiag)
+!!$    if (adiag(i) /= czero) then
+!!$      adinv(i) = cone / adiag(i)
+!!$    else
+!!$      adinv(i) = cone
+!!$    end if
+!!$  end do
+!!$
+!!$
+!!$
+!!$  ! 1. Allocate Ptilde in sparse matrix form 
+!!$  call op_prol%mv_to(tmpcoo)
+!!$  call ptilde%mv_from(tmpcoo)
+!!$  call ptilde%cscnv(info,type='csr')
+!!$
+!!$  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & ' Initial copies done.'
+!!$
+!!$  call da%scal(adinv,info)
+!!$
+!!$  call psb_spspmm(da,ptilde,dap,info)
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  call dap%clone(atmp,info)
+!!$
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  call psb_spspmm(da,atmp,dadap,info)
+!!$  call atmp%free()
+!!$
+!!$  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
+!!$  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
+!!$  call dap%mv_to(csc_dap)
+!!$  call dadap%mv_to(csc_dadap)
+!!$
+!!$  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
+!!$  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$  ! !$  write(0,*) trim(name),' OMP :',omp
+!!$  ! !$  write(0,*) trim(name),' ODEN:',oden
+!!$
+!!$  omp = omp/oden
+!!$
+!!$  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done NUMBMM 1'
+!!$
+!!$  call am3%mv_to(acsr3)
+!!$  ! Compute omega_int
+!!$  ommx = czero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = czero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
+!!$    end do
+!!$ ! !     if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero
+!!$    if(psb_minreal(omf(i)) < szero) omf(i) = czero
+!!$  end do
+!!$
+!!$  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
+!!$
+!!$  if (filter_mat) then
+!!$    !
+!!$    ! Build the filtered matrix Af from A
+!!$    ! 
+!!$    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
+!!$
+!!$    do i=1,nrow
+!!$      tmp = czero
+!!$      jd  = -1 
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) jd = j 
+!!$        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
+!!$          tmp=tmp+acsrf%val(j)
+!!$          acsrf%val(j)=czero
+!!$        endif
+!!$      enddo
+!!$      if (jd == -1) then 
+!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$      else
+!!$        acsrf%val(jd)=acsrf%val(jd)-tmp
+!!$      end if
+!!$    enddo
+!!$    ! Take out zeroed terms 
+!!$    call acsrf%clean_zeros(info)
+!!$
+!!$    !
+!!$    ! Build the smoothed prolongator using the filtered matrix
+!!$    ! 
+!!$    do i=1,acsrf%get_nrows()
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) then 
+!!$          acsrf%val(j) = cone - omf(i)*acsrf%val(j) 
+!!$        else
+!!$          acsrf%val(j) = - omf(i)*acsrf%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$
+!!$    call af%mv_from(acsrf)
+!!$    !
+!!$    ! op_prol = (I-w*D*Af)Ptilde
+!!$    ! Doing it this way means to consider diag(Af_i)
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(af,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done SPSPMM 1'
+!!$  else
+!!$    !
+!!$    ! Build the smoothed prolongator using the original matrix
+!!$    !
+!!$    do i=1,acsr3%get_nrows()
+!!$      do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$        if (acsr3%ja(j) == i) then 
+!!$          acsr3%val(j) = cone - omf(i)*acsr3%val(j) 
+!!$        else
+!!$          acsr3%val(j) = - omf(i)*acsr3%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    call am3%mv_from(acsr3)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$    !
+!!$    ! 
+!!$    ! op_prol = (I-w*D*A)Ptilde
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(am3,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done NUMBMM 1'
+!!$
+!!$  end if
+!!$
+!!$
+!!$  !
+!!$  ! Ok, let's start over with the restrictor
+!!$  ! 
+!!$  call ptilde%transc(rtilde)
+!!$  call la%cscnv(atmp,info,type='csr')
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.true.,rowscale=.true.)
+!!$  nrt  = am4%get_nrows() 
+!!$  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
+!!$  call atmp2%cscnv(info,type='CSR')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
+!!$  call am4%free()
+!!$  call atmp2%free()
+!!$
+!!$  ! This is to compute the transpose. It ONLY works if the
+!!$  ! original A has a symmetric pattern.
+!!$  call atmp%transc(atmp2) 
+!!$  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
+!!$  call dat%cscnv(info,type='csr')
+!!$  call dat%scal(adinv,info)
+!!$
+!!$  ! Now for the product. 
+!!$  call psb_spspmm(dat,ptilde,datp,info)
+!!$
+!!$  call datp%clone(atmp2,info)
+!!$  call psb_sphalo(atmp2,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$
+!!$  call psb_symbmm(dat,atmp2,datdatp,info)
+!!$  call psb_numbmm(dat,atmp2,datdatp)
+!!$  call atmp2%free()
+!!$
+!!$  call datp%mv_to(csc_datp)    
+!!$  call datdatp%mv_to(csc_datdatp)    
+!!$
+!!$  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
+!!$  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$
+!!$
+!!$  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
+!!$  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
+!!$  omp = omp/oden
+!!$  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
+!!$  ! Compute omega_int
+!!$  ommx = czero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = czero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  ! Going over the columns of atmp means going over the rows
+!!$  ! of A^T. Hopefully ;-) 
+!!$  call atmp%cp_to(acsc)
+!!$
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j= acsc%icp(i),acsc%icp(i+1)-1
+!!$      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
+!!$    end do
+!!$ ! !    if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = czero
+!!$    if(psb_minreal(omf(i)) < szero) omf(i) = czero
+!!$  end do
+!!$  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
+!!$  call psb_halo(omf,desc_a,info)
+!!$  call acsc%free() 
+!!$
+!!$
+!!$  call atmp%mv_to(acsr1)
+!!$
+!!$  do i=1,acsr1%get_nrows()
+!!$    do j=acsr1%irp(i),acsr1%irp(i+1)-1
+!!$      if (acsr1%ja(j) == i) then 
+!!$        acsr1%val(j) = cone - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      else
+!!$        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      end if
+!!$    end do
+!!$  end do
+!!$  call atmp%mv_from(acsr1)
+!!$
+!!$  call rtilde%mv_to(tmpcoo)
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call rtilde%mv_from(tmpcoo)
+!!$  call rtilde%cscnv(info,type='csr')
+!!$
+!!$  call psb_spspmm(rtilde,atmp,op_restr,info)
+!!$
+!!$  !
+!!$  ! Now we have to gather the halo of op_prol, and add it to itself
+!!$  ! to multiply it by A,
+!!$  !
+!!$  call op_prol%clone(tmp_prol,info) 
+!!$  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  !
+!!$  ! Now we have to fix this.  The only rows of B that are correct 
+!!$  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
+!!$  !
+!!$  call op_restr%mv_to(tmpcoo)
+!!$
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call op_restr%mv_from(tmpcoo)
+!!$  call op_restr%cscnv(info,type='csr')
+!!$
+!!$
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'starting sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(la,tmp_prol,am3,info)
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done SPSPMM 2'
+!!$
+!!$  call psb_sphalo(am3,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         & a_err='Extend am3')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(op_restr,am3,ac,info)
+!!$  if (info == psb_success_) call am3%free()
+!!$  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
+!!$
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         &a_err='Build ac = op_restr x am3')
+!!$    goto 9999
+!!$  end if
 
 
 

amgprec/impl/aggregator/amg_caggrmat_smth_bld.f90

@@ -116,7 +116,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)      :: parms 
-  type(psb_cspmat_type), intent(out)      :: op_prol,ac,op_restr
+  type(psb_cspmat_type), intent(inout)      :: op_prol,ac,op_restr
   type(psb_lcspmat_type), intent(inout)    :: t_prol
   type(psb_desc_type), intent(inout)       :: desc_ac
   integer(psb_ipk_), intent(out)           :: info

amgprec/impl/aggregator/amg_d_dec_aggregator_tprol.f90

@@ -83,8 +83,8 @@ subroutine  amg_d_dec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_d_dec_aggregator_type), target, intent(inout) :: ag
   type(amg_dml_parms), intent(inout)  :: parms
   type(amg_daggr_data), intent(in)    :: ag_data
-  type(psb_dspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_dspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)    :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_ldspmat_type), intent(out)  :: t_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_d_parmatch_aggregator_tprol.F90

@@ -58,7 +58,7 @@ subroutine  amg_d_parmatch_aggregator_build_tprol(ag,parms,ag_data,&
   type(amg_dml_parms), intent(inout)   :: parms
   type(amg_daggr_data), intent(in)     :: ag_data
   type(psb_dspmat_type), intent(inout) :: a
-  type(psb_desc_type), intent(inout)   :: desc_a
+  type(psb_desc_type), intent(inout)      :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_ldspmat_type), intent(out)  :: t_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_d_parmatch_smth_bld.F90

@@ -122,7 +122,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)    :: parms
   type(psb_ldspmat_type), intent(inout)  :: t_prol
-  type(psb_dspmat_type), intent(out)    :: op_prol,ac,op_restr
+  type(psb_dspmat_type), intent(inout)   :: op_prol,ac,op_restr
   type(psb_desc_type), intent(inout)    :: desc_ac
   integer(psb_ipk_), intent(out)           :: info
 

amgprec/impl/aggregator/amg_d_parmatch_spmm_bld.F90

@@ -108,7 +108,7 @@ subroutine amg_d_parmatch_spmm_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   ! Arguments
   type(psb_dspmat_type), intent(in)       :: a
-  type(psb_desc_type), intent(in)         :: desc_a
+  type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)      :: parms
   type(psb_ldspmat_type), intent(inout)   :: t_prol

amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_inner.F90

@@ -108,11 +108,11 @@ subroutine amg_d_parmatch_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,&
 
   ! Arguments
   type(psb_d_csr_sparse_mat), intent(inout) :: a_csr
-  type(psb_desc_type), intent(in)           :: desc_a
+  type(psb_desc_type), intent(inout)          :: desc_a
   integer(psb_lpk_), intent(inout)          :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)        :: parms
   type(psb_ldspmat_type), intent(inout)     :: t_prol
-  type(psb_dspmat_type), intent(out)        :: ac, op_prol, op_restr
+  type(psb_dspmat_type), intent(inout)      :: ac, op_prol, op_restr
   type(psb_desc_type), intent(out)          :: desc_ac
   integer(psb_ipk_), intent(out)            :: info
 

amgprec/impl/aggregator/amg_d_parmatch_spmm_bld_ov.F90

@@ -108,7 +108,7 @@ subroutine amg_d_parmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,&
 
   ! Arguments
   type(psb_dspmat_type), intent(inout)    :: a
-  type(psb_desc_type), intent(in)         :: desc_a
+  type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)      :: parms
   type(psb_ldspmat_type), intent(inout)   :: t_prol

amgprec/impl/aggregator/amg_d_symdec_aggregator_tprol.f90

@@ -86,8 +86,8 @@ subroutine  amg_d_symdec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_d_symdec_aggregator_type), target, intent(inout) :: ag
   type(amg_dml_parms), intent(inout)  :: parms
   type(amg_daggr_data), intent(in)    :: ag_data
-  type(psb_dspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_dspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)   :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_ldspmat_type), intent(out)  :: op_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_daggrmat_minnrg_bld.f90

@@ -105,7 +105,7 @@
 !
 !
 subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,info)
+     & ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_d_inner_mod, amg_protect_name => amg_daggrmat_minnrg_bld
@@ -117,8 +117,8 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)         :: parms 
-  type(psb_ldspmat_type), intent(inout)       :: op_prol
-  type(psb_ldspmat_type), intent(out)         :: ac,op_restr
+  type(psb_ldspmat_type), intent(inout)       :: t_prol
+  type(psb_dspmat_type), intent(inout)       :: op_prol, ac,op_restr
   type(psb_desc_type), intent(inout)            :: desc_ac
   integer(psb_ipk_), intent(out)                :: info
 
@@ -171,6 +171,8 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  !NEEDS TO BE REWORKED !!
+  
   ! naggr: number of local aggregates
   ! nrow: local rows. 
   ! 
@@ -183,361 +185,361 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
     goto 9999      
   end if
 
-  ! Get the diagonal D
-  adiag =  a%get_diag(info)
-  if (info == psb_success_) &
-       & call psb_realloc(ncol,adiag,info)    
-  if (info == psb_success_) &
-       & call psb_halo(adiag,desc_a,info)
-  if (info == psb_success_) call a%cp_to_l(la)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
-    goto 9999
-  end if
-
-  do i=1,size(adiag)
-    if (adiag(i) /= dzero) then
-      adinv(i) = done / adiag(i)
-    else
-      adinv(i) = done
-    end if
-  end do
-
-
-
-  ! 1. Allocate Ptilde in sparse matrix form 
-  call op_prol%mv_to(tmpcoo)
-  call ptilde%mv_from(tmpcoo)
-  call ptilde%cscnv(info,type='csr')
-
-  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
-  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & ' Initial copies done.'
-
-  call da%scal(adinv,info)
-
-  call psb_spspmm(da,ptilde,dap,info)
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
-    goto 9999
-  end if
-
-  call dap%clone(atmp,info)
-
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  call psb_spspmm(da,atmp,dadap,info)
-  call atmp%free()
-
-  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
-  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
-  call dap%mv_to(csc_dap)
-  call dadap%mv_to(csc_dadap)
-
-  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
-  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-  ! !$  write(0,*) trim(name),' OMP :',omp
-  ! !$  write(0,*) trim(name),' ODEN:',oden
-
-  omp = omp/oden
-
-  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done NUMBMM 1'
-
-  call am3%mv_to(acsr3)
-  ! Compute omega_int
-  ommx = dzero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = dzero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  do i=1, nrow
-    omf(i) = ommx
-    do j=acsr3%irp(i),acsr3%irp(i+1)-1
-      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero
-    if(psb_minreal(omf(i)) < dzero) omf(i) = dzero
-  end do
-
-  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
-
-  if (filter_mat) then
-    !
-    ! Build the filtered matrix Af from A
-    ! 
-    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
-
-    do i=1,nrow
-      tmp = dzero
-      jd  = -1 
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) jd = j 
-        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
-          tmp=tmp+acsrf%val(j)
-          acsrf%val(j)=dzero
-        endif
-      enddo
-      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
-        acsrf%val(jd)=acsrf%val(jd)-tmp
-      end if
-    enddo
-    ! Take out zeroed terms 
-    call acsrf%clean_zeros(info)
-
-    !
-    ! Build the smoothed prolongator using the filtered matrix
-    ! 
-    do i=1,acsrf%get_nrows()
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) then 
-          acsrf%val(j) = done - omf(i)*acsrf%val(j) 
-        else
-          acsrf%val(j) = - omf(i)*acsrf%val(j) 
-        end if
-      end do
-    end do
-
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-
-    call af%mv_from(acsrf)
-    !
-    ! op_prol = (I-w*D*Af)Ptilde
-    ! Doing it this way means to consider diag(Af_i)
-    ! 
-    !
-    call psb_spspmm(af,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done SPSPMM 1'
-  else
-    !
-    ! Build the smoothed prolongator using the original matrix
-    !
-    do i=1,acsr3%get_nrows()
-      do j=acsr3%irp(i),acsr3%irp(i+1)-1
-        if (acsr3%ja(j) == i) then 
-          acsr3%val(j) = done - omf(i)*acsr3%val(j) 
-        else
-          acsr3%val(j) = - omf(i)*acsr3%val(j) 
-        end if
-      end do
-    end do
-
-    call am3%mv_from(acsr3)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-    !
-    ! 
-    ! op_prol = (I-w*D*A)Ptilde
-    ! 
-    !
-    call psb_spspmm(am3,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done NUMBMM 1'
-
-  end if
-
-
-  !
-  ! Ok, let's start over with the restrictor
-  ! 
-  call ptilde%transc(rtilde)
-  call la%cscnv(atmp,info,type='csr')
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.true.,rowscale=.true.)
-  nrt  = am4%get_nrows() 
-  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
-  call atmp2%cscnv(info,type='CSR')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
-  call am4%free()
-  call atmp2%free()
-
-  ! This is to compute the transpose. It ONLY works if the
-  ! original A has a symmetric pattern.
-  call atmp%transc(atmp2) 
-  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
-  call dat%cscnv(info,type='csr')
-  call dat%scal(adinv,info)
-
-  ! Now for the product. 
-  call psb_spspmm(dat,ptilde,datp,info)
-
-  call datp%clone(atmp2,info)
-  call psb_sphalo(atmp2,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-
-  call psb_symbmm(dat,atmp2,datdatp,info)
-  call psb_numbmm(dat,atmp2,datdatp)
-  call atmp2%free()
-
-  call datp%mv_to(csc_datp)    
-  call datdatp%mv_to(csc_datdatp)    
-
-  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
-  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-
-
-  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
-  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
-  omp = omp/oden
-  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
-  ! Compute omega_int
-  ommx = dzero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = dzero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  ! Going over the columns of atmp means going over the rows
-  ! of A^T. Hopefully ;-) 
-  call atmp%cp_to(acsc)
-
-  do i=1, nrow
-    omf(i) = ommx
-    do j= acsc%icp(i),acsc%icp(i+1)-1
-      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero
-    if(psb_minreal(omf(i)) < dzero) omf(i) = dzero
-  end do
-  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
-  call psb_halo(omf,desc_a,info)
-  call acsc%free() 
-
-
-  call atmp%mv_to(acsr1)
-
-  do i=1,acsr1%get_nrows()
-    do j=acsr1%irp(i),acsr1%irp(i+1)-1
-      if (acsr1%ja(j) == i) then 
-        acsr1%val(j) = done - acsr1%val(j)*omf(acsr1%ja(j))
-      else
-        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
-      end if
-    end do
-  end do
-  call atmp%mv_from(acsr1)
-
-  call rtilde%mv_to(tmpcoo)
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call rtilde%mv_from(tmpcoo)
-  call rtilde%cscnv(info,type='csr')
-
-  call psb_spspmm(rtilde,atmp,op_restr,info)
-
-  !
-  ! Now we have to gather the halo of op_prol, and add it to itself
-  ! to multiply it by A,
-  !
-  call op_prol%clone(tmp_prol,info) 
-  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
-    goto 9999
-  end if
-
-  !
-  ! Now we have to fix this.  The only rows of B that are correct 
-  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
-  !
-  call op_restr%mv_to(tmpcoo)
-
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call op_restr%mv_from(tmpcoo)
-  call op_restr%cscnv(info,type='csr')
-
-
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'starting sphalo/ rwxtd'
-
-  call psb_spspmm(la,tmp_prol,am3,info)
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done SPSPMM 2'
-
-  call psb_sphalo(am3,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         & a_err='Extend am3')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done sphalo/ rwxtd'
-
-  call psb_spspmm(op_restr,am3,ac,info)
-  if (info == psb_success_) call am3%free()
-  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
-
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         &a_err='Build ac = op_restr x am3')
-    goto 9999
-  end if
+!!$  ! Get the diagonal D
+!!$  adiag =  a%get_diag(info)
+!!$  if (info == psb_success_) &
+!!$       & call psb_realloc(ncol,adiag,info)    
+!!$  if (info == psb_success_) &
+!!$       & call psb_halo(adiag,desc_a,info)
+!!$  if (info == psb_success_) call a%cp_to_l(la)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  do i=1,size(adiag)
+!!$    if (adiag(i) /= dzero) then
+!!$      adinv(i) = done / adiag(i)
+!!$    else
+!!$      adinv(i) = done
+!!$    end if
+!!$  end do
+!!$
+!!$
+!!$
+!!$  ! 1. Allocate Ptilde in sparse matrix form 
+!!$  call op_prol%mv_to(tmpcoo)
+!!$  call ptilde%mv_from(tmpcoo)
+!!$  call ptilde%cscnv(info,type='csr')
+!!$
+!!$  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & ' Initial copies done.'
+!!$
+!!$  call da%scal(adinv,info)
+!!$
+!!$  call psb_spspmm(da,ptilde,dap,info)
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  call dap%clone(atmp,info)
+!!$
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  call psb_spspmm(da,atmp,dadap,info)
+!!$  call atmp%free()
+!!$
+!!$  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
+!!$  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
+!!$  call dap%mv_to(csc_dap)
+!!$  call dadap%mv_to(csc_dadap)
+!!$
+!!$  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
+!!$  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$  ! !$  write(0,*) trim(name),' OMP :',omp
+!!$  ! !$  write(0,*) trim(name),' ODEN:',oden
+!!$
+!!$  omp = omp/oden
+!!$
+!!$  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done NUMBMM 1'
+!!$
+!!$  call am3%mv_to(acsr3)
+!!$  ! Compute omega_int
+!!$  ommx = dzero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = dzero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
+!!$    end do
+!!$ ! !     if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero
+!!$    if(psb_minreal(omf(i)) < dzero) omf(i) = dzero
+!!$  end do
+!!$
+!!$  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
+!!$
+!!$  if (filter_mat) then
+!!$    !
+!!$    ! Build the filtered matrix Af from A
+!!$    ! 
+!!$    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
+!!$
+!!$    do i=1,nrow
+!!$      tmp = dzero
+!!$      jd  = -1 
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) jd = j 
+!!$        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
+!!$          tmp=tmp+acsrf%val(j)
+!!$          acsrf%val(j)=dzero
+!!$        endif
+!!$      enddo
+!!$      if (jd == -1) then 
+!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$      else
+!!$        acsrf%val(jd)=acsrf%val(jd)-tmp
+!!$      end if
+!!$    enddo
+!!$    ! Take out zeroed terms 
+!!$    call acsrf%clean_zeros(info)
+!!$
+!!$    !
+!!$    ! Build the smoothed prolongator using the filtered matrix
+!!$    ! 
+!!$    do i=1,acsrf%get_nrows()
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) then 
+!!$          acsrf%val(j) = done - omf(i)*acsrf%val(j) 
+!!$        else
+!!$          acsrf%val(j) = - omf(i)*acsrf%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$
+!!$    call af%mv_from(acsrf)
+!!$    !
+!!$    ! op_prol = (I-w*D*Af)Ptilde
+!!$    ! Doing it this way means to consider diag(Af_i)
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(af,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done SPSPMM 1'
+!!$  else
+!!$    !
+!!$    ! Build the smoothed prolongator using the original matrix
+!!$    !
+!!$    do i=1,acsr3%get_nrows()
+!!$      do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$        if (acsr3%ja(j) == i) then 
+!!$          acsr3%val(j) = done - omf(i)*acsr3%val(j) 
+!!$        else
+!!$          acsr3%val(j) = - omf(i)*acsr3%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    call am3%mv_from(acsr3)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$    !
+!!$    ! 
+!!$    ! op_prol = (I-w*D*A)Ptilde
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(am3,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done NUMBMM 1'
+!!$
+!!$  end if
+!!$
+!!$
+!!$  !
+!!$  ! Ok, let's start over with the restrictor
+!!$  ! 
+!!$  call ptilde%transc(rtilde)
+!!$  call la%cscnv(atmp,info,type='csr')
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.true.,rowscale=.true.)
+!!$  nrt  = am4%get_nrows() 
+!!$  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
+!!$  call atmp2%cscnv(info,type='CSR')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
+!!$  call am4%free()
+!!$  call atmp2%free()
+!!$
+!!$  ! This is to compute the transpose. It ONLY works if the
+!!$  ! original A has a symmetric pattern.
+!!$  call atmp%transc(atmp2) 
+!!$  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
+!!$  call dat%cscnv(info,type='csr')
+!!$  call dat%scal(adinv,info)
+!!$
+!!$  ! Now for the product. 
+!!$  call psb_spspmm(dat,ptilde,datp,info)
+!!$
+!!$  call datp%clone(atmp2,info)
+!!$  call psb_sphalo(atmp2,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$
+!!$  call psb_symbmm(dat,atmp2,datdatp,info)
+!!$  call psb_numbmm(dat,atmp2,datdatp)
+!!$  call atmp2%free()
+!!$
+!!$  call datp%mv_to(csc_datp)    
+!!$  call datdatp%mv_to(csc_datdatp)    
+!!$
+!!$  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
+!!$  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$
+!!$
+!!$  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
+!!$  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
+!!$  omp = omp/oden
+!!$  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
+!!$  ! Compute omega_int
+!!$  ommx = dzero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = dzero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  ! Going over the columns of atmp means going over the rows
+!!$  ! of A^T. Hopefully ;-) 
+!!$  call atmp%cp_to(acsc)
+!!$
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j= acsc%icp(i),acsc%icp(i+1)-1
+!!$      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
+!!$    end do
+!!$ ! !    if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = dzero
+!!$    if(psb_minreal(omf(i)) < dzero) omf(i) = dzero
+!!$  end do
+!!$  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
+!!$  call psb_halo(omf,desc_a,info)
+!!$  call acsc%free() 
+!!$
+!!$
+!!$  call atmp%mv_to(acsr1)
+!!$
+!!$  do i=1,acsr1%get_nrows()
+!!$    do j=acsr1%irp(i),acsr1%irp(i+1)-1
+!!$      if (acsr1%ja(j) == i) then 
+!!$        acsr1%val(j) = done - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      else
+!!$        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      end if
+!!$    end do
+!!$  end do
+!!$  call atmp%mv_from(acsr1)
+!!$
+!!$  call rtilde%mv_to(tmpcoo)
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call rtilde%mv_from(tmpcoo)
+!!$  call rtilde%cscnv(info,type='csr')
+!!$
+!!$  call psb_spspmm(rtilde,atmp,op_restr,info)
+!!$
+!!$  !
+!!$  ! Now we have to gather the halo of op_prol, and add it to itself
+!!$  ! to multiply it by A,
+!!$  !
+!!$  call op_prol%clone(tmp_prol,info) 
+!!$  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  !
+!!$  ! Now we have to fix this.  The only rows of B that are correct 
+!!$  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
+!!$  !
+!!$  call op_restr%mv_to(tmpcoo)
+!!$
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call op_restr%mv_from(tmpcoo)
+!!$  call op_restr%cscnv(info,type='csr')
+!!$
+!!$
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'starting sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(la,tmp_prol,am3,info)
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done SPSPMM 2'
+!!$
+!!$  call psb_sphalo(am3,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         & a_err='Extend am3')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(op_restr,am3,ac,info)
+!!$  if (info == psb_success_) call am3%free()
+!!$  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
+!!$
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         &a_err='Build ac = op_restr x am3')
+!!$    goto 9999
+!!$  end if
 
 
 

amgprec/impl/aggregator/amg_daggrmat_smth_bld.f90

@@ -116,7 +116,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)      :: parms 
-  type(psb_dspmat_type), intent(out)      :: op_prol,ac,op_restr
+  type(psb_dspmat_type), intent(inout)      :: op_prol,ac,op_restr
   type(psb_ldspmat_type), intent(inout)    :: t_prol
   type(psb_desc_type), intent(inout)       :: desc_ac
   integer(psb_ipk_), intent(out)           :: info

amgprec/impl/aggregator/amg_s_dec_aggregator_tprol.f90

@@ -83,8 +83,8 @@ subroutine  amg_s_dec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_s_dec_aggregator_type), target, intent(inout) :: ag
   type(amg_sml_parms), intent(inout)  :: parms
   type(amg_saggr_data), intent(in)    :: ag_data
-  type(psb_sspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_sspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)    :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lsspmat_type), intent(out)  :: t_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_s_parmatch_aggregator_tprol.F90

@@ -58,7 +58,7 @@ subroutine  amg_s_parmatch_aggregator_build_tprol(ag,parms,ag_data,&
   type(amg_sml_parms), intent(inout)   :: parms
   type(amg_saggr_data), intent(in)     :: ag_data
   type(psb_sspmat_type), intent(inout) :: a
-  type(psb_desc_type), intent(inout)   :: desc_a
+  type(psb_desc_type), intent(inout)      :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lsspmat_type), intent(out)  :: t_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_s_parmatch_smth_bld.F90

@@ -122,7 +122,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)    :: parms
   type(psb_lsspmat_type), intent(inout)  :: t_prol
-  type(psb_sspmat_type), intent(out)    :: op_prol,ac,op_restr
+  type(psb_sspmat_type), intent(inout)   :: op_prol,ac,op_restr
   type(psb_desc_type), intent(inout)    :: desc_ac
   integer(psb_ipk_), intent(out)           :: info
 

amgprec/impl/aggregator/amg_s_parmatch_spmm_bld.F90

@@ -108,7 +108,7 @@ subroutine amg_s_parmatch_spmm_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   ! Arguments
   type(psb_sspmat_type), intent(in)       :: a
-  type(psb_desc_type), intent(in)         :: desc_a
+  type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)      :: parms
   type(psb_lsspmat_type), intent(inout)   :: t_prol

amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_inner.F90

@@ -108,11 +108,11 @@ subroutine amg_s_parmatch_spmm_bld_inner(a_csr,desc_a,ilaggr,nlaggr,parms,&
 
   ! Arguments
   type(psb_s_csr_sparse_mat), intent(inout) :: a_csr
-  type(psb_desc_type), intent(in)           :: desc_a
+  type(psb_desc_type), intent(inout)          :: desc_a
   integer(psb_lpk_), intent(inout)          :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)        :: parms
   type(psb_lsspmat_type), intent(inout)     :: t_prol
-  type(psb_sspmat_type), intent(out)        :: ac, op_prol, op_restr
+  type(psb_sspmat_type), intent(inout)      :: ac, op_prol, op_restr
   type(psb_desc_type), intent(out)          :: desc_ac
   integer(psb_ipk_), intent(out)            :: info
 

amgprec/impl/aggregator/amg_s_parmatch_spmm_bld_ov.F90

@@ -108,7 +108,7 @@ subroutine amg_s_parmatch_spmm_bld_ov(a,desc_a,ilaggr,nlaggr,parms,&
 
   ! Arguments
   type(psb_sspmat_type), intent(inout)    :: a
-  type(psb_desc_type), intent(in)         :: desc_a
+  type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)        :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)      :: parms
   type(psb_lsspmat_type), intent(inout)   :: t_prol

amgprec/impl/aggregator/amg_s_symdec_aggregator_tprol.f90

@@ -86,8 +86,8 @@ subroutine  amg_s_symdec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_s_symdec_aggregator_type), target, intent(inout) :: ag
   type(amg_sml_parms), intent(inout)  :: parms
   type(amg_saggr_data), intent(in)    :: ag_data
-  type(psb_sspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_sspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)   :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lsspmat_type), intent(out)  :: op_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_saggrmat_minnrg_bld.f90

@@ -105,7 +105,7 @@
 !
 !
 subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,info)
+     & ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_s_inner_mod, amg_protect_name => amg_saggrmat_minnrg_bld
@@ -117,8 +117,8 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)         :: parms 
-  type(psb_lsspmat_type), intent(inout)       :: op_prol
-  type(psb_lsspmat_type), intent(out)         :: ac,op_restr
+  type(psb_lsspmat_type), intent(inout)       :: t_prol
+  type(psb_sspmat_type), intent(inout)       :: op_prol, ac,op_restr
   type(psb_desc_type), intent(inout)            :: desc_ac
   integer(psb_ipk_), intent(out)                :: info
 
@@ -171,6 +171,8 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  !NEEDS TO BE REWORKED !!
+  
   ! naggr: number of local aggregates
   ! nrow: local rows. 
   ! 
@@ -183,361 +185,361 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
     goto 9999      
   end if
 
-  ! Get the diagonal D
-  adiag =  a%get_diag(info)
-  if (info == psb_success_) &
-       & call psb_realloc(ncol,adiag,info)    
-  if (info == psb_success_) &
-       & call psb_halo(adiag,desc_a,info)
-  if (info == psb_success_) call a%cp_to_l(la)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
-    goto 9999
-  end if
-
-  do i=1,size(adiag)
-    if (adiag(i) /= szero) then
-      adinv(i) = sone / adiag(i)
-    else
-      adinv(i) = sone
-    end if
-  end do
-
-
-
-  ! 1. Allocate Ptilde in sparse matrix form 
-  call op_prol%mv_to(tmpcoo)
-  call ptilde%mv_from(tmpcoo)
-  call ptilde%cscnv(info,type='csr')
-
-  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
-  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & ' Initial copies done.'
-
-  call da%scal(adinv,info)
-
-  call psb_spspmm(da,ptilde,dap,info)
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
-    goto 9999
-  end if
-
-  call dap%clone(atmp,info)
-
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  call psb_spspmm(da,atmp,dadap,info)
-  call atmp%free()
-
-  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
-  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
-  call dap%mv_to(csc_dap)
-  call dadap%mv_to(csc_dadap)
-
-  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
-  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-  ! !$  write(0,*) trim(name),' OMP :',omp
-  ! !$  write(0,*) trim(name),' ODEN:',oden
-
-  omp = omp/oden
-
-  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done NUMBMM 1'
-
-  call am3%mv_to(acsr3)
-  ! Compute omega_int
-  ommx = szero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = szero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  do i=1, nrow
-    omf(i) = ommx
-    do j=acsr3%irp(i),acsr3%irp(i+1)-1
-      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero
-    if(psb_minreal(omf(i)) < szero) omf(i) = szero
-  end do
-
-  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
-
-  if (filter_mat) then
-    !
-    ! Build the filtered matrix Af from A
-    ! 
-    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
-
-    do i=1,nrow
-      tmp = szero
-      jd  = -1 
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) jd = j 
-        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
-          tmp=tmp+acsrf%val(j)
-          acsrf%val(j)=szero
-        endif
-      enddo
-      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
-        acsrf%val(jd)=acsrf%val(jd)-tmp
-      end if
-    enddo
-    ! Take out zeroed terms 
-    call acsrf%clean_zeros(info)
-
-    !
-    ! Build the smoothed prolongator using the filtered matrix
-    ! 
-    do i=1,acsrf%get_nrows()
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) then 
-          acsrf%val(j) = sone - omf(i)*acsrf%val(j) 
-        else
-          acsrf%val(j) = - omf(i)*acsrf%val(j) 
-        end if
-      end do
-    end do
-
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-
-    call af%mv_from(acsrf)
-    !
-    ! op_prol = (I-w*D*Af)Ptilde
-    ! Doing it this way means to consider diag(Af_i)
-    ! 
-    !
-    call psb_spspmm(af,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done SPSPMM 1'
-  else
-    !
-    ! Build the smoothed prolongator using the original matrix
-    !
-    do i=1,acsr3%get_nrows()
-      do j=acsr3%irp(i),acsr3%irp(i+1)-1
-        if (acsr3%ja(j) == i) then 
-          acsr3%val(j) = sone - omf(i)*acsr3%val(j) 
-        else
-          acsr3%val(j) = - omf(i)*acsr3%val(j) 
-        end if
-      end do
-    end do
-
-    call am3%mv_from(acsr3)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-    !
-    ! 
-    ! op_prol = (I-w*D*A)Ptilde
-    ! 
-    !
-    call psb_spspmm(am3,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done NUMBMM 1'
-
-  end if
-
-
-  !
-  ! Ok, let's start over with the restrictor
-  ! 
-  call ptilde%transc(rtilde)
-  call la%cscnv(atmp,info,type='csr')
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.true.,rowscale=.true.)
-  nrt  = am4%get_nrows() 
-  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
-  call atmp2%cscnv(info,type='CSR')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
-  call am4%free()
-  call atmp2%free()
-
-  ! This is to compute the transpose. It ONLY works if the
-  ! original A has a symmetric pattern.
-  call atmp%transc(atmp2) 
-  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
-  call dat%cscnv(info,type='csr')
-  call dat%scal(adinv,info)
-
-  ! Now for the product. 
-  call psb_spspmm(dat,ptilde,datp,info)
-
-  call datp%clone(atmp2,info)
-  call psb_sphalo(atmp2,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-
-  call psb_symbmm(dat,atmp2,datdatp,info)
-  call psb_numbmm(dat,atmp2,datdatp)
-  call atmp2%free()
-
-  call datp%mv_to(csc_datp)    
-  call datdatp%mv_to(csc_datdatp)    
-
-  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
-  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-
-
-  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
-  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
-  omp = omp/oden
-  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
-  ! Compute omega_int
-  ommx = szero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = szero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  ! Going over the columns of atmp means going over the rows
-  ! of A^T. Hopefully ;-) 
-  call atmp%cp_to(acsc)
-
-  do i=1, nrow
-    omf(i) = ommx
-    do j= acsc%icp(i),acsc%icp(i+1)-1
-      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero
-    if(psb_minreal(omf(i)) < szero) omf(i) = szero
-  end do
-  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
-  call psb_halo(omf,desc_a,info)
-  call acsc%free() 
-
-
-  call atmp%mv_to(acsr1)
-
-  do i=1,acsr1%get_nrows()
-    do j=acsr1%irp(i),acsr1%irp(i+1)-1
-      if (acsr1%ja(j) == i) then 
-        acsr1%val(j) = sone - acsr1%val(j)*omf(acsr1%ja(j))
-      else
-        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
-      end if
-    end do
-  end do
-  call atmp%mv_from(acsr1)
-
-  call rtilde%mv_to(tmpcoo)
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call rtilde%mv_from(tmpcoo)
-  call rtilde%cscnv(info,type='csr')
-
-  call psb_spspmm(rtilde,atmp,op_restr,info)
-
-  !
-  ! Now we have to gather the halo of op_prol, and add it to itself
-  ! to multiply it by A,
-  !
-  call op_prol%clone(tmp_prol,info) 
-  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
-    goto 9999
-  end if
-
-  !
-  ! Now we have to fix this.  The only rows of B that are correct 
-  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
-  !
-  call op_restr%mv_to(tmpcoo)
-
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call op_restr%mv_from(tmpcoo)
-  call op_restr%cscnv(info,type='csr')
-
-
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'starting sphalo/ rwxtd'
-
-  call psb_spspmm(la,tmp_prol,am3,info)
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done SPSPMM 2'
-
-  call psb_sphalo(am3,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         & a_err='Extend am3')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done sphalo/ rwxtd'
-
-  call psb_spspmm(op_restr,am3,ac,info)
-  if (info == psb_success_) call am3%free()
-  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
-
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         &a_err='Build ac = op_restr x am3')
-    goto 9999
-  end if
+!!$  ! Get the diagonal D
+!!$  adiag =  a%get_diag(info)
+!!$  if (info == psb_success_) &
+!!$       & call psb_realloc(ncol,adiag,info)    
+!!$  if (info == psb_success_) &
+!!$       & call psb_halo(adiag,desc_a,info)
+!!$  if (info == psb_success_) call a%cp_to_l(la)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  do i=1,size(adiag)
+!!$    if (adiag(i) /= szero) then
+!!$      adinv(i) = sone / adiag(i)
+!!$    else
+!!$      adinv(i) = sone
+!!$    end if
+!!$  end do
+!!$
+!!$
+!!$
+!!$  ! 1. Allocate Ptilde in sparse matrix form 
+!!$  call op_prol%mv_to(tmpcoo)
+!!$  call ptilde%mv_from(tmpcoo)
+!!$  call ptilde%cscnv(info,type='csr')
+!!$
+!!$  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & ' Initial copies done.'
+!!$
+!!$  call da%scal(adinv,info)
+!!$
+!!$  call psb_spspmm(da,ptilde,dap,info)
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  call dap%clone(atmp,info)
+!!$
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  call psb_spspmm(da,atmp,dadap,info)
+!!$  call atmp%free()
+!!$
+!!$  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
+!!$  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
+!!$  call dap%mv_to(csc_dap)
+!!$  call dadap%mv_to(csc_dadap)
+!!$
+!!$  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
+!!$  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$  ! !$  write(0,*) trim(name),' OMP :',omp
+!!$  ! !$  write(0,*) trim(name),' ODEN:',oden
+!!$
+!!$  omp = omp/oden
+!!$
+!!$  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done NUMBMM 1'
+!!$
+!!$  call am3%mv_to(acsr3)
+!!$  ! Compute omega_int
+!!$  ommx = szero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = szero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
+!!$    end do
+!!$ ! !     if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero
+!!$    if(psb_minreal(omf(i)) < szero) omf(i) = szero
+!!$  end do
+!!$
+!!$  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
+!!$
+!!$  if (filter_mat) then
+!!$    !
+!!$    ! Build the filtered matrix Af from A
+!!$    ! 
+!!$    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
+!!$
+!!$    do i=1,nrow
+!!$      tmp = szero
+!!$      jd  = -1 
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) jd = j 
+!!$        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
+!!$          tmp=tmp+acsrf%val(j)
+!!$          acsrf%val(j)=szero
+!!$        endif
+!!$      enddo
+!!$      if (jd == -1) then 
+!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$      else
+!!$        acsrf%val(jd)=acsrf%val(jd)-tmp
+!!$      end if
+!!$    enddo
+!!$    ! Take out zeroed terms 
+!!$    call acsrf%clean_zeros(info)
+!!$
+!!$    !
+!!$    ! Build the smoothed prolongator using the filtered matrix
+!!$    ! 
+!!$    do i=1,acsrf%get_nrows()
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) then 
+!!$          acsrf%val(j) = sone - omf(i)*acsrf%val(j) 
+!!$        else
+!!$          acsrf%val(j) = - omf(i)*acsrf%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$
+!!$    call af%mv_from(acsrf)
+!!$    !
+!!$    ! op_prol = (I-w*D*Af)Ptilde
+!!$    ! Doing it this way means to consider diag(Af_i)
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(af,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done SPSPMM 1'
+!!$  else
+!!$    !
+!!$    ! Build the smoothed prolongator using the original matrix
+!!$    !
+!!$    do i=1,acsr3%get_nrows()
+!!$      do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$        if (acsr3%ja(j) == i) then 
+!!$          acsr3%val(j) = sone - omf(i)*acsr3%val(j) 
+!!$        else
+!!$          acsr3%val(j) = - omf(i)*acsr3%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    call am3%mv_from(acsr3)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$    !
+!!$    ! 
+!!$    ! op_prol = (I-w*D*A)Ptilde
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(am3,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done NUMBMM 1'
+!!$
+!!$  end if
+!!$
+!!$
+!!$  !
+!!$  ! Ok, let's start over with the restrictor
+!!$  ! 
+!!$  call ptilde%transc(rtilde)
+!!$  call la%cscnv(atmp,info,type='csr')
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.true.,rowscale=.true.)
+!!$  nrt  = am4%get_nrows() 
+!!$  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
+!!$  call atmp2%cscnv(info,type='CSR')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
+!!$  call am4%free()
+!!$  call atmp2%free()
+!!$
+!!$  ! This is to compute the transpose. It ONLY works if the
+!!$  ! original A has a symmetric pattern.
+!!$  call atmp%transc(atmp2) 
+!!$  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
+!!$  call dat%cscnv(info,type='csr')
+!!$  call dat%scal(adinv,info)
+!!$
+!!$  ! Now for the product. 
+!!$  call psb_spspmm(dat,ptilde,datp,info)
+!!$
+!!$  call datp%clone(atmp2,info)
+!!$  call psb_sphalo(atmp2,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$
+!!$  call psb_symbmm(dat,atmp2,datdatp,info)
+!!$  call psb_numbmm(dat,atmp2,datdatp)
+!!$  call atmp2%free()
+!!$
+!!$  call datp%mv_to(csc_datp)    
+!!$  call datdatp%mv_to(csc_datdatp)    
+!!$
+!!$  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
+!!$  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$
+!!$
+!!$  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
+!!$  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
+!!$  omp = omp/oden
+!!$  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
+!!$  ! Compute omega_int
+!!$  ommx = szero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = szero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  ! Going over the columns of atmp means going over the rows
+!!$  ! of A^T. Hopefully ;-) 
+!!$  call atmp%cp_to(acsc)
+!!$
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j= acsc%icp(i),acsc%icp(i+1)-1
+!!$      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
+!!$    end do
+!!$ ! !    if(min(real(omf(i)),aimag(omf(i))) < szero) omf(i) = szero
+!!$    if(psb_minreal(omf(i)) < szero) omf(i) = szero
+!!$  end do
+!!$  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
+!!$  call psb_halo(omf,desc_a,info)
+!!$  call acsc%free() 
+!!$
+!!$
+!!$  call atmp%mv_to(acsr1)
+!!$
+!!$  do i=1,acsr1%get_nrows()
+!!$    do j=acsr1%irp(i),acsr1%irp(i+1)-1
+!!$      if (acsr1%ja(j) == i) then 
+!!$        acsr1%val(j) = sone - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      else
+!!$        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      end if
+!!$    end do
+!!$  end do
+!!$  call atmp%mv_from(acsr1)
+!!$
+!!$  call rtilde%mv_to(tmpcoo)
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call rtilde%mv_from(tmpcoo)
+!!$  call rtilde%cscnv(info,type='csr')
+!!$
+!!$  call psb_spspmm(rtilde,atmp,op_restr,info)
+!!$
+!!$  !
+!!$  ! Now we have to gather the halo of op_prol, and add it to itself
+!!$  ! to multiply it by A,
+!!$  !
+!!$  call op_prol%clone(tmp_prol,info) 
+!!$  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  !
+!!$  ! Now we have to fix this.  The only rows of B that are correct 
+!!$  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
+!!$  !
+!!$  call op_restr%mv_to(tmpcoo)
+!!$
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call op_restr%mv_from(tmpcoo)
+!!$  call op_restr%cscnv(info,type='csr')
+!!$
+!!$
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'starting sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(la,tmp_prol,am3,info)
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done SPSPMM 2'
+!!$
+!!$  call psb_sphalo(am3,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         & a_err='Extend am3')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(op_restr,am3,ac,info)
+!!$  if (info == psb_success_) call am3%free()
+!!$  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
+!!$
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         &a_err='Build ac = op_restr x am3')
+!!$    goto 9999
+!!$  end if
 
 
 

amgprec/impl/aggregator/amg_saggrmat_smth_bld.f90

@@ -116,7 +116,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
   type(amg_sml_parms), intent(inout)      :: parms 
-  type(psb_sspmat_type), intent(out)      :: op_prol,ac,op_restr
+  type(psb_sspmat_type), intent(inout)      :: op_prol,ac,op_restr
   type(psb_lsspmat_type), intent(inout)    :: t_prol
   type(psb_desc_type), intent(inout)       :: desc_ac
   integer(psb_ipk_), intent(out)           :: info

amgprec/impl/aggregator/amg_z_dec_aggregator_tprol.f90

@@ -83,8 +83,8 @@ subroutine  amg_z_dec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_z_dec_aggregator_type), target, intent(inout) :: ag
   type(amg_dml_parms), intent(inout)  :: parms
   type(amg_daggr_data), intent(in)    :: ag_data
-  type(psb_zspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_zspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)    :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lzspmat_type), intent(out)  :: t_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_z_symdec_aggregator_tprol.f90

@@ -86,8 +86,8 @@ subroutine  amg_z_symdec_aggregator_build_tprol(ag,parms,ag_data,&
   class(amg_z_symdec_aggregator_type), target, intent(inout) :: ag
   type(amg_dml_parms), intent(inout)  :: parms
   type(amg_daggr_data), intent(in)    :: ag_data
-  type(psb_zspmat_type), intent(in)   :: a
-  type(psb_desc_type), intent(in)     :: desc_a
+  type(psb_zspmat_type), intent(inout) :: a
+  type(psb_desc_type), intent(inout)   :: desc_a
   integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
   type(psb_lzspmat_type), intent(out)  :: op_prol
   integer(psb_ipk_), intent(out)      :: info

amgprec/impl/aggregator/amg_zaggrmat_minnrg_bld.f90

@@ -105,7 +105,7 @@
 !
 !
 subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
-     & ac,desc_ac,op_prol,op_restr,info)
+     & ac,desc_ac,op_prol,op_restr,t_prol,info)
   use psb_base_mod
   use amg_base_prec_type
   use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_minnrg_bld
@@ -117,8 +117,8 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)            :: desc_a
   integer(psb_lpk_), intent(inout)              :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)         :: parms 
-  type(psb_lzspmat_type), intent(inout)       :: op_prol
-  type(psb_lzspmat_type), intent(out)         :: ac,op_restr
+  type(psb_lzspmat_type), intent(inout)       :: t_prol
+  type(psb_zspmat_type), intent(inout)       :: op_prol, ac,op_restr
   type(psb_desc_type), intent(inout)            :: desc_ac
   integer(psb_ipk_), intent(out)                :: info
 
@@ -171,6 +171,8 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 
   filter_mat = (parms%aggr_filter == amg_filter_mat_)
 
+  !NEEDS TO BE REWORKED !!
+  
   ! naggr: number of local aggregates
   ! nrow: local rows. 
   ! 
@@ -183,361 +185,361 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
     goto 9999      
   end if
 
-  ! Get the diagonal D
-  adiag =  a%get_diag(info)
-  if (info == psb_success_) &
-       & call psb_realloc(ncol,adiag,info)    
-  if (info == psb_success_) &
-       & call psb_halo(adiag,desc_a,info)
-  if (info == psb_success_) call a%cp_to_l(la)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
-    goto 9999
-  end if
-
-  do i=1,size(adiag)
-    if (adiag(i) /= zzero) then
-      adinv(i) = zone / adiag(i)
-    else
-      adinv(i) = zone
-    end if
-  end do
-
-
-
-  ! 1. Allocate Ptilde in sparse matrix form 
-  call op_prol%mv_to(tmpcoo)
-  call ptilde%mv_from(tmpcoo)
-  call ptilde%cscnv(info,type='csr')
-
-  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
-  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & ' Initial copies done.'
-
-  call da%scal(adinv,info)
-
-  call psb_spspmm(da,ptilde,dap,info)
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
-    goto 9999
-  end if
-
-  call dap%clone(atmp,info)
-
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  call psb_spspmm(da,atmp,dadap,info)
-  call atmp%free()
-
-  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
-  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
-  call dap%mv_to(csc_dap)
-  call dadap%mv_to(csc_dadap)
-
-  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
-  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-  ! !$  write(0,*) trim(name),' OMP :',omp
-  ! !$  write(0,*) trim(name),' ODEN:',oden
-
-  omp = omp/oden
-
-  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done NUMBMM 1'
-
-  call am3%mv_to(acsr3)
-  ! Compute omega_int
-  ommx = zzero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = zzero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  do i=1, nrow
-    omf(i) = ommx
-    do j=acsr3%irp(i),acsr3%irp(i+1)-1
-      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero
-    if(psb_minreal(omf(i)) < dzero) omf(i) = zzero
-  end do
-
-  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
-
-  if (filter_mat) then
-    !
-    ! Build the filtered matrix Af from A
-    ! 
-    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
-
-    do i=1,nrow
-      tmp = zzero
-      jd  = -1 
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) jd = j 
-        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
-          tmp=tmp+acsrf%val(j)
-          acsrf%val(j)=zzero
-        endif
-      enddo
-      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
-      else
-        acsrf%val(jd)=acsrf%val(jd)-tmp
-      end if
-    enddo
-    ! Take out zeroed terms 
-    call acsrf%clean_zeros(info)
-
-    !
-    ! Build the smoothed prolongator using the filtered matrix
-    ! 
-    do i=1,acsrf%get_nrows()
-      do j=acsrf%irp(i),acsrf%irp(i+1)-1
-        if (acsrf%ja(j) == i) then 
-          acsrf%val(j) = zone - omf(i)*acsrf%val(j) 
-        else
-          acsrf%val(j) = - omf(i)*acsrf%val(j) 
-        end if
-      end do
-    end do
-
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-
-    call af%mv_from(acsrf)
-    !
-    ! op_prol = (I-w*D*Af)Ptilde
-    ! Doing it this way means to consider diag(Af_i)
-    ! 
-    !
-    call psb_spspmm(af,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done SPSPMM 1'
-  else
-    !
-    ! Build the smoothed prolongator using the original matrix
-    !
-    do i=1,acsr3%get_nrows()
-      do j=acsr3%irp(i),acsr3%irp(i+1)-1
-        if (acsr3%ja(j) == i) then 
-          acsr3%val(j) = zone - omf(i)*acsr3%val(j) 
-        else
-          acsr3%val(j) = - omf(i)*acsr3%val(j) 
-        end if
-      end do
-    end do
-
-    call am3%mv_from(acsr3)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done gather, going for SYMBMM 1'
-    !
-    ! 
-    ! op_prol = (I-w*D*A)Ptilde
-    ! 
-    !
-    call psb_spspmm(am3,ptilde,op_prol,info)
-    if (debug_level >= psb_debug_outer_) &
-         & write(debug_unit,*) me,' ',trim(name),&
-         & 'Done NUMBMM 1'
-
-  end if
-
-
-  !
-  ! Ok, let's start over with the restrictor
-  ! 
-  call ptilde%transc(rtilde)
-  call la%cscnv(atmp,info,type='csr')
-  call psb_sphalo(atmp,desc_a,am4,info,&
-       & colcnv=.true.,rowscale=.true.)
-  nrt  = am4%get_nrows() 
-  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
-  call atmp2%cscnv(info,type='CSR')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
-  call am4%free()
-  call atmp2%free()
-
-  ! This is to compute the transpose. It ONLY works if the
-  ! original A has a symmetric pattern.
-  call atmp%transc(atmp2) 
-  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
-  call dat%cscnv(info,type='csr')
-  call dat%scal(adinv,info)
-
-  ! Now for the product. 
-  call psb_spspmm(dat,ptilde,datp,info)
-
-  call datp%clone(atmp2,info)
-  call psb_sphalo(atmp2,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
-  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-
-  call psb_symbmm(dat,atmp2,datdatp,info)
-  call psb_numbmm(dat,atmp2,datdatp)
-  call atmp2%free()
-
-  call datp%mv_to(csc_datp)    
-  call datdatp%mv_to(csc_datdatp)    
-
-  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
-  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
-  call psb_sum(ctxt,omp)
-  call psb_sum(ctxt,oden)
-
-
-  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
-  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
-  omp = omp/oden
-  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
-  ! Compute omega_int
-  ommx = zzero
-  do i=1, ncol
-    if (ilaggr(i) >0) then 
-      omi(i) = omp(ilaggr(i))
-    else
-      omi(i) = zzero
-    end if
-    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
-  end do
-  ! Compute omega_fine
-  ! Going over the columns of atmp means going over the rows
-  ! of A^T. Hopefully ;-) 
-  call atmp%cp_to(acsc)
-
-  do i=1, nrow
-    omf(i) = ommx
-    do j= acsc%icp(i),acsc%icp(i+1)-1
-      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
-    end do
-!!$    if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero
-    if(psb_minreal(omf(i)) < dzero) omf(i) = zzero
-  end do
-  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
-  call psb_halo(omf,desc_a,info)
-  call acsc%free() 
-
-
-  call atmp%mv_to(acsr1)
-
-  do i=1,acsr1%get_nrows()
-    do j=acsr1%irp(i),acsr1%irp(i+1)-1
-      if (acsr1%ja(j) == i) then 
-        acsr1%val(j) = zone - acsr1%val(j)*omf(acsr1%ja(j))
-      else
-        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
-      end if
-    end do
-  end do
-  call atmp%mv_from(acsr1)
-
-  call rtilde%mv_to(tmpcoo)
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call rtilde%mv_from(tmpcoo)
-  call rtilde%cscnv(info,type='csr')
-
-  call psb_spspmm(rtilde,atmp,op_restr,info)
-
-  !
-  ! Now we have to gather the halo of op_prol, and add it to itself
-  ! to multiply it by A,
-  !
-  call op_prol%clone(tmp_prol,info) 
-  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
-    goto 9999
-  end if
-
-  !
-  ! Now we have to fix this.  The only rows of B that are correct 
-  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
-  !
-  call op_restr%mv_to(tmpcoo)
-
-  nzl = tmpcoo%get_nzeros()
-  i=0
-  do k=1, nzl
-    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
-      i = i+1
-      tmpcoo%val(i) = tmpcoo%val(k)
-      tmpcoo%ia(i)  = tmpcoo%ia(k)
-      tmpcoo%ja(i)  = tmpcoo%ja(k)
-    end if
-  end do
-  call tmpcoo%set_nzeros(i)
-  call op_restr%mv_from(tmpcoo)
-  call op_restr%cscnv(info,type='csr')
-
-
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'starting sphalo/ rwxtd'
-
-  call psb_spspmm(la,tmp_prol,am3,info)
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done SPSPMM 2'
-
-  call psb_sphalo(am3,desc_a,am4,info,&
-       & colcnv=.false.,rowscale=.true.)
-  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
-  if (info == psb_success_) call am4%free()
-
-  if(info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         & a_err='Extend am3')
-    goto 9999
-  end if
-  if (debug_level >= psb_debug_outer_) &
-       & write(debug_unit,*) me,' ',trim(name),&
-       & 'Done sphalo/ rwxtd'
-
-  call psb_spspmm(op_restr,am3,ac,info)
-  if (info == psb_success_) call am3%free()
-  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
-
-  if (info /= psb_success_) then
-    call psb_errpush(psb_err_internal_error_,name,&
-         &a_err='Build ac = op_restr x am3')
-    goto 9999
-  end if
+!!$  ! Get the diagonal D
+!!$  adiag =  a%get_diag(info)
+!!$  if (info == psb_success_) &
+!!$       & call psb_realloc(ncol,adiag,info)    
+!!$  if (info == psb_success_) &
+!!$       & call psb_halo(adiag,desc_a,info)
+!!$  if (info == psb_success_) call a%cp_to_l(la)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  do i=1,size(adiag)
+!!$    if (adiag(i) /= zzero) then
+!!$      adinv(i) = zone / adiag(i)
+!!$    else
+!!$      adinv(i) = zone
+!!$    end if
+!!$  end do
+!!$
+!!$
+!!$
+!!$  ! 1. Allocate Ptilde in sparse matrix form 
+!!$  call op_prol%mv_to(tmpcoo)
+!!$  call ptilde%mv_from(tmpcoo)
+!!$  call ptilde%cscnv(info,type='csr')
+!!$
+!!$  if (info == psb_success_) call la%cscnv(am3,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info == psb_success_) call la%cscnv(da,info,type='csr',dupl=psb_dupl_add_)
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spcnv')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & ' Initial copies done.'
+!!$
+!!$  call da%scal(adinv,info)
+!!$
+!!$  call psb_spspmm(da,ptilde,dap,info)
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  call dap%clone(atmp,info)
+!!$
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  call psb_spspmm(da,atmp,dadap,info)
+!!$  call atmp%free()
+!!$
+!!$  !  !$  write(0,*) 'Columns of AP',psb_sp_get_ncols(ap)
+!!$  !  !$  write(0,*) 'Columns of ADAP',psb_sp_get_ncols(adap)
+!!$  call dap%mv_to(csc_dap)
+!!$  call dadap%mv_to(csc_dadap)
+!!$
+!!$  call csc_mat_col_prod(csc_dap,csc_dadap,omp,info)
+!!$  call csc_mat_col_prod(csc_dadap,csc_dadap,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$  ! !$  write(0,*) trim(name),' OMP :',omp
+!!$  ! !$  write(0,*) trim(name),' ODEN:',oden
+!!$
+!!$  omp = omp/oden
+!!$
+!!$  ! !$  write(0,*) 'Check on output prolongator ',omp(1:min(size(omp),10))
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done NUMBMM 1'
+!!$
+!!$  call am3%mv_to(acsr3)
+!!$  ! Compute omega_int
+!!$  ommx = zzero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = zzero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$      if(abs(omi(acsr3%ja(j))) .lt. abs(omf(i))) omf(i)=omi(acsr3%ja(j))
+!!$    end do
+!!$ ! !     if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero
+!!$    if(psb_minreal(omf(i)) < dzero) omf(i) = zzero
+!!$  end do
+!!$
+!!$  omf(1:nrow) = omf(1:nrow) * adinv(1:nrow)
+!!$
+!!$  if (filter_mat) then
+!!$    !
+!!$    ! Build the filtered matrix Af from A
+!!$    ! 
+!!$    call la%cscnv(acsrf,info,dupl=psb_dupl_add_)
+!!$
+!!$    do i=1,nrow
+!!$      tmp = zzero
+!!$      jd  = -1 
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) jd = j 
+!!$        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
+!!$          tmp=tmp+acsrf%val(j)
+!!$          acsrf%val(j)=zzero
+!!$        endif
+!!$      enddo
+!!$      if (jd == -1) then 
+!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$      else
+!!$        acsrf%val(jd)=acsrf%val(jd)-tmp
+!!$      end if
+!!$    enddo
+!!$    ! Take out zeroed terms 
+!!$    call acsrf%clean_zeros(info)
+!!$
+!!$    !
+!!$    ! Build the smoothed prolongator using the filtered matrix
+!!$    ! 
+!!$    do i=1,acsrf%get_nrows()
+!!$      do j=acsrf%irp(i),acsrf%irp(i+1)-1
+!!$        if (acsrf%ja(j) == i) then 
+!!$          acsrf%val(j) = zone - omf(i)*acsrf%val(j) 
+!!$        else
+!!$          acsrf%val(j) = - omf(i)*acsrf%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$
+!!$    call af%mv_from(acsrf)
+!!$    !
+!!$    ! op_prol = (I-w*D*Af)Ptilde
+!!$    ! Doing it this way means to consider diag(Af_i)
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(af,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done SPSPMM 1'
+!!$  else
+!!$    !
+!!$    ! Build the smoothed prolongator using the original matrix
+!!$    !
+!!$    do i=1,acsr3%get_nrows()
+!!$      do j=acsr3%irp(i),acsr3%irp(i+1)-1
+!!$        if (acsr3%ja(j) == i) then 
+!!$          acsr3%val(j) = zone - omf(i)*acsr3%val(j) 
+!!$        else
+!!$          acsr3%val(j) = - omf(i)*acsr3%val(j) 
+!!$        end if
+!!$      end do
+!!$    end do
+!!$
+!!$    call am3%mv_from(acsr3)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done gather, going for SYMBMM 1'
+!!$    !
+!!$    ! 
+!!$    ! op_prol = (I-w*D*A)Ptilde
+!!$    ! 
+!!$    !
+!!$    call psb_spspmm(am3,ptilde,op_prol,info)
+!!$    if (debug_level >= psb_debug_outer_) &
+!!$         & write(debug_unit,*) me,' ',trim(name),&
+!!$         & 'Done NUMBMM 1'
+!!$
+!!$  end if
+!!$
+!!$
+!!$  !
+!!$  ! Ok, let's start over with the restrictor
+!!$  ! 
+!!$  call ptilde%transc(rtilde)
+!!$  call la%cscnv(atmp,info,type='csr')
+!!$  call psb_sphalo(atmp,desc_a,am4,info,&
+!!$       & colcnv=.true.,rowscale=.true.)
+!!$  nrt  = am4%get_nrows() 
+!!$  call am4%csclip(atmp2,info,lone,nrt,lone,ncol)
+!!$  call atmp2%cscnv(info,type='CSR')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp,info,b=atmp2)      
+!!$  call am4%free()
+!!$  call atmp2%free()
+!!$
+!!$  ! This is to compute the transpose. It ONLY works if the
+!!$  ! original A has a symmetric pattern.
+!!$  call atmp%transc(atmp2) 
+!!$  call atmp2%csclip(dat,info,lone,nrow,lone,ncol)
+!!$  call dat%cscnv(info,type='csr')
+!!$  call dat%scal(adinv,info)
+!!$
+!!$  ! Now for the product. 
+!!$  call psb_spspmm(dat,ptilde,datp,info)
+!!$
+!!$  call datp%clone(atmp2,info)
+!!$  call psb_sphalo(atmp2,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.,outfmt='CSR  ')
+!!$  if (info == psb_success_) call psb_rwextd(ncol,atmp2,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$
+!!$  call psb_symbmm(dat,atmp2,datdatp,info)
+!!$  call psb_numbmm(dat,atmp2,datdatp)
+!!$  call atmp2%free()
+!!$
+!!$  call datp%mv_to(csc_datp)    
+!!$  call datdatp%mv_to(csc_datdatp)    
+!!$
+!!$  call csc_mat_col_prod(csc_datp,csc_datdatp,omp,info)
+!!$  call csc_mat_col_prod(csc_datdatp,csc_datdatp,oden,info)
+!!$  call psb_sum(ctxt,omp)
+!!$  call psb_sum(ctxt,oden)
+!!$
+!!$
+!!$  ! !$  write(debug_unit,*) trim(name),' OMP_R :',omp
+!!$  ! ! $  write(debug_unit,*) trim(name),' ODEN_R:',oden
+!!$  omp = omp/oden
+!!$  ! !$  write(0,*) 'Check on output restrictor',omp(1:min(size(omp),10))
+!!$  ! Compute omega_int
+!!$  ommx = zzero
+!!$  do i=1, ncol
+!!$    if (ilaggr(i) >0) then 
+!!$      omi(i) = omp(ilaggr(i))
+!!$    else
+!!$      omi(i) = zzero
+!!$    end if
+!!$    if(abs(omi(i)) .gt. abs(ommx)) ommx = omi(i)
+!!$  end do
+!!$  ! Compute omega_fine
+!!$  ! Going over the columns of atmp means going over the rows
+!!$  ! of A^T. Hopefully ;-) 
+!!$  call atmp%cp_to(acsc)
+!!$
+!!$  do i=1, nrow
+!!$    omf(i) = ommx
+!!$    do j= acsc%icp(i),acsc%icp(i+1)-1
+!!$      if(abs(omi(acsc%ia(j))) .lt. abs(omf(i))) omf(i)=omi(acsc%ia(j))
+!!$    end do
+!!$ ! !    if(min(real(omf(i)),aimag(omf(i))) < dzero) omf(i) = zzero
+!!$    if(psb_minreal(omf(i)) < dzero) omf(i) = zzero
+!!$  end do
+!!$  omf(1:nrow) = omf(1:nrow)*adinv(1:nrow)
+!!$  call psb_halo(omf,desc_a,info)
+!!$  call acsc%free() 
+!!$
+!!$
+!!$  call atmp%mv_to(acsr1)
+!!$
+!!$  do i=1,acsr1%get_nrows()
+!!$    do j=acsr1%irp(i),acsr1%irp(i+1)-1
+!!$      if (acsr1%ja(j) == i) then 
+!!$        acsr1%val(j) = zone - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      else
+!!$        acsr1%val(j) =      - acsr1%val(j)*omf(acsr1%ja(j))
+!!$      end if
+!!$    end do
+!!$  end do
+!!$  call atmp%mv_from(acsr1)
+!!$
+!!$  call rtilde%mv_to(tmpcoo)
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call rtilde%mv_from(tmpcoo)
+!!$  call rtilde%cscnv(info,type='csr')
+!!$
+!!$  call psb_spspmm(rtilde,atmp,op_restr,info)
+!!$
+!!$  !
+!!$  ! Now we have to gather the halo of op_prol, and add it to itself
+!!$  ! to multiply it by A,
+!!$  !
+!!$  call op_prol%clone(tmp_prol,info) 
+!!$  if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
+!!$    goto 9999
+!!$  end if
+!!$
+!!$  !
+!!$  ! Now we have to fix this.  The only rows of B that are correct 
+!!$  ! are those corresponding to "local" aggregates, i.e. indices in ilaggr(:)
+!!$  !
+!!$  call op_restr%mv_to(tmpcoo)
+!!$
+!!$  nzl = tmpcoo%get_nzeros()
+!!$  i=0
+!!$  do k=1, nzl
+!!$    if ((naggrm1 < tmpcoo%ia(k)) .and. (tmpcoo%ia(k) <= naggrp1)) then
+!!$      i = i+1
+!!$      tmpcoo%val(i) = tmpcoo%val(k)
+!!$      tmpcoo%ia(i)  = tmpcoo%ia(k)
+!!$      tmpcoo%ja(i)  = tmpcoo%ja(k)
+!!$    end if
+!!$  end do
+!!$  call tmpcoo%set_nzeros(i)
+!!$  call op_restr%mv_from(tmpcoo)
+!!$  call op_restr%cscnv(info,type='csr')
+!!$
+!!$
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'starting sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(la,tmp_prol,am3,info)
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done SPSPMM 2'
+!!$
+!!$  call psb_sphalo(am3,desc_a,am4,info,&
+!!$       & colcnv=.false.,rowscale=.true.)
+!!$  if (info == psb_success_) call psb_rwextd(ncol,am3,info,b=am4)      
+!!$  if (info == psb_success_) call am4%free()
+!!$
+!!$  if(info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         & a_err='Extend am3')
+!!$    goto 9999
+!!$  end if
+!!$  if (debug_level >= psb_debug_outer_) &
+!!$       & write(debug_unit,*) me,' ',trim(name),&
+!!$       & 'Done sphalo/ rwxtd'
+!!$
+!!$  call psb_spspmm(op_restr,am3,ac,info)
+!!$  if (info == psb_success_) call am3%free()
+!!$  if (info == psb_success_) call ac%cscnv(info,type='coo',dupl=psb_dupl_add_)
+!!$
+!!$  if (info /= psb_success_) then
+!!$    call psb_errpush(psb_err_internal_error_,name,&
+!!$         &a_err='Build ac = op_restr x am3')
+!!$    goto 9999
+!!$  end if
 
 
 

amgprec/impl/aggregator/amg_zaggrmat_smth_bld.f90

@@ -116,7 +116,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
   type(psb_desc_type), intent(inout)       :: desc_a
   integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
   type(amg_dml_parms), intent(inout)      :: parms 
-  type(psb_zspmat_type), intent(out)      :: op_prol,ac,op_restr
+  type(psb_zspmat_type), intent(inout)      :: op_prol,ac,op_restr
   type(psb_lzspmat_type), intent(inout)    :: t_prol
   type(psb_desc_type), intent(inout)       :: desc_ac
   integer(psb_ipk_), intent(out)           :: info

amgprec/impl/amg_ccprecset.F90

@@ -462,6 +462,14 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
           select case (psb_toupper(string))
           case('BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_SLU_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+            call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+            call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
           case('L1-BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_SLU_)
@@ -612,6 +620,14 @@ subroutine amg_ccprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
         select case (psb_toupper(trim(string)))
         case('BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_SLU_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+          call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+          call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
         case('L1-BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_SLU_)

amgprec/impl/amg_cfile_prec_descr.f90

@@ -65,7 +65,7 @@
 !         0: normal
 !        >1: increased details 
 !
-subroutine amg_cfile_prec_descr(prec,iout,root, verbosity)
+subroutine amg_cfile_prec_descr(prec,info,iout,root, verbosity)
   use psb_base_mod
   use amg_c_prec_mod, amg_protect_name => amg_cfile_prec_descr
   use amg_c_inner_mod
@@ -74,13 +74,14 @@ subroutine amg_cfile_prec_descr(prec,iout,root, verbosity)
   implicit none 
   ! Arguments
   class(amg_cprec_type), intent(in)      :: prec
+  integer(psb_ipk_), intent(out) :: info
   integer(psb_ipk_), intent(in), optional :: iout
   integer(psb_ipk_), intent(in), optional :: root
   integer(psb_ipk_), intent(in), optional   :: verbosity
 
 
   ! Local variables
-  integer(psb_ipk_)   :: ilev, nlev, ilmin, info, nswps
+  integer(psb_ipk_)   :: ilev, nlev, ilmin, nswps
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   logical             :: is_symgs

amgprec/impl/amg_cprecinit.F90

@@ -1,15 +1,15 @@
-!   
-!   
+!
+!
 !                             AMG4PSBLAS version 1.0
 !    Algebraic Multigrid Package
 !               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2021 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Fabio Durastante
+!
 !    Redistribution and use in source and binary forms, with or without
 !    modification, are permitted provided that the following conditions
 !    are met:
@@ -21,7 +21,7 @@
 !      3. The name of the AMG4PSBLAS group or the names of its contributors may
 !         not be used to endorse or promote products derived from this
 !         software without specific written permission.
-!   
+!
 !    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@@ -33,8 +33,8 @@
 !    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
+!
+!
 ! File: amg_cprecinit.f90
 !
 ! Subroutine: amg_cprecinit
@@ -42,21 +42,21 @@
 !
 !  This routine allocates and initializes the preconditioner data structure,
 !  according to the preconditioner type chosen by the user.
-!  
+!
 !  A default preconditioner is set for each preconditioner type
 !  specified by the user:
 !
 !    'NOPREC'         - no preconditioner
 !
-!    'DIAG', 'JACOBI' - diagonal/Jacobi           
+!    'DIAG', 'JACOBI' - diagonal/Jacobi
 !
 !    'L1-DIAG', 'L1-JACOBI' - diagonal/Jacobi  with L1 norm correction
 !
 !    'GS', 'FBGS'     - Hybrid Gauss-Seidel, also symmetrized
-!                       
+!
 !    'BJAC'           - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks
-!                       
+!
 !    'L1-BJAC'        - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks and L1 correction for off-diag blocks
 !
@@ -70,12 +70,12 @@
 !                       applied as post-smoother at each level, but the
 !                       coarsest one; four sweeps of the block-Jacobi solver,
 !                       with LU from UMFPACK on the blocks, are applied at
-!                       the coarsest level, on the distributed coarse matrix. 
+!                       the coarsest level, on the distributed coarse matrix.
 !                       The smoothed aggregation algorithm with threshold 0
 !                       is used to build the coarse matrix.
 !
 !  For the multilevel preconditioners, the levels are numbered in increasing
-!  order starting from the finest one, i.e. level 1 is the finest level. 
+!  order starting from the finest one, i.e. level 1 is the finest level.
 !
 !
 ! Arguments:
@@ -87,7 +87,7 @@
 !               lowercase strings).
 !    info    -  integer, output.
 !               Error code.
-!  
+!
 subroutine amg_cprecinit(ctxt,prec,ptype,info)
 
   use psb_base_mod
@@ -113,98 +113,105 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
 
   ! Local variables
   integer(psb_ipk_)            :: nlev_, ilev_
+  integer(psb_ipk_)            :: err_act
+  integer(psb_ipk_)            :: debug_level, debug_unit
   real(psb_spk_)               :: thr
   character(len=*), parameter  :: name='amg_precinit'
   info = psb_success_
+  call psb_erractionsave(err_act)
+  if (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_; goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
 
-  if (allocated(prec%precv)) then 
-    call prec%free(info) 
-    if (info /= psb_success_) then 
-      ! Do we want to do something? 
+  if (allocated(prec%precv)) then
+    call prec%free(info)
+    if (info /= psb_success_) then
+      ! Do we want to do something?
     endif
   endif
   prec%ctxt = ctxt
-  prec%ag_data%min_coarse_size             = -1
-  prec%ag_data%min_coarse_size_per_process = -1
+  call prec%ag_data%default()
 
   select case(psb_toupper(trim(ptype)))
-  case ('NOPREC','NONE') 
+  case ('NOPREC','NONE')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_base_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_base_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('JAC','DIAG','JACOBI') 
+  case ('JAC','DIAG','JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') 
+  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('GS','FWGS') 
+  case ('GS','FWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('BWGS') 
+  case ('BWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('FBGS') 
+  case ('FBGS')
     nlev_ = 1
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
     call prec%precv(ilev_)%default()
 
-  case ('BJAC') 
+  case ('BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-BJAC','L1_BJAC') 
+  case ('L1-BJAC','L1_BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
   case ('AS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_c_as_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_c_as_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_c_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
 
@@ -213,20 +220,24 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
     nlev_ = prec%ag_data%max_levs
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
+    if (info /= psb_success_ ) then
+      call psb_errpush(info,name,a_err='Error from hierarchy init')
+      goto 9999
+    endif
 
-    do ilev_ = 1, nlev_  
+    do ilev_ = 1, nlev_
       call prec%precv(ilev_)%default()
     end do
     call prec%set('ML_CYCLE','VCYCLE',info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
 #if  defined(HAVE_MUMPS_)
-    call prec%set('COARSE_SOLVE','MUMPS',info)    
+    call prec%set('COARSE_SOLVE','MUMPS',info)
 #elif defined(HAVE_SLU_)
     call prec%set('COARSE_SOLVE','SLU',info)
 #else
     call prec%set('COARSE_SOLVE','ILU',info)
 #endif
-    
+
   case default
     write(psb_err_unit,*) name,&
          &': Warning: Unknown preconditioner type request "',ptype,'"'
@@ -234,5 +245,10 @@ subroutine amg_cprecinit(ctxt,prec,ptype,info)
 
   end select
 
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
 
 end subroutine amg_cprecinit

amgprec/impl/amg_cprecset.F90

@@ -45,7 +45,7 @@ subroutine amg_cprecsetsm(p,val,info,ilev,ilmax,pos)
   implicit none
 
   ! Arguments
-  class(amg_cprec_type), intent(inout)         :: p
+  class(amg_cprec_type), target, intent(inout):: p
   class(amg_c_base_smoother_type), intent(in) :: val
   integer(psb_ipk_), intent(out)              :: info
   integer(psb_ipk_), optional, intent(in)     :: ilev,ilmax

amgprec/impl/amg_dcprecset.F90

@@ -474,6 +474,16 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
           select case (psb_toupper(string))
           case('BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_UMF_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_umf_,info,pos=pos)
+#elif defined(HAVE_SLU_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+            call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+            call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
           case('L1-BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_UMF_)
@@ -638,6 +648,16 @@ subroutine amg_dcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
         select case (psb_toupper(trim(string)))
         case('BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_UMF_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_umf_,info,pos=pos)
+#elif defined(HAVE_SLU_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+          call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+          call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
         case('L1-BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_UMF_)

amgprec/impl/amg_dfile_prec_descr.f90

@@ -65,7 +65,7 @@
 !         0: normal
 !        >1: increased details 
 !
-subroutine amg_dfile_prec_descr(prec,iout,root, verbosity)
+subroutine amg_dfile_prec_descr(prec,info,iout,root, verbosity)
   use psb_base_mod
   use amg_d_prec_mod, amg_protect_name => amg_dfile_prec_descr
   use amg_d_inner_mod
@@ -74,13 +74,14 @@ subroutine amg_dfile_prec_descr(prec,iout,root, verbosity)
   implicit none 
   ! Arguments
   class(amg_dprec_type), intent(in)      :: prec
+  integer(psb_ipk_), intent(out) :: info
   integer(psb_ipk_), intent(in), optional :: iout
   integer(psb_ipk_), intent(in), optional :: root
   integer(psb_ipk_), intent(in), optional   :: verbosity
 
 
   ! Local variables
-  integer(psb_ipk_)   :: ilev, nlev, ilmin, info, nswps
+  integer(psb_ipk_)   :: ilev, nlev, ilmin, nswps
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   logical             :: is_symgs

amgprec/impl/amg_dprecinit.F90

@@ -1,15 +1,15 @@
-!   
-!   
+!
+!
 !                             AMG4PSBLAS version 1.0
 !    Algebraic Multigrid Package
 !               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2021 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Fabio Durastante
+!
 !    Redistribution and use in source and binary forms, with or without
 !    modification, are permitted provided that the following conditions
 !    are met:
@@ -21,7 +21,7 @@
 !      3. The name of the AMG4PSBLAS group or the names of its contributors may
 !         not be used to endorse or promote products derived from this
 !         software without specific written permission.
-!   
+!
 !    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@@ -33,8 +33,8 @@
 !    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
+!
+!
 ! File: amg_dprecinit.f90
 !
 ! Subroutine: amg_dprecinit
@@ -42,21 +42,21 @@
 !
 !  This routine allocates and initializes the preconditioner data structure,
 !  according to the preconditioner type chosen by the user.
-!  
+!
 !  A default preconditioner is set for each preconditioner type
 !  specified by the user:
 !
 !    'NOPREC'         - no preconditioner
 !
-!    'DIAG', 'JACOBI' - diagonal/Jacobi           
+!    'DIAG', 'JACOBI' - diagonal/Jacobi
 !
 !    'L1-DIAG', 'L1-JACOBI' - diagonal/Jacobi  with L1 norm correction
 !
 !    'GS', 'FBGS'     - Hybrid Gauss-Seidel, also symmetrized
-!                       
+!
 !    'BJAC'           - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks
-!                       
+!
 !    'L1-BJAC'        - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks and L1 correction for off-diag blocks
 !
@@ -70,12 +70,12 @@
 !                       applied as post-smoother at each level, but the
 !                       coarsest one; four sweeps of the block-Jacobi solver,
 !                       with LU from UMFPACK on the blocks, are applied at
-!                       the coarsest level, on the distributed coarse matrix. 
+!                       the coarsest level, on the distributed coarse matrix.
 !                       The smoothed aggregation algorithm with threshold 0
 !                       is used to build the coarse matrix.
 !
 !  For the multilevel preconditioners, the levels are numbered in increasing
-!  order starting from the finest one, i.e. level 1 is the finest level. 
+!  order starting from the finest one, i.e. level 1 is the finest level.
 !
 !
 ! Arguments:
@@ -87,7 +87,7 @@
 !               lowercase strings).
 !    info    -  integer, output.
 !               Error code.
-!  
+!
 subroutine amg_dprecinit(ctxt,prec,ptype,info)
 
   use psb_base_mod
@@ -116,98 +116,105 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
 
   ! Local variables
   integer(psb_ipk_)            :: nlev_, ilev_
+  integer(psb_ipk_)            :: err_act
+  integer(psb_ipk_)            :: debug_level, debug_unit
   real(psb_dpk_)               :: thr
   character(len=*), parameter  :: name='amg_precinit'
   info = psb_success_
+  call psb_erractionsave(err_act)
+  if (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_; goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
 
-  if (allocated(prec%precv)) then 
-    call prec%free(info) 
-    if (info /= psb_success_) then 
-      ! Do we want to do something? 
+  if (allocated(prec%precv)) then
+    call prec%free(info)
+    if (info /= psb_success_) then
+      ! Do we want to do something?
     endif
   endif
   prec%ctxt = ctxt
-  prec%ag_data%min_coarse_size             = -1
-  prec%ag_data%min_coarse_size_per_process = -1
+  call prec%ag_data%default()
 
   select case(psb_toupper(trim(ptype)))
-  case ('NOPREC','NONE') 
+  case ('NOPREC','NONE')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_base_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_base_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('JAC','DIAG','JACOBI') 
+  case ('JAC','DIAG','JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') 
+  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('GS','FWGS') 
+  case ('GS','FWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('BWGS') 
+  case ('BWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('FBGS') 
+  case ('FBGS')
     nlev_ = 1
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
     call prec%precv(ilev_)%default()
 
-  case ('BJAC') 
+  case ('BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-BJAC','L1_BJAC') 
+  case ('L1-BJAC','L1_BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
   case ('AS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_d_as_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_d_as_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_d_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
 
@@ -216,22 +223,26 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
     nlev_ = prec%ag_data%max_levs
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
+    if (info /= psb_success_ ) then
+      call psb_errpush(info,name,a_err='Error from hierarchy init')
+      goto 9999
+    endif
 
-    do ilev_ = 1, nlev_  
+    do ilev_ = 1, nlev_
       call prec%precv(ilev_)%default()
     end do
     call prec%set('ML_CYCLE','VCYCLE',info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
-#if defined(HAVE_UMF_) 
+#if defined(HAVE_UMF_)
     call prec%set('COARSE_SOLVE','UMF',info)
 #elif defined(HAVE_MUMPS_)
-    call prec%set('COARSE_SOLVE','MUMPS',info)    
+    call prec%set('COARSE_SOLVE','MUMPS',info)
 #elif defined(HAVE_SLU_)
     call prec%set('COARSE_SOLVE','SLU',info)
 #else
     call prec%set('COARSE_SOLVE','ILU',info)
 #endif
-    
+
   case default
     write(psb_err_unit,*) name,&
          &': Warning: Unknown preconditioner type request "',ptype,'"'
@@ -239,5 +250,10 @@ subroutine amg_dprecinit(ctxt,prec,ptype,info)
 
   end select
 
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
 
 end subroutine amg_dprecinit

amgprec/impl/amg_dprecset.F90

@@ -45,7 +45,7 @@ subroutine amg_dprecsetsm(p,val,info,ilev,ilmax,pos)
   implicit none
 
   ! Arguments
-  class(amg_dprec_type), intent(inout)         :: p
+  class(amg_dprec_type), target, intent(inout):: p
   class(amg_d_base_smoother_type), intent(in) :: val
   integer(psb_ipk_), intent(out)              :: info
   integer(psb_ipk_), optional, intent(in)     :: ilev,ilmax

amgprec/impl/amg_scprecset.F90

@@ -462,6 +462,14 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
           select case (psb_toupper(string))
           case('BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_SLU_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+            call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+            call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
           case('L1-BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_SLU_)
@@ -612,6 +620,14 @@ subroutine amg_scprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
         select case (psb_toupper(trim(string)))
         case('BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_SLU_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+          call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+          call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
         case('L1-BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_SLU_)

amgprec/impl/amg_sfile_prec_descr.f90

@@ -65,7 +65,7 @@
 !         0: normal
 !        >1: increased details 
 !
-subroutine amg_sfile_prec_descr(prec,iout,root, verbosity)
+subroutine amg_sfile_prec_descr(prec,info,iout,root, verbosity)
   use psb_base_mod
   use amg_s_prec_mod, amg_protect_name => amg_sfile_prec_descr
   use amg_s_inner_mod
@@ -74,13 +74,14 @@ subroutine amg_sfile_prec_descr(prec,iout,root, verbosity)
   implicit none 
   ! Arguments
   class(amg_sprec_type), intent(in)      :: prec
+  integer(psb_ipk_), intent(out) :: info
   integer(psb_ipk_), intent(in), optional :: iout
   integer(psb_ipk_), intent(in), optional :: root
   integer(psb_ipk_), intent(in), optional   :: verbosity
 
 
   ! Local variables
-  integer(psb_ipk_)   :: ilev, nlev, ilmin, info, nswps
+  integer(psb_ipk_)   :: ilev, nlev, ilmin, nswps
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   logical             :: is_symgs

amgprec/impl/amg_sprecinit.F90

@@ -1,15 +1,15 @@
-!   
-!   
+!
+!
 !                             AMG4PSBLAS version 1.0
 !    Algebraic Multigrid Package
 !               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2021 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Fabio Durastante
+!
 !    Redistribution and use in source and binary forms, with or without
 !    modification, are permitted provided that the following conditions
 !    are met:
@@ -21,7 +21,7 @@
 !      3. The name of the AMG4PSBLAS group or the names of its contributors may
 !         not be used to endorse or promote products derived from this
 !         software without specific written permission.
-!   
+!
 !    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@@ -33,8 +33,8 @@
 !    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
+!
+!
 ! File: amg_sprecinit.f90
 !
 ! Subroutine: amg_sprecinit
@@ -42,21 +42,21 @@
 !
 !  This routine allocates and initializes the preconditioner data structure,
 !  according to the preconditioner type chosen by the user.
-!  
+!
 !  A default preconditioner is set for each preconditioner type
 !  specified by the user:
 !
 !    'NOPREC'         - no preconditioner
 !
-!    'DIAG', 'JACOBI' - diagonal/Jacobi           
+!    'DIAG', 'JACOBI' - diagonal/Jacobi
 !
 !    'L1-DIAG', 'L1-JACOBI' - diagonal/Jacobi  with L1 norm correction
 !
 !    'GS', 'FBGS'     - Hybrid Gauss-Seidel, also symmetrized
-!                       
+!
 !    'BJAC'           - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks
-!                       
+!
 !    'L1-BJAC'        - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks and L1 correction for off-diag blocks
 !
@@ -70,12 +70,12 @@
 !                       applied as post-smoother at each level, but the
 !                       coarsest one; four sweeps of the block-Jacobi solver,
 !                       with LU from UMFPACK on the blocks, are applied at
-!                       the coarsest level, on the distributed coarse matrix. 
+!                       the coarsest level, on the distributed coarse matrix.
 !                       The smoothed aggregation algorithm with threshold 0
 !                       is used to build the coarse matrix.
 !
 !  For the multilevel preconditioners, the levels are numbered in increasing
-!  order starting from the finest one, i.e. level 1 is the finest level. 
+!  order starting from the finest one, i.e. level 1 is the finest level.
 !
 !
 ! Arguments:
@@ -87,7 +87,7 @@
 !               lowercase strings).
 !    info    -  integer, output.
 !               Error code.
-!  
+!
 subroutine amg_sprecinit(ctxt,prec,ptype,info)
 
   use psb_base_mod
@@ -113,98 +113,105 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
 
   ! Local variables
   integer(psb_ipk_)            :: nlev_, ilev_
+  integer(psb_ipk_)            :: err_act
+  integer(psb_ipk_)            :: debug_level, debug_unit
   real(psb_spk_)               :: thr
   character(len=*), parameter  :: name='amg_precinit'
   info = psb_success_
+  call psb_erractionsave(err_act)
+  if (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_; goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
 
-  if (allocated(prec%precv)) then 
-    call prec%free(info) 
-    if (info /= psb_success_) then 
-      ! Do we want to do something? 
+  if (allocated(prec%precv)) then
+    call prec%free(info)
+    if (info /= psb_success_) then
+      ! Do we want to do something?
     endif
   endif
   prec%ctxt = ctxt
-  prec%ag_data%min_coarse_size             = -1
-  prec%ag_data%min_coarse_size_per_process = -1
+  call prec%ag_data%default()
 
   select case(psb_toupper(trim(ptype)))
-  case ('NOPREC','NONE') 
+  case ('NOPREC','NONE')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_base_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_base_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('JAC','DIAG','JACOBI') 
+  case ('JAC','DIAG','JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') 
+  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('GS','FWGS') 
+  case ('GS','FWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('BWGS') 
+  case ('BWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('FBGS') 
+  case ('FBGS')
     nlev_ = 1
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
     call prec%precv(ilev_)%default()
 
-  case ('BJAC') 
+  case ('BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-BJAC','L1_BJAC') 
+  case ('L1-BJAC','L1_BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
   case ('AS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_s_as_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_s_as_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_s_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
 
@@ -213,20 +220,24 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
     nlev_ = prec%ag_data%max_levs
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
+    if (info /= psb_success_ ) then
+      call psb_errpush(info,name,a_err='Error from hierarchy init')
+      goto 9999
+    endif
 
-    do ilev_ = 1, nlev_  
+    do ilev_ = 1, nlev_
       call prec%precv(ilev_)%default()
     end do
     call prec%set('ML_CYCLE','VCYCLE',info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
 #if  defined(HAVE_MUMPS_)
-    call prec%set('COARSE_SOLVE','MUMPS',info)    
+    call prec%set('COARSE_SOLVE','MUMPS',info)
 #elif defined(HAVE_SLU_)
     call prec%set('COARSE_SOLVE','SLU',info)
 #else
     call prec%set('COARSE_SOLVE','ILU',info)
 #endif
-    
+
   case default
     write(psb_err_unit,*) name,&
          &': Warning: Unknown preconditioner type request "',ptype,'"'
@@ -234,5 +245,10 @@ subroutine amg_sprecinit(ctxt,prec,ptype,info)
 
   end select
 
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
 
 end subroutine amg_sprecinit

amgprec/impl/amg_sprecset.F90

@@ -45,7 +45,7 @@ subroutine amg_sprecsetsm(p,val,info,ilev,ilmax,pos)
   implicit none
 
   ! Arguments
-  class(amg_sprec_type), intent(inout)         :: p
+  class(amg_sprec_type), target, intent(inout):: p
   class(amg_s_base_smoother_type), intent(in) :: val
   integer(psb_ipk_), intent(out)              :: info
   integer(psb_ipk_), optional, intent(in)     :: ilev,ilmax

amgprec/impl/amg_zcprecset.F90

@@ -474,6 +474,16 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
           select case (psb_toupper(string))
           case('BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_UMF_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_umf_,info,pos=pos)
+#elif defined(HAVE_SLU_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+            call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+            call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+            call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
           case('L1-BJAC')
             call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_UMF_)
@@ -638,6 +648,16 @@ subroutine amg_zcprecsetc(p,what,string,info,ilev,ilmax,pos,idx)
         select case (psb_toupper(trim(string)))
         case('BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_bjac_,info,pos=pos)
+#if defined(HAVE_UMF_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_umf_,info,pos=pos)
+#elif defined(HAVE_SLU_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_slu_,info,pos=pos)
+#elif defined(HAVE_MUMPS_)
+          call p%precv(nlev_)%set('SUB_SOLVE',amg_mumps_,info,pos=pos)
+#else
+          call p%precv(nlev_)%set('SUB_SOLVE',psb_ilu_n_,info,pos=pos)
+#endif
+          call p%precv(nlev_)%set('COARSE_MAT',amg_distr_mat_,info)
         case('L1-BJAC')
           call p%precv(nlev_)%set('SMOOTHER_TYPE',amg_l1_bjac_,info,pos=pos)
 #if defined(HAVE_UMF_)

amgprec/impl/amg_zfile_prec_descr.f90

@@ -65,7 +65,7 @@
 !         0: normal
 !        >1: increased details 
 !
-subroutine amg_zfile_prec_descr(prec,iout,root, verbosity)
+subroutine amg_zfile_prec_descr(prec,info,iout,root, verbosity)
   use psb_base_mod
   use amg_z_prec_mod, amg_protect_name => amg_zfile_prec_descr
   use amg_z_inner_mod
@@ -74,13 +74,14 @@ subroutine amg_zfile_prec_descr(prec,iout,root, verbosity)
   implicit none 
   ! Arguments
   class(amg_zprec_type), intent(in)      :: prec
+  integer(psb_ipk_), intent(out) :: info
   integer(psb_ipk_), intent(in), optional :: iout
   integer(psb_ipk_), intent(in), optional :: root
   integer(psb_ipk_), intent(in), optional   :: verbosity
 
 
   ! Local variables
-  integer(psb_ipk_)   :: ilev, nlev, ilmin, info, nswps
+  integer(psb_ipk_)   :: ilev, nlev, ilmin, nswps
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_)   :: me, np
   logical             :: is_symgs

amgprec/impl/amg_zprecinit.F90

@@ -1,15 +1,15 @@
-!   
-!   
+!
+!
 !                             AMG4PSBLAS version 1.0
 !    Algebraic Multigrid Package
 !               based on PSBLAS (Parallel Sparse BLAS version 3.7)
-!    
-!    (C) Copyright 2021 
-!  
-!        Salvatore Filippone  
-!        Pasqua D'Ambra   
-!        Fabio Durastante        
-!   
+!
+!    (C) Copyright 2021
+!
+!        Salvatore Filippone
+!        Pasqua D'Ambra
+!        Fabio Durastante
+!
 !    Redistribution and use in source and binary forms, with or without
 !    modification, are permitted provided that the following conditions
 !    are met:
@@ -21,7 +21,7 @@
 !      3. The name of the AMG4PSBLAS group or the names of its contributors may
 !         not be used to endorse or promote products derived from this
 !         software without specific written permission.
-!   
+!
 !    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@@ -33,8 +33,8 @@
 !    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
-!   
-!  
+!
+!
 ! File: amg_zprecinit.f90
 !
 ! Subroutine: amg_zprecinit
@@ -42,21 +42,21 @@
 !
 !  This routine allocates and initializes the preconditioner data structure,
 !  according to the preconditioner type chosen by the user.
-!  
+!
 !  A default preconditioner is set for each preconditioner type
 !  specified by the user:
 !
 !    'NOPREC'         - no preconditioner
 !
-!    'DIAG', 'JACOBI' - diagonal/Jacobi           
+!    'DIAG', 'JACOBI' - diagonal/Jacobi
 !
 !    'L1-DIAG', 'L1-JACOBI' - diagonal/Jacobi  with L1 norm correction
 !
 !    'GS', 'FBGS'     - Hybrid Gauss-Seidel, also symmetrized
-!                       
+!
 !    'BJAC'           - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks
-!                       
+!
 !    'L1-BJAC'        - block Jacobi preconditioner, with ILU(0)
 !                       on the local blocks and L1 correction for off-diag blocks
 !
@@ -70,12 +70,12 @@
 !                       applied as post-smoother at each level, but the
 !                       coarsest one; four sweeps of the block-Jacobi solver,
 !                       with LU from UMFPACK on the blocks, are applied at
-!                       the coarsest level, on the distributed coarse matrix. 
+!                       the coarsest level, on the distributed coarse matrix.
 !                       The smoothed aggregation algorithm with threshold 0
 !                       is used to build the coarse matrix.
 !
 !  For the multilevel preconditioners, the levels are numbered in increasing
-!  order starting from the finest one, i.e. level 1 is the finest level. 
+!  order starting from the finest one, i.e. level 1 is the finest level.
 !
 !
 ! Arguments:
@@ -87,7 +87,7 @@
 !               lowercase strings).
 !    info    -  integer, output.
 !               Error code.
-!  
+!
 subroutine amg_zprecinit(ctxt,prec,ptype,info)
 
   use psb_base_mod
@@ -116,98 +116,105 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
 
   ! Local variables
   integer(psb_ipk_)            :: nlev_, ilev_
+  integer(psb_ipk_)            :: err_act
+  integer(psb_ipk_)            :: debug_level, debug_unit
   real(psb_dpk_)               :: thr
   character(len=*), parameter  :: name='amg_precinit'
   info = psb_success_
+  call psb_erractionsave(err_act)
+  if (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_; goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
 
-  if (allocated(prec%precv)) then 
-    call prec%free(info) 
-    if (info /= psb_success_) then 
-      ! Do we want to do something? 
+  if (allocated(prec%precv)) then
+    call prec%free(info)
+    if (info /= psb_success_) then
+      ! Do we want to do something?
     endif
   endif
   prec%ctxt = ctxt
-  prec%ag_data%min_coarse_size             = -1
-  prec%ag_data%min_coarse_size_per_process = -1
+  call prec%ag_data%default()
 
   select case(psb_toupper(trim(ptype)))
-  case ('NOPREC','NONE') 
+  case ('NOPREC','NONE')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_base_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_base_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_id_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('JAC','DIAG','JACOBI') 
+  case ('JAC','DIAG','JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI') 
+  case ('L1-DIAG','L1-JACOBI','L1_DIAG','L1_JACOBI')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_l1_diag_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('GS','FWGS') 
+  case ('GS','FWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_gs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('BWGS') 
+  case ('BWGS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_bwgs_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('FBGS') 
+  case ('FBGS')
     nlev_ = 1
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
     call prec%precv(ilev_)%default()
 
-  case ('BJAC') 
+  case ('BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
-  case ('L1-BJAC','L1_BJAC') 
+  case ('L1-BJAC','L1_BJAC')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_l1_jac_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
   case ('AS')
     nlev_ = 1
     ilev_ = 1
-    allocate(prec%precv(nlev_),stat=info) 
-    allocate(amg_z_as_smoother_type :: prec%precv(ilev_)%sm, stat=info) 
+    allocate(prec%precv(nlev_),stat=info)
+    allocate(amg_z_as_smoother_type :: prec%precv(ilev_)%sm, stat=info)
     if (info /= psb_success_) return
-    allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info) 
+    allocate(amg_z_ilu_solver_type :: prec%precv(ilev_)%sm%sv, stat=info)
     call prec%precv(ilev_)%default()
 
 
@@ -216,22 +223,26 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
     nlev_ = prec%ag_data%max_levs
     ilev_ = 1
     allocate(prec%precv(nlev_),stat=info)
+    if (info /= psb_success_ ) then
+      call psb_errpush(info,name,a_err='Error from hierarchy init')
+      goto 9999
+    endif
 
-    do ilev_ = 1, nlev_  
+    do ilev_ = 1, nlev_
       call prec%precv(ilev_)%default()
     end do
     call prec%set('ML_CYCLE','VCYCLE',info)
     call prec%set('SMOOTHER_TYPE','FBGS',info)
-#if defined(HAVE_UMF_) 
+#if defined(HAVE_UMF_)
     call prec%set('COARSE_SOLVE','UMF',info)
 #elif defined(HAVE_MUMPS_)
-    call prec%set('COARSE_SOLVE','MUMPS',info)    
+    call prec%set('COARSE_SOLVE','MUMPS',info)
 #elif defined(HAVE_SLU_)
     call prec%set('COARSE_SOLVE','SLU',info)
 #else
     call prec%set('COARSE_SOLVE','ILU',info)
 #endif
-    
+
   case default
     write(psb_err_unit,*) name,&
          &': Warning: Unknown preconditioner type request "',ptype,'"'
@@ -239,5 +250,10 @@ subroutine amg_zprecinit(ctxt,prec,ptype,info)
 
   end select
 
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+  return
 
 end subroutine amg_zprecinit

amgprec/impl/amg_zprecset.F90

@@ -45,7 +45,7 @@ subroutine amg_zprecsetsm(p,val,info,ilev,ilmax,pos)
   implicit none
 
   ! Arguments
-  class(amg_zprec_type), intent(inout)         :: p
+  class(amg_zprec_type), target, intent(inout):: p
   class(amg_z_base_smoother_type), intent(in) :: val
   integer(psb_ipk_), intent(out)              :: info
   integer(psb_ipk_), optional, intent(in)     :: ilev,ilmax

amgprec/impl/smoother/amg_c_as_smoother_free.f90

@@ -61,6 +61,7 @@ subroutine amg_c_as_smoother_free(sm,info)
     end if
   end if
   call sm%nd%free()
+  call sm%desc_data%free(info)
 
   call psb_erractionrestore(err_act)
   return

amgprec/impl/smoother/amg_d_as_smoother_free.f90

@@ -61,6 +61,7 @@ subroutine amg_d_as_smoother_free(sm,info)
     end if
   end if
   call sm%nd%free()
+  call sm%desc_data%free(info)
 
   call psb_erractionrestore(err_act)
   return

amgprec/impl/smoother/amg_s_as_smoother_free.f90

@@ -61,6 +61,7 @@ subroutine amg_s_as_smoother_free(sm,info)
     end if
   end if
   call sm%nd%free()
+  call sm%desc_data%free(info)
 
   call psb_erractionrestore(err_act)
   return

amgprec/impl/smoother/amg_z_as_smoother_free.f90

@@ -61,6 +61,7 @@ subroutine amg_z_as_smoother_free(sm,info)
     end if
   end if
   call sm%nd%free()
+  call sm%desc_data%free(info)
 
   call psb_erractionrestore(err_act)
   return

amgprec/impl/solver/amg_c_base_solver_csetr.f90

@@ -42,7 +42,7 @@ subroutine amg_c_base_solver_csetr(sv,what,val,info,idx)
   Implicit None
   ! Arguments
   class(amg_c_base_solver_type), intent(inout) :: sv 
-  integer(psb_ipk_), intent(in)                  :: what 
+  character(len=*), intent(in)                   :: what 
   real(psb_spk_), intent(in)                      :: val
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_), intent(in), optional        :: idx

amgprec/impl/solver/amg_c_bwgs_solver_bld.f90

@@ -44,7 +44,7 @@ subroutine amg_c_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   ! Arguments
   type(psb_cspmat_type), intent(in), target           :: a
-  Type(psb_desc_type), Intent(in)                     :: desc_a 
+  Type(psb_desc_type), Intent(inout)                  :: desc_a 
   class(amg_c_bwgs_solver_type), intent(inout)         :: sv
   integer(psb_ipk_), intent(out)                      :: info
   type(psb_cspmat_type), intent(in), target, optional :: b

amgprec/impl/solver/amg_c_krm_solver_impl.f90

@@ -77,7 +77,7 @@
 !  This is the implementation file corresponding to amg_c_krm_solver_mod.
 !
 !
-subroutine amg_c_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
+subroutine amg_c_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   use psb_base_mod
   use amg_c_krm_solver, amg_protect_name => amg_c_krm_solver_bld
@@ -85,13 +85,14 @@ subroutine amg_c_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
   Implicit None
 
   ! Arguments
-  type(psb_cspmat_type), intent(inout), target        :: a
+  type(psb_cspmat_type), intent(in), 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
   type(psb_cspmat_type), intent(in), target, optional :: b
   class(psb_c_base_sparse_mat), intent(in), optional  :: amold
   class(psb_c_base_vect_type), intent(in), optional   :: vmold
+  class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
   integer(psb_lpk_) :: lnr

amgprec/impl/solver/amg_d_base_solver_csetr.f90

@@ -42,7 +42,7 @@ subroutine amg_d_base_solver_csetr(sv,what,val,info,idx)
   Implicit None
   ! Arguments
   class(amg_d_base_solver_type), intent(inout) :: sv 
-  integer(psb_ipk_), intent(in)                  :: what 
+  character(len=*), intent(in)                   :: what 
   real(psb_dpk_), intent(in)                      :: val
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_), intent(in), optional        :: idx

amgprec/impl/solver/amg_d_bwgs_solver_bld.f90

@@ -44,7 +44,7 @@ subroutine amg_d_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   ! Arguments
   type(psb_dspmat_type), intent(in), target           :: a
-  Type(psb_desc_type), Intent(in)                     :: desc_a 
+  Type(psb_desc_type), Intent(inout)                  :: desc_a 
   class(amg_d_bwgs_solver_type), intent(inout)         :: sv
   integer(psb_ipk_), intent(out)                      :: info
   type(psb_dspmat_type), intent(in), target, optional :: b

amgprec/impl/solver/amg_d_krm_solver_impl.f90

@@ -77,7 +77,7 @@
 !  This is the implementation file corresponding to amg_d_krm_solver_mod.
 !
 !
-subroutine amg_d_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
+subroutine amg_d_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   use psb_base_mod
   use amg_d_krm_solver, amg_protect_name => amg_d_krm_solver_bld
@@ -85,13 +85,14 @@ subroutine amg_d_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
   Implicit None
 
   ! Arguments
-  type(psb_dspmat_type), intent(inout), target        :: a
+  type(psb_dspmat_type), intent(in), 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
   type(psb_dspmat_type), intent(in), target, optional :: b
   class(psb_d_base_sparse_mat), intent(in), optional  :: amold
   class(psb_d_base_vect_type), intent(in), optional   :: vmold
+  class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
   integer(psb_lpk_) :: lnr

amgprec/impl/solver/amg_s_base_solver_csetr.f90

@@ -42,7 +42,7 @@ subroutine amg_s_base_solver_csetr(sv,what,val,info,idx)
   Implicit None
   ! Arguments
   class(amg_s_base_solver_type), intent(inout) :: sv 
-  integer(psb_ipk_), intent(in)                  :: what 
+  character(len=*), intent(in)                   :: what 
   real(psb_spk_), intent(in)                      :: val
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_), intent(in), optional        :: idx

amgprec/impl/solver/amg_s_bwgs_solver_bld.f90

@@ -44,7 +44,7 @@ subroutine amg_s_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   ! Arguments
   type(psb_sspmat_type), intent(in), target           :: a
-  Type(psb_desc_type), Intent(in)                     :: desc_a 
+  Type(psb_desc_type), Intent(inout)                  :: desc_a 
   class(amg_s_bwgs_solver_type), intent(inout)         :: sv
   integer(psb_ipk_), intent(out)                      :: info
   type(psb_sspmat_type), intent(in), target, optional :: b

amgprec/impl/solver/amg_s_krm_solver_impl.f90

@@ -77,7 +77,7 @@
 !  This is the implementation file corresponding to amg_s_krm_solver_mod.
 !
 !
-subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
+subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   use psb_base_mod
   use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_bld
@@ -85,13 +85,14 @@ subroutine amg_s_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
   Implicit None
 
   ! Arguments
-  type(psb_sspmat_type), intent(inout), target        :: a
+  type(psb_sspmat_type), intent(in), 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
   type(psb_sspmat_type), intent(in), target, optional :: b
   class(psb_s_base_sparse_mat), intent(in), optional  :: amold
   class(psb_s_base_vect_type), intent(in), optional   :: vmold
+  class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
   integer(psb_lpk_) :: lnr

amgprec/impl/solver/amg_z_base_solver_csetr.f90

@@ -42,7 +42,7 @@ subroutine amg_z_base_solver_csetr(sv,what,val,info,idx)
   Implicit None
   ! Arguments
   class(amg_z_base_solver_type), intent(inout) :: sv 
-  integer(psb_ipk_), intent(in)                  :: what 
+  character(len=*), intent(in)                   :: what 
   real(psb_dpk_), intent(in)                      :: val
   integer(psb_ipk_), intent(out)                 :: info
   integer(psb_ipk_), intent(in), optional        :: idx

amgprec/impl/solver/amg_z_bwgs_solver_bld.f90

@@ -44,7 +44,7 @@ subroutine amg_z_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   ! Arguments
   type(psb_zspmat_type), intent(in), target           :: a
-  Type(psb_desc_type), Intent(in)                     :: desc_a 
+  Type(psb_desc_type), Intent(inout)                  :: desc_a 
   class(amg_z_bwgs_solver_type), intent(inout)         :: sv
   integer(psb_ipk_), intent(out)                      :: info
   type(psb_zspmat_type), intent(in), target, optional :: b

amgprec/impl/solver/amg_z_krm_solver_impl.f90

@@ -77,7 +77,7 @@
 !  This is the implementation file corresponding to amg_z_krm_solver_mod.
 !
 !
-subroutine amg_z_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
+subroutine amg_z_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
 
   use psb_base_mod
   use amg_z_krm_solver, amg_protect_name => amg_z_krm_solver_bld
@@ -85,13 +85,14 @@ subroutine amg_z_krm_solver_bld(a,desc_a,sv,info,b,amold,vmold)
   Implicit None
 
   ! Arguments
-  type(psb_zspmat_type), intent(inout), target        :: a
+  type(psb_zspmat_type), intent(in), 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
   type(psb_zspmat_type), intent(in), target, optional :: b
   class(psb_z_base_sparse_mat), intent(in), optional  :: amold
   class(psb_z_base_vect_type), intent(in), optional   :: vmold
+  class(psb_i_base_vect_type), intent(in), optional   :: imold
   ! Local variables
   integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
   integer(psb_lpk_) :: lnr

cbind/amgprec/amg_dprec_cbind_mod.F90

@@ -384,7 +384,7 @@ contains
 
    integer(psb_c_ipk_) :: res
    type(psb_c_object_type) :: ph
-   integer               :: info
+   integer(psb_c_ipk_)     :: info
    type(amg_dprec_type), pointer :: precp
 
    res = -1
@@ -396,7 +396,7 @@ contains
    end if
 
 
-   call precp%descr()
+   call precp%descr(info)
    call flush(psb_out_unit)
 
    info = 0

cbind/amgprec/amg_zprec_cbind_mod.F90

@@ -384,7 +384,7 @@ contains
 
    integer(psb_c_ipk_) :: res
    type(psb_c_object_type) :: ph
-   integer               :: info
+   integer(psb_c_ipk_)     :: info
    type(amg_zprec_type), pointer :: precp
 
    res = -1
@@ -396,7 +396,7 @@ contains
    end if
 
 
-   call precp%descr()
+   call precp%descr(info)
    call flush(psb_out_unit)
 
    info = 0

docs/html/index.html

@@ -31,7 +31,7 @@ class="cmr-12">University of Rome Tor-Vergata and IAC-CNR</span><br
 class="newline" /> <span 
 class="cmr-12">Software version: 1.0</span><br 
 class="newline" /><span 
-class="cmr-12">April 12, 2021</span>
+class="cmr-12">May 11th, 2021</span>
                                                                                
 
                                                                                

docs/html/userhtml.css

@@ -149,6 +149,7 @@ div.abstract {width:100%;}
 .Ovalbox-thick { padding-left:3pt; padding-right:3pt; border:solid thick; }
 .shadowbox { padding-left:3pt; padding-right:3pt; border:solid thin; border-right:solid thick; border-bottom:solid thick; }
 .doublebox { padding-left:3pt; padding-right:3pt; border-style:double; border:solid thick; }
+.rotatebox{display: inline-block;}
 .figure img.graphics {margin-left:10%;}
 .lstlisting .label{margin-right:0.5em; }
 div.lstlisting{font-family: monospace; white-space: nowrap; margin-top:0.5em; margin-bottom:0.5em; }

docs/html/userhtml.html

@@ -31,7 +31,7 @@ class="cmr-12">University of Rome Tor-Vergata and IAC-CNR</span><br
 class="newline" /> <span 
 class="cmr-12">Software version: 1.0</span><br 
 class="newline" /><span 
-class="cmr-12">April 12, 2021</span>
+class="cmr-12">May 11th, 2021</span>
                                                                                
 
                                                                                

docs/html/userhtmlli5.html

@@ -211,7 +211,7 @@ class="cmr-12">&#x00A0;Pothen,  </span><span
 class="cmti-12">Distributed-memory parallel algorithms for matching and</span>
      <span 
 class="cmti-12">coloring</span><span 
-class="cmr-12">, in PCO11 New Trends in Parallel Computing and Optimization,</span>
+class="cmr-12">, in PCO&#8217;11 New Trends in Parallel Computing and Optimization,</span>
      <span 
 class="cmr-12">IEEE  International  Symposium  on  Parallel  and  Distributed  Processing</span>
      <span 

docs/html/userhtmlse8.html

@@ -122,7 +122,7 @@ class="cmr-12">abide by its terms:</span>
 &#x00A0;<br />
 </div>
 <!--l. 87--><p class="nopar" > <span 
-class="cmr-12">AMG4PSBLAS is distributed together with (a small part) of the graph-matching</span>
+class="cmr-12">AMG4PSBLAS is distributed together with (a small part of) the graph-matching</span>
                                                                                
 
                                                                                
@@ -133,7 +133,7 @@ class="cmr-12">[</span><a
 href="userhtmlli5.html#XMatchBoxP"><span 
 class="cmr-12">9</span></a><span 
 class="cmr-12">]</span></span><span 
-class="cmr-12">. Per the license requirements, we reproduce the relative part</span>
+class="cmr-12">. Per the license requirements, we reproduce the relevant part</span>
 <span 
 class="cmr-12">here.</span>
                                                                                

docs/html/userhtmlse9.html

@@ -91,7 +91,7 @@ class="cmr-12">Trolling, insulting or derogatory comments, and personal or polit
 class="cmr-12">Public or private harassment</span>
      </li>
      <li class="itemize"><span 
-class="cmr-12">Publishing others private information, such as a physical or email address,</span>
+class="cmr-12">Publishing others&#8217; private information, such as a physical or email address,</span>
      <span 
 class="cmr-12">without their explicit permission</span>
      </li>
@@ -234,7 +234,7 @@ class="cmr-12">_of</span><span
 class="cmr-12">_conduct</span>
 <span 
 class="cmr-12">.html</span></a><span 
-class="cmr-12">. Community Impact Guidelines were inspired by Mozillas code of conduct</span>
+class="cmr-12">. Community Impact Guidelines were inspired by Mozilla&#8217;s code of conduct</span>
 <span 
 class="cmr-12">enforcement ladder. For answers to common questions about this code of conduct, see</span>
 <span 

docs/html/userhtmlsu3.html

@@ -4227,9 +4227,9 @@ class="cmtt-12">github</span><span
 class="cmtt-12">.</span><span 
 class="cmtt-12">com</span><span 
 class="cmtt-12">/</span><span 
-class="cmtt-12">sfilippone</span><span 
+class="cmtt-12">psctoolkit</span><span 
 class="cmtt-12">/</span><span 
-class="cmtt-12">amg4psblas</span><span 
+class="cmtt-12">psctoolkit</span><span 
 class="cmtt-12">/</span><span 
 class="cmtt-12">issues</span><span 
 class="cmtt-12">&#x003E;.</span>

docs/html/userhtmlsu4.html

@@ -36,8 +36,8 @@ class="cmr-12">If you find any bugs in our codes, please report them through our
 <span 
 class="cmr-12">on</span><br 
 class="newline" /> <a 
-href="https://github.com/psctoolkit/amg4psblas/issues" class="url" ><span 
-class="cmtt-12">https://github.com/psctoolkit/amg4psblas/issues</span></a><br 
+href="https://github.com/psctoolkit/psctoolkit/issues" class="url" ><span 
+class="cmtt-12">https://github.com/psctoolkit/psctoolkit/issues</span></a><br 
 class="newline" />
 <!--l. 195--><p class="indent" >   <span 
 class="cmr-12">To enable us to track the bug, please provide a log from the failing application, the</span>

docs/html/userhtmlsu5.html

@@ -29,40 +29,41 @@ class="cmr-12">3.5   </span></span> <a
  id="x13-120003.5"></a><span 
 class="cmr-12">Example and test programs</span></h4>
 <!--l. 200--><p class="noindent" ><span 
-class="cmr-12">The package contains the </span><span class="obeylines-h"><span class="verb"><span 
-class="cmtt-12">examples</span></span></span> <span 
-class="cmr-12">and </span><span class="obeylines-h"><span class="verb"><span 
-class="cmtt-12">tests</span></span></span> <span 
-class="cmr-12">directories; both of them are</span>
-<span 
-class="cmr-12">further divided into </span><span class="obeylines-h"><span class="verb"><span 
+class="cmr-12">The package contains a </span><span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">samples</span></span></span> <span 
+class="cmr-12">directory, divided in two subdirs </span><span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">simple</span></span></span> <span 
+class="cmr-12">and</span>
+<span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">advanced</span></span></span><span 
+class="cmr-12">; both of them are further divided into </span><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">fileread</span></span></span> <span 
 class="cmr-12">and </span><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">pdegen</span></span></span> <span 
-class="cmr-12">subdirectories. Their purpose is as</span>
+class="cmr-12">subdirectories.</span>
 <span 
-class="cmr-12">follows:</span>
+class="cmr-12">Their purpose is as follows:</span>
      <dl class="description"><dt class="description">
 <span 
-class="cmtt-12">examples</span> </dt><dd 
+class="cmtt-12">simple</span> </dt><dd 
 class="description"><span 
-class="cmr-12">contains a set of simple example programs with a predefined choice</span>
+class="cmr-12">contains a set of simple example programs with a predefined choice of</span>
      <span 
-class="cmr-12">of preconditioners, selectable via integer values. These are intended to get</span>
+class="cmr-12">preconditioners,  selectable  via  integer  values.  These  are  intended  to  get</span>
      <span 
 class="cmr-12">acquainted with the multilevel preconditioners available in AMG4PSBLAS.</span>
      </dd><dt class="description">
 <span 
-class="cmtt-12">tests</span> </dt><dd 
+class="cmtt-12">advanced</span> </dt><dd 
 class="description"><span 
-class="cmr-12">contains a set of more sophisticated examples that will allow the user, via</span>
+class="cmr-12">contains a set of more sophisticated examples that will allow the user,</span>
      <span 
-class="cmr-12">the input files in the </span><span class="obeylines-h"><span class="verb"><span 
+class="cmr-12">via the input files in the </span><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">runs</span></span></span> <span 
-class="cmr-12">subdirectories, to experiment with the full range</span>
+class="cmr-12">subdirectories, to experiment with the full</span>
      <span 
-class="cmr-12">of preconditioners implemented in the package.</span></dd></dl>
-<!--l. 213--><p class="noindent" ><span 
+class="cmr-12">range of preconditioners implemented in the package.</span></dd></dl>
+<!--l. 214--><p class="noindent" ><span 
 class="cmr-12">The </span><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">fileread</span></span></span> <span 
 class="cmr-12">directories contain sample programs that read sparse matrices from files,</span>

docs/html/userhtmlsu6.html

@@ -71,38 +71,39 @@ class="cmr-12">The part of the code dealing with reading and assembling the spar
 <span 
 class="cmr-12">right-hand side vector and the deallocation of the relevant data structures, performed</span>
 <span 
-class="cmr-12">through the PSBLAS routines for sparse matrix and vector management, is not</span>
+class="cmr-12">through the PSBLAS routines for sparse matrix and vector management,</span>
 <span 
-class="cmr-12">reported here for the sake of conciseness. The complete code can be found in the</span>
+class="cmr-12">is not reported here for the sake of conciseness. The complete code can be</span>
 <span 
-class="cmr-12">example program file </span><span class="obeylines-h"><span class="verb"><span 
+class="cmr-12">found in the example program file </span><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">amg_dexample_ml.f90</span></span></span><span 
-class="cmr-12">, in the directory </span><span class="obeylines-h"><span class="verb"><span 
-class="cmtt-12">examples/fileread</span></span></span> <span 
-class="cmr-12">of</span>
-<span 
-class="cmr-12">the AMG4PSBLAS implementation (see Section</span><span 
+class="cmr-12">, in the directory</span>
+<span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">samples/simple/file</span></span></span><span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">read</span></span></span> <span 
+class="cmr-12">of the AMG4PSBLAS implementation (see Section</span><span 
 class="cmr-12">&#x00A0;</span><a 
 href="userhtmlsu5.html#x13-120003.5"><span 
 class="cmr-12">3.5</span><!--tex4ht:ref: sec:ex_and_test --></a><span 
-class="cmr-12">). A sample test problem along</span>
+class="cmr-12">). A</span>
 <span 
-class="cmr-12">with the relevant input data is available in </span><span class="obeylines-h"><span class="verb"><span 
-class="cmtt-12">examples/fileread/runs</span></span></span><span 
-class="cmr-12">. For details on</span>
+class="cmr-12">sample test problem along with the relevant input data is available in</span>
+<span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">samples/simple/fileread/runs</span></span></span><span 
+class="cmr-12">. For details on the use of the PSBLAS routines, see</span>
 <span 
-class="cmr-12">the use of the PSBLAS routines, see the PSBLAS User&#8217;s Guide</span><span 
+class="cmr-12">the PSBLAS User&#8217;s Guide</span><span 
 class="cmr-12">&#x00A0;</span><span class="cite"><span 
 class="cmr-12">[</span><a 
 href="userhtmlli5.html#XPSBLASGUIDE"><span 
 class="cmr-12">20</span></a><span 
 class="cmr-12">]</span></span><span 
 class="cmr-12">.</span>
-<!--l. 138--><p class="indent" >   <span 
-class="cmr-12">The setup and application of the default multilevel preconditioner for the real single</span>
                                                                                
 
                                                                                
+<!--l. 138--><p class="indent" >   <span 
+class="cmr-12">The setup and application of the default multilevel preconditioner for the real single</span>
 <span 
 class="cmr-12">precision and the complex, single and double precision, versions are obtained</span>
 <span 
@@ -114,7 +115,8 @@ class="cmr-12">for</span>
 <span 
 class="cmr-12">details). If these versions are installed, the corresponding codes are available in</span>
 <span class="obeylines-h"><span class="verb"><span 
-class="cmtt-12">examples/fileread/</span></span></span><span 
+class="cmtt-12">samples/simple/file</span></span></span><span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">read</span></span></span><span 
 class="cmr-12">.</span>
                                                                                
 
@@ -278,7 +280,7 @@ class="cmr-12">For all the previous preconditioners, example programs where the
 class="cmr-12">and the right-hand side are generated by discretizing a PDE with Dirichlet</span>
 <span 
 class="cmr-12">boundary conditions are also available in the directory </span><span class="obeylines-h"><span class="verb"><span 
-class="cmtt-12">examples/pdegen</span></span></span><span 
+class="cmtt-12">samples/simple/pdegen</span></span></span><span 
 class="cmr-12">.</span>
                                                                                
 

docs/html/userhtmlsu7.html

@@ -145,7 +145,7 @@ class="cmr-12">GPU environment</span>
                                                                                
 <div class="center" 
 >
-<!--l. 552--><p class="noindent" >
+<!--l. 553--><p class="noindent" >
 <div class="minipage"><div class="verbatim" id="verbatim-12">
 &#x00A0;&#x00A0;call&#x00A0;desc_a%cnv(mold=igmold)
 &#x00A0;<br />&#x00A0;&#x00A0;call&#x00A0;a%cscnv(info,mold=agmold)
@@ -172,7 +172,7 @@ class="cmr-12">GPU environment</span>
 &#x00A0;<br />
 &#x00A0;<br />&#x00A0;
 </div>
-<!--l. 579--><p class="nopar" ></div></div>
+<!--l. 580--><p class="nopar" ></div></div>
 <br /> <div class="caption" 
 ><span class="id">Listing 7: </span><span  
 class="content">setup of a GPU-enabled test program part three.</span></div><!--tex4ht:label?: x16-15003r7 -->
@@ -180,25 +180,30 @@ class="content">setup of a GPU-enabled test program part three.</span></div><!--
 
                                                                                
    </div><hr class="endfloat" />
-<!--l. 587--><p class="indent" >   <span 
-class="cmr-12">It is very important to employ solvers that are suited to the GPU, i.e. solvers that</span>
+<!--l. 588--><p class="indent" >   <span 
+class="cmr-12">It is very important to employ smoothers and coarsest solvers that are suited to the</span>
 <span 
-class="cmr-12">do NOT employ triangular system solve kernels. Solvers that satisfy this constraint</span>
+class="cmr-12">GPU, i.e. methods that do NOT employ triangular system solve kernels. Methods that</span>
 <span 
-class="cmr-12">include:</span>
+class="cmr-12">satisfy this constraint include:</span>
      <ul class="itemize1">
      <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">JACOBI</span></span></span>
      </li>
      <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
+class="cmtt-12">BJAC</span></span></span> <span 
+class="cmr-12">with the following methods on the local blocks:</span>
+          <ul class="itemize2">
+          <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">INVK</span></span></span>
-     </li>
-     <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
+          </li>
+          <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">INVT</span></span></span>
-     </li>
-     <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
+          </li>
+          <li class="itemize"><span class="obeylines-h"><span class="verb"><span 
 class="cmtt-12">AINV</span></span></span></li></ul>
-<!--l. 596--><p class="noindent" ><span 
+     </li></ul>
+<!--l. 600--><p class="noindent" ><span 
 class="cmr-12">and their </span><span 
 class="cmmi-12">&#x2113;</span><sub><span 
 class="cmr-8">1</span></sub> <span 

docs/html/userhtmlsu9.html

@@ -324,8 +324,8 @@ class="td11"><!--l. 124--><p class="noindent" >      </td><td  style="white-spac
 class="td11"><!--l. 124--><p class="noindent" ><span 
 class="cmr-12">An auxiliary input argument that can be passed to the underlying</span>
 <span 
-class="cmr-12">objects.</span>                                                                             </td>
-    </tr></table></div>
+class="cmr-12">objects.</span>                                                                             </td> 
+</tr></table></div>
 <!--l. 129--><p class="noindent" >
 <!--l. 134--><p class="indent" >   <span 
 class="cmr-12">A variety of preconditioners can be obtained by setting the appropriate</span>

docs/src/building.tex

@@ -190,22 +190,23 @@ make install
 \subsection{Bug reporting}
 If you find any bugs in our codes, please report them through our
 issues page on \\[2mm]
-\url{https://github.com/psctoolkit/amg4psblas/issues}\\
+\url{https://github.com/psctoolkit/psctoolkit/issues}\\
 
 To enable us to track the bug, please provide a log from the failing
 application, the test conditions, and ideally a self-contained test
 program reproducing the issue.
 
 \subsection{Example and test programs\label{sec:ex_and_test}}
-The package contains the \verb|examples| and \verb|tests| directories;
+The package contains a \verb|samples| directory, divided in two
+subdirs  \verb|simple| and \verb|advanced|;
 both of them are further divided into \verb|fileread| and
 \verb|pdegen| subdirectories. Their purpose is as follows:
 \begin{description}
-\item[\tt examples] contains a set of simple example programs with a
+\item[\tt simple] contains a set of simple example programs with a
   predefined choice of preconditioners, selectable via integer
   values. These are intended to get  acquainted with the
   multilevel preconditioners available in AMG4PSBLAS.
-\item[\tt tests] contains a set of more sophisticated examples that
+\item[\tt advanced] contains a set of more sophisticated examples that
   will allow the user, via the input files in the \verb|runs|
   subdirectories, to experiment with the full range of preconditioners
   implemented in the package.

docs/src/configureout.txt

@@ -159,4 +159,4 @@ Some influential environment variables:
 Use these variables to override the choices made by `configure' or to help
 it to find libraries and programs with nonstandard names/locations.
 
-Report bugs to <https://github.com/sfilippone/amg4psblas/issues>.
+Report bugs to <https://github.com/psctoolkit/psctoolkit/issues>.

docs/src/gettingstarted.tex

@@ -129,9 +129,9 @@ relevant data structures, performed
 through the PSBLAS routines for sparse matrix and vector management, is not reported
 here for the sake of conciseness.
 The complete code can be found in the example program file \verb|amg_dexample_ml.f90|,
-in the directory \verb|examples/fileread| of the AMG4PSBLAS implementation (see
+in the directory \verb|samples/simple/file|\-\verb|read| of the AMG4PSBLAS implementation (see
 Section~\ref{sec:ex_and_test}). A sample test problem along with the relevant
-input data is available in \verb|examples/fileread/runs|.
+input data is available in \verb|samples/simple/fileread/runs|.
 For details on the use of the PSBLAS routines, see the PSBLAS User's
 Guide~\cite{PSBLASGUIDE}.
 
@@ -139,7 +139,7 @@ The setup and application of the default multilevel preconditioner
 for the real single precision and the complex, single and double
 precision, versions are obtained with straightforward modifications of the previous
 example (see Section~\ref{sec:userinterface} for details). If these versions are installed,
-the corresponding codes are available in \verb|examples/fileread/|.
+the corresponding codes are available in \verb|samples/simple/file|\-\verb|read|.
 
 \begin{listing}[tbp]
 \begin{center}
@@ -300,7 +300,7 @@ The corresponding example program is available in the file
 
 For all the previous preconditioners, example programs where the sparse matrix and
 the right-hand side are generated by discretizing a PDE with Dirichlet
-boundary conditions are also available in the directory \verb|examples/pdegen|.
+boundary conditions are also available in the directory \verb|samples/simple/pdegen|.
 \vspace{-1em}\begin{listing}[tbh]
 \ifpdf%
 \begin{minted}[breaklines=true,bgcolor=bg,fontsize=\small]{fortran}
@@ -535,7 +535,8 @@ Krylov method. At the end of the code, we close the GPU environment
   call prec%allocate_wrk(info)
   t1 = psb_wtime()
   call psb_krylov(s_choice%kmethd,a,prec,b,x,s_choice%eps,&
-       & desc_a,info,itmax=s_choice%itmax,iter=iter,err=err,itrace=s_choice%itrace,&
+       & desc_a,info,itmax=s_choice%itmax,iter=iter,err=err,&
+       & itrace=s_choice%itrace,&
        & istop=s_choice%istopc,irst=s_choice%irst)
   call prec%deallocate_wrk(info)
   call psb_barrier(ctxt)
@@ -584,15 +585,18 @@ Krylov method. At the end of the code, we close the GPU environment
 \caption{setup of a GPU-enabled test program part three.\label{fig:gpu-ex3}}
 \end{listing}
 
-It is very important to employ solvers that are suited
-to the GPU, i.e.  solvers that do NOT  employ triangular
-system solve kernels.  Solvers that satisfy this constraint include:
+It is very important to employ smoothers and coarsest solvers that are suited
+to the GPU, i.e.  methods that do NOT employ triangular
+system solve kernels. Methods that satisfy this constraint include:
 \begin{itemize}
 \item \verb|JACOBI|
+\item \verb|BJAC| with the following methods on the local blocks:
+\begin{itemize}
 \item \verb|INVK|
 \item \verb|INVT|
 \item \verb|AINV|
 \end{itemize}
+\end{itemize}
 and their $\ell_1$ variants. 
 
 %%% Local Variables:

docs/src/license.tex

@@ -87,9 +87,9 @@ terms: {\small
 \end{verbatim}
    }
 \pagebreak
-AMG4PSBLAS is distributed together with (a small part) of the graph-matching
+AMG4PSBLAS is distributed together with (a small part of) the graph-matching
 library MatchBox-P~\cite{MatchBoxP}. Per the license requirements, we reproduce
-the relative part here.
+the relevant part here.
 {\small
 \begin{verbatim}
 // ***********************************************************************

docs/src/userguide.tex

@@ -154,7 +154,7 @@ Preconditioners Package based on PSBLAS}
 \flushright
 \large Software version: 1.0\\
 %\todaym
-\large April 12, 2021
+\large May 11th, 2021
 \end{minipage}}
 %\addtolength{\textwidth}{\centeroffset}
 \vspace{\stretch{2}}

docs/src/userhtml.tex

@@ -114,7 +114,7 @@
 %\today
 Software version: 1.0\\
 %\today
- April 12, 2021
+ May 11th, 2021
 \clearpage
 \ \\
 \thispagestyle{empty}

examples/gpu/amg_dexample_gpu.f90

@@ -39,23 +39,18 @@
 !
 ! This sample program solves a linear system obtained by discretizing a
 ! PDE with Dirichlet BCs. The solver is CG, coupled with one of the
-! following multi-level preconditioner, as explained in Section 4.1 of
+! following multi-level preconditioner, as explained in Section 4.2 of
 ! the AMG4PSBLAS User's and Reference Guide:
 !
-! - choice = 1, the default multi-level preconditioner solver, i.e., 
-! V-cycle with decoupled smoothed aggregation, 1 hybrid forward/backward
-! GS sweep as pre/post-smoother and UMFPACK as coarsest-level
-! solver (Sec. 4.1, Listing 1)
+! - choice = 1, a V-cycle with decoupled smoothed aggregation, 4 Jacobi
+! sweeps as pre/post-smoother and 8 Jacobi sweeps as coarsest-level
+! solver with replicated coarsest matrix
 !
-! - choice = 2, a V-cycle preconditioner with 1 block-Jacobi sweep
-! (with ILU(0) on the blocks) as pre- and post-smoother, and 8 block-Jacobi
-! sweeps (with ILU(0) on the blocks) as coarsest-level solver (Sec. 4.1, Listing 2)
-!
-! - choice = 3,  W-cycle preconditioner based on the coupled aggregation relying
-! on matching, with maximum size of aggregates equal to 8 and smoothed prolongators,
-! 2 hybrid forward/backward GS sweeps as pre/post-smoother, a distributed coarsest
-! matrix, and preconditioned Flexible Conjugate Gradient as coarsest-level solver
-! (Sec. 4.1, Listing 3)
+! - choice = 2, a W-cycle based on the coupled aggregation relying on matching,
+! with maximum size of aggregates equal to 8 and smoothed prolongators,
+! 2 sweeps of Block-Jacobi ipre/post-smoother using approximate inverse INVK and
+! 4 sweeps of Block-Jacobi with INVK as coarsest-level solver on distributed
+! coarsest matrix 
 !
 ! The matrix and the rhs are read from files (if an rhs is not available, the
 ! unit rhs is set).
@@ -183,8 +178,9 @@ program amg_dexample_gpu
 
  case(1)
 
-    ! initialize a V-cycle preconditioner with 4 Jacobi sweep 
-    ! and 8 Jacobi sweeps as coarsest-level solver
+    ! initialize a V-cycle preconditioner, relying on decoupled smoothed aggregation 
+    ! with 4 Jacobi sweeps as pre/post-smoother
+    ! and 8 Jacobi sweeps as coarsest-level solver on replicated coarsest matrix
 
     call P%init(ctxt,'ML',info)
     call P%set('SMOOTHER_TYPE','JACOBI',info)
@@ -195,19 +191,22 @@ program amg_dexample_gpu
 
   case(2)
 
-   ! initialize a V-cycle preconditioner based on the coupled aggregation relying on matching,
+   ! initialize a W-cycle preconditioner based on the coupled aggregation relying on matching,
    ! with maximum size of aggregates equal to 8 and smoothed prolongators,
-   ! Block-Jacobi smoother using approximate inverse INVK and
-   ! and 4 sweeps of INVK on he coarsest level 
+   ! 2 sweeps of Block-Jacobi pre/post-smoother using approximate inverse INVK and
+   ! 4 sweeps of Block-Jacobi with INVK on the coarsest level distributed matrix
 
     call P%init(ctxt,'ML',info)
     call P%set('PAR_AGGR_ALG','COUPLED',info)
     call P%set('AGGR_TYPE','MATCHBOXP',info)
     call P%set('AGGR_SIZE',8,info)
     call P%set('ML_CYCLE','WCYCLE',info)
+    call P%set('SMOOTHER_TYPE','BJAC',info)
     call P%set('SMOOTHER_SWEEPS',2,info)
     call P%set('SUB_SOLVE','INVK',info)  
-    call P%set('COARSE_SOLVE','INVK',info)
+    call P%set('COARSE_SOLVE','BJAC',info)
+    call P%set('COARSE_SUBSOLVE','INVK',info)
+    call P%set('COARSE_SWEEPS',4,info)
     call P%set('COARSE_MAT','DIST',info)
     kmethod = 'CG'
 

samples/advanced/fileread/Makefile

@@ -1,4 +1,4 @@
-AMGDIR=../..
+AMGDIR=../../..
 AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules

samples/advanced/fileread/amg_cf_sample.f90

@@ -523,7 +523,7 @@ program amg_cf_sample
   call psb_sum(ctxt,amatsize)
   call psb_sum(ctxt,descsize)
   call psb_sum(ctxt,precsize)
-  call prec%descr(iout=psb_out_unit)
+  call prec%descr(info,iout=psb_out_unit)
   if (iam == psb_root_) then 
     write(psb_out_unit,'("Matrix: ",a)')mtrx_file
     write(psb_out_unit,'("Computed solution on ",i8," processors")')np

samples/advanced/fileread/amg_df_sample.f90

@@ -523,7 +523,7 @@ program amg_df_sample
   call psb_sum(ctxt,amatsize)
   call psb_sum(ctxt,descsize)
   call psb_sum(ctxt,precsize)
-  call prec%descr(iout=psb_out_unit)
+  call prec%descr(info,iout=psb_out_unit)
   if (iam == psb_root_) then 
     write(psb_out_unit,'("Matrix: ",a)')mtrx_file
     write(psb_out_unit,'("Computed solution on ",i8," processors")')np

samples/advanced/fileread/amg_sf_sample.f90

@@ -523,7 +523,7 @@ program amg_sf_sample
   call psb_sum(ctxt,amatsize)
   call psb_sum(ctxt,descsize)
   call psb_sum(ctxt,precsize)
-  call prec%descr(iout=psb_out_unit)
+  call prec%descr(info,iout=psb_out_unit)
   if (iam == psb_root_) then 
     write(psb_out_unit,'("Matrix: ",a)')mtrx_file
     write(psb_out_unit,'("Computed solution on ",i8," processors")')np

samples/advanced/fileread/amg_zf_sample.f90

@@ -523,7 +523,7 @@ program amg_zf_sample
   call psb_sum(ctxt,amatsize)
   call psb_sum(ctxt,descsize)
   call psb_sum(ctxt,precsize)
-  call prec%descr(iout=psb_out_unit)
+  call prec%descr(info,iout=psb_out_unit)
   if (iam == psb_root_) then 
     write(psb_out_unit,'("Matrix: ",a)')mtrx_file
     write(psb_out_unit,'("Computed solution on ",i8," processors")')np

samples/advanced/fileread/data_input.f90

@@ -1,14 +1,14 @@
 !   
 !   
-!                             MLD2P4  version 2.2
-!    MultiLevel Domain Decomposition Parallel Preconditioners Package
-!               based on PSBLAS (Parallel Sparse BLAS version 3.5)
+!                             AMG4PSBLAS version 1.0
+!    Algebraic Multigrid Package
+!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
 !    
-!    (C) Copyright 2008-2018 
+!    (C) Copyright 2021 
 !  
 !        Salvatore Filippone  
 !        Pasqua D'Ambra   
-!        Daniela di Serafino   
+!        Fabio Durastante        
 !   
 !    Redistribution and use in source and binary forms, with or without
 !    modification, are permitted provided that the following conditions
@@ -18,14 +18,14 @@
 !      2. Redistributions in binary form must reproduce the above copyright
 !         notice, this list of conditions, and the following disclaimer in the
 !         documentation and/or other materials provided with the distribution.
-!      3. The name of the MLD2P4 group or the names of its contributors may
+!      3. The name of the AMG4PSBLAS group or the names of its contributors may
 !         not be used to endorse or promote products derived from this
 !         software without specific written permission.
 !   
 !    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
 !    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 !    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 !    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS