Adjusted CBIND test program.

Fixes for krylov->linsolve. Test program to be completed.
Switch from KRYLOV to LINSOLVE
2026-10-07 07:04:59 +00:00 · 2024-11-11 13:06:48 +01:00 · 2024-11-11 09:11:07 +01:00 · 2024-11-10 14:00:51 +01:00 · 2024-11-08 18:25:13 +01:00 · 2024-10-16 10:56:53 +02:00
104 changed files with 655 additions and 417 deletions
@@ -49,7 +49,7 @@ PSBLAS_INCLUDES=@PSBLAS_INCLUDES@
 PSBLAS_LIBS=@PSBLAS_LIBS@	
 PSBBASEMODNAME=psb_base_mod
 PSBPRECMODNAME=psb_prec_mod
-PSBMETHDMODNAME=psb_krylov_mod
+PSBMETHDMODNAME=psb_linsolve_mod
 PSBUTILMODNAME=psb_util_mod


@@ -3,9 +3,9 @@ include Make.inc

 all:  objs lib

-objs: amgp cbnd
+objs: libdir amgp cbnd

-lib: libdir objs
+lib:  objs
 	cd amgprec && $(MAKE) lib
 	cd cbind && $(MAKE)  lib

@@ -46,6 +46,7 @@ cleanlib:

 veryclean: cleanlib
 	(cd amgprec && $(MAKE) veryclean)
+	(cd cbind && $(MAKE) veryclean)
 	(cd samples/simple/fileread && $(MAKE) clean)
 	(cd samples/simple/pdegen && $(MAKE) clean)
 	(cd samples/advanced/fileread && $(MAKE) clean)
@@ -56,3 +57,4 @@ check: all

 clean:
 	(cd amgprec && $(MAKE) clean)
+	(cd cbind && $(MAKE) clean)
@@ -323,10 +323,10 @@ module amg_base_prec_type
  !
  ! Legal values for entry: amg_poly_variant_
  !
-  integer(psb_ipk_), parameter :: amg_poly_lottes_      = 0
-  integer(psb_ipk_), parameter :: amg_poly_lottes_beta_ = 1
-  integer(psb_ipk_), parameter :: amg_poly_new_         = 2
-  integer(psb_ipk_), parameter :: amg_poly_dbg_         = 8
+  integer(psb_ipk_), parameter :: amg_cheb_4_     = 0
+  integer(psb_ipk_), parameter :: amg_cheb_4_opt_ = 1
+  integer(psb_ipk_), parameter :: amg_cheb_1_opt_ = 2
+  integer(psb_ipk_), parameter :: amg_poly_dbg_   = 8
  
  integer(psb_ipk_), parameter :: amg_poly_rho_est_power_ = 0

@@ -570,12 +570,12 @@ contains
      val = amg_as_
    case('POLY')
      val = amg_poly_
-    case('POLY_LOTTES')
-      val = amg_poly_lottes_
-    case('POLY_LOTTES_BETA')
-      val = amg_poly_lottes_beta_
-    case('POLY_NEW')
-      val = amg_poly_new_
+    case('CHEB_4')
+      val = amg_cheb_4_
+    case('CHEB_4_OPT')
+      val = amg_cheb_4_opt_
+    case('CHEB_1_OPT')
+      val = amg_cheb_1_opt_
    case('POLY_DBG')
      val = amg_poly_dbg_
    case('POLY_RHO_EST_POWER')
@@ -322,7 +322,7 @@ contains
    !
    sm%pdegree      = 1
    sm%rho_ba       = -done
-    sm%variant      = amg_poly_lottes_
+    sm%variant      = amg_cheb_4_
    sm%rho_estimate = amg_poly_rho_est_power_
    sm%rho_estimate_iterations = 20
    if (allocated(sm%sv)) then
@@ -272,7 +272,9 @@ contains
          write(0,*) 'Impossible: mate(k) > nc'
          cycle
        else
+
          if (ilaggr(k) == ilaggr_neginit) then
+
            wk   = w(k)
            widx = w(idx)
            wmax = max(abs(wk),abs(widx))
@@ -322,7 +322,7 @@ contains
    !
    sm%pdegree      = 1
    sm%rho_ba       = -sone
-    sm%variant      = amg_poly_lottes_
+    sm%variant      = amg_cheb_4_
    sm%rho_estimate = amg_poly_rho_est_power_
    sm%rho_estimate_iterations = 20
    if (allocated(sm%sv)) then
@@ -82,6 +82,8 @@ const int BundleTag = 9;  // Predefined tag

 static vector<MilanLongInt> DEFAULT_VECTOR;

+#if !defined(SERIAL_MPI)
+
 // MPI type map
 template <typename T>
 MPI_Datatype TypeMap();
@@ -93,6 +95,7 @@ template <>
 inline MPI_Datatype TypeMap<double>() { return MPI_DOUBLE; }
 template <>
 inline MPI_Datatype TypeMap<float>() { return MPI_FLOAT; }
+#endif

 #ifdef __cplusplus
 extern "C"
@@ -300,7 +300,7 @@ subroutine amg_caggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if
@@ -230,7 +230,7 @@ subroutine amg_caggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -246,7 +246,7 @@ subroutine amg_d_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -300,7 +300,7 @@ subroutine amg_daggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if
@@ -230,7 +230,7 @@ subroutine amg_daggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -246,7 +246,7 @@ subroutine amg_s_parmatch_smth_bld(ag,a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -300,7 +300,7 @@ subroutine amg_saggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if
@@ -230,7 +230,7 @@ subroutine amg_saggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -300,7 +300,7 @@ subroutine amg_zaggrmat_minnrg_bld(a,desc_a,ilaggr,nlaggr,parms,&
 !!$        endif
 !!$      enddo
 !!$      if (jd == -1) then 
-!!$        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+!!$                write(0,*) name,': Warning: there is no diagonal element', i
 !!$      else
 !!$        acsrf%val(jd)=acsrf%val(jd)-tmp
 !!$      end if
@@ -230,7 +230,7 @@ subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&

      enddo
      if (jd == -1) then 
-        write(0,*) 'Wrong input: we need the diagonal!!!!', i
+        write(0,*) name,': Warning: there is no diagonal element', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
@@ -1,5 +1,6 @@
 #include "MatchBoxPC.h"
 // TODO comment
+#if !defined(SERIAL_MPI)

 void clean(MilanLongInt NLVer,
           MilanInt myRank,
@@ -88,3 +89,4 @@ void clean(MilanLongInt NLVer,
        }
    }
 }
+#endif
@@ -8,6 +8,7 @@
 * @param edgeLocWeight
 * @return
 */
+
 MilanLongInt firstComputeCandidateMateD(MilanLongInt adj1,
                                       MilanLongInt adj2,
                                       MilanLongInt *verLocInd,
@@ -134,3 +135,4 @@ MilanLongInt computeCandidateMateS(MilanLongInt adj1,

    return w;
 }
+
@@ -1,4 +1,5 @@
 #include "MatchBoxPC.h"
+
 void PARALLEL_COMPUTE_CANDIDATE_MATE_BD(MilanLongInt NLVer,
                                              MilanLongInt *verLocPtr,
                                              MilanLongInt *verLocInd,
@@ -52,3 +53,4 @@ void PARALLEL_COMPUTE_CANDIDATE_MATE_BS(MilanLongInt NLVer,
        }
    }
 }
+
@@ -366,3 +366,4 @@ void PARALLEL_PROCESS_EXPOSED_VERTEX_BS(MilanLongInt NLVer,
 	
    } // End of parallel region
 }
+
@@ -582,3 +582,4 @@ void processMatchedVerticesS(
 #endif
    } // End of parallel region
 }
+
@@ -588,3 +588,4 @@ void processMatchedVerticesAndSendMessagesS(
    cout << myRank<<" Done sending messages"<<endl;
 #endif
 }
+
@@ -1,5 +1,6 @@
 #include "MatchBoxPC.h"
 //#define DEBUG_HANG_
+#if !defined(SERIAL_MPI)

 void processMessagesD(
    MilanLongInt NLVer,
@@ -629,3 +630,4 @@ void processMessagesS(

    return;
 }
+#endif
@@ -30,3 +30,4 @@ void queuesTransfer(vector<MilanLongInt> &U,
    privateQOwner.clear();

 }
+
@@ -451,11 +451,16 @@ subroutine amg_c_hierarchy_bld(a,desc_a,prec,info)
    if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
      prec%precv(1)%base_desc => prec%precv(1)%desc_ac
    end if
-    do i=2, iszv 
+    do i=2, iszv
      prec%precv(i)%base_a       => prec%precv(i)%ac
      prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
    end do
  end if

@@ -451,11 +451,16 @@ subroutine amg_d_hierarchy_bld(a,desc_a,prec,info)
    if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
      prec%precv(1)%base_desc => prec%precv(1)%desc_ac
    end if
-    do i=2, iszv 
+    do i=2, iszv
      prec%precv(i)%base_a       => prec%precv(i)%ac
      prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
    end do
  end if

@@ -451,11 +451,16 @@ subroutine amg_s_hierarchy_bld(a,desc_a,prec,info)
    if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
      prec%precv(1)%base_desc => prec%precv(1)%desc_ac
    end if
-    do i=2, iszv 
+    do i=2, iszv
      prec%precv(i)%base_a       => prec%precv(i)%ac
      prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
    end do
  end if

@@ -451,11 +451,16 @@ subroutine amg_z_hierarchy_bld(a,desc_a,prec,info)
    if (.not.associated(prec%precv(1)%base_desc,desc_a)) then
      prec%precv(1)%base_desc => prec%precv(1)%desc_ac
    end if
-    do i=2, iszv 
+    do i=2, iszv
      prec%precv(i)%base_a       => prec%precv(i)%ac
      prec%precv(i)%base_desc    => prec%precv(i)%desc_ac
-      prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
-      prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
+      ! This is needed when the linmap object has been built
+      ! reusing the base_desc descriptor through a pointer.
+      ! With PSBLAS 4 we will have a better solution
+      if (associated(prec%precv(i)%linmap%p_desc_U)) &
+           & prec%precv(i)%linmap%p_desc_U => prec%precv(i-1)%base_desc
+      if (associated(prec%precv(i)%linmap%p_desc_V))&
+           & prec%precv(i)%linmap%p_desc_V => prec%precv(i)%base_desc
    end do
  end if

@@ -241,6 +241,8 @@ subroutine amg_c_base_onelev_csetc(lv,what,val,info,pos,idx)
      if (info == 0) deallocate(lv%aggr,stat=info)
      if (info /= 0) then
        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
        return
      end if
    end if
@@ -250,8 +252,16 @@ subroutine amg_c_base_onelev_csetc(lv,what,val,info,pos,idx)
      allocate(amg_c_dec_aggregator_type :: lv%aggr, stat=info)
    case('SYMDEC')
      allocate(amg_c_symdec_aggregator_type :: lv%aggr, stat=info)
-    case default
+#if !defined(SERIAL_MPI)
+#endif
+  case default
      info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
    end select
    if (info == psb_success_) call lv%aggr%default()

@@ -127,8 +127,7 @@ subroutine amg_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
        ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
      end if
    end if
  else
@@ -151,8 +150,7 @@ subroutine amg_c_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
      if (tprol_) then 
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
      end if
    end if
  end if
@@ -42,7 +42,9 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
  use amg_d_base_aggregator_mod
  use amg_d_dec_aggregator_mod
  use amg_d_symdec_aggregator_mod
+#if !defined(SERIAL_MPI)
  use amg_d_parmatch_aggregator_mod
+#endif
  use amg_d_poly_smoother
  use amg_d_jac_smoother
  use amg_d_as_smoother
@@ -267,6 +269,8 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
      if (info == 0) deallocate(lv%aggr,stat=info)
      if (info /= 0) then
        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
        return
      end if
    end if
@@ -276,10 +280,18 @@ subroutine amg_d_base_onelev_csetc(lv,what,val,info,pos,idx)
      allocate(amg_d_dec_aggregator_type :: lv%aggr, stat=info)
    case('SYMDEC')
      allocate(amg_d_symdec_aggregator_type :: lv%aggr, stat=info)
+#if !defined(SERIAL_MPI)
    case('COUP','COUPLED')
      allocate(amg_d_parmatch_aggregator_type :: lv%aggr, stat=info)
-    case default
+#endif
+  case default
      info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
    end select
    if (info == psb_success_) call lv%aggr%default()

@@ -127,8 +127,7 @@ subroutine amg_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
        ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
      end if
    end if
  else
@@ -151,8 +150,7 @@ subroutine amg_d_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
      if (tprol_) then 
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
      end if
    end if
  end if
@@ -42,7 +42,9 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
  use amg_s_base_aggregator_mod
  use amg_s_dec_aggregator_mod
  use amg_s_symdec_aggregator_mod
+#if !defined(SERIAL_MPI)
  use amg_s_parmatch_aggregator_mod
+#endif
  use amg_s_poly_smoother
  use amg_s_jac_smoother
  use amg_s_as_smoother
@@ -247,6 +249,8 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
      if (info == 0) deallocate(lv%aggr,stat=info)
      if (info /= 0) then
        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
        return
      end if
    end if
@@ -256,10 +260,18 @@ subroutine amg_s_base_onelev_csetc(lv,what,val,info,pos,idx)
      allocate(amg_s_dec_aggregator_type :: lv%aggr, stat=info)
    case('SYMDEC')
      allocate(amg_s_symdec_aggregator_type :: lv%aggr, stat=info)
+#if !defined(SERIAL_MPI)
    case('COUP','COUPLED')
      allocate(amg_s_parmatch_aggregator_type :: lv%aggr, stat=info)
-    case default
+#endif
+  case default
      info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
    end select
    if (info == psb_success_) call lv%aggr%default()

@@ -127,8 +127,7 @@ subroutine amg_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
        ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
      end if
    end if
  else
@@ -151,8 +150,7 @@ subroutine amg_s_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
      if (tprol_) then 
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
      end if
    end if
  end if
@@ -261,6 +261,8 @@ subroutine amg_z_base_onelev_csetc(lv,what,val,info,pos,idx)
      if (info == 0) deallocate(lv%aggr,stat=info)
      if (info /= 0) then
        info = psb_err_internal_error_
+        call psb_errpush(info,name,a_err='aggregator deallocation?')
+        goto 9999
        return
      end if
    end if
@@ -270,8 +272,16 @@ subroutine amg_z_base_onelev_csetc(lv,what,val,info,pos,idx)
      allocate(amg_z_dec_aggregator_type :: lv%aggr, stat=info)
    case('SYMDEC')
      allocate(amg_z_symdec_aggregator_type :: lv%aggr, stat=info)
-    case default
+#if !defined(SERIAL_MPI)
+#endif
+  case default
      info =  psb_err_internal_error_
+#if !defined(SERIAL_MPI)
+      call psb_errpush(info,name,a_err='Unsupported PAR_AGGR_ALG')
+#else
+      call psb_errpush(info,name,a_err='PAR_AGGR_ALG unsupported (SERIAL_MPI on)')
+#endif
+      goto 9999
    end select
    if (info == psb_success_) call lv%aggr%default()

@@ -127,8 +127,7 @@ subroutine amg_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
        ivr = lv%linmap%p_desc_U%get_global_indices(owned=.false.)
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head,ivr=ivr)
+        call lv%tprol%print(fname,head=head,ivr=ivr)
      end if
    end if
  else
@@ -151,8 +150,7 @@ subroutine amg_z_base_onelev_dump(lv,level,info,prefix,head,ac,rp,&
      if (tprol_) then 
        write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_tprol.mtx'
        !
-        !  This is not implemented yet.
-        !call lv%tprol%print(fname,head=head)
+        call lv%tprol%print(fname,head=head)
      end if
    end if
  end if
@@ -40,7 +40,7 @@ subroutine amg_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_c_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_c_jac_smoother, amg_protect_name => amg_c_jac_smoother_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_d_jac_smoother, amg_protect_name => amg_d_jac_smoother_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_d_poly_smoother, amg_protect_name => amg_d_poly_smoother_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -140,7 +140,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
    call tz%zero()

    select case(sm%variant)
-    case(amg_poly_lottes_)
+    case(amg_cheb_4_)
      if (do_timings) call psb_tic(poly_1)
      block
        real(psb_dpk_)      :: cz, cr
@@ -155,7 +155,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
          cz = (2*i*done-3)/(2*i*done+done)
          cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
          if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
+          call psb_upd_xyz(cr,cz,done,done,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
          if (do_timings) call psb_toc(poly_vect)
          if (do_timings) call psb_tic(poly_mv)
          call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
@@ -167,12 +167,12 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
        cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
        cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
        if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,done,done,ty,tz,tx,desc_data,info)
        if (do_timings) call psb_toc(poly_vect)
      end block
      if (do_timings) call psb_toc(poly_1)

-    case(amg_poly_lottes_beta_)
+    case(amg_cheb_4_opt_)
      if (do_timings) call psb_tic(poly_2)
      block
        real(psb_dpk_)      :: cz, cr
@@ -195,7 +195,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
          cz = (2*i*done-3)/(2*i*done+done)
          cr = (8*i*done-4)/((2*i*done+done)*sm%rho_ba)
          if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
+          call psb_upd_xyz(cr,cz,sm%poly_beta(i),done,ty,tz,tx,desc_data,info)
          if (do_timings) call psb_toc(poly_vect)
          if (do_timings) call psb_tic(poly_mv)
          call psb_spmm(-done,sm%pa,tz,done,r,desc_data,info,work=aux,trans=trans_)
@@ -205,11 +205,11 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
        cz = (2*sm%pdegree*done-3)/(2*sm%pdegree*done+done)
        cr = (8*sm%pdegree*done-4)/((2*sm%pdegree*done+done)*sm%rho_ba)
        if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,sm%poly_beta(sm%pdegree),done,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,sm%poly_beta(sm%pdegree),done,ty,tz,tx,desc_data,info)
        if (do_timings) call psb_toc(poly_vect)
      end block
      if (do_timings) call psb_toc(poly_2)
-    case(amg_poly_new_)
+    case(amg_cheb_1_opt_)
      if (do_timings) call psb_tic(poly_3)
      block
        real(psb_dpk_)      :: sigma, theta, delta, rho_old, rho
@@ -226,7 +226,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
        if (do_timings) call psb_toc(poly_sv)
        call psb_geaxpby((done/sm%rho_ba),ty,dzero,r,desc_data,info)
        if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
+        call psb_upd_xyz((done/theta),dzero,done,done,r,tz,tx,desc_data,info)
        if (do_timings) call psb_toc(poly_vect)

        ! tz == d
@@ -244,7 +244,7 @@ subroutine amg_d_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
          ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
          rho = done/(2*sigma - rho_old)
          if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
+          call psb_upd_xyz((2*rho/delta),(rho*rho_old),done,done,r,tz,tx,desc_data,info)
          if (do_timings) call psb_toc(poly_vect)
          rho_old = rho
        end do
@@ -75,9 +75,9 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
  nrow_a = a%get_nrows()
  nztota = a%get_nzeros()
  select case(sm%variant)
-  case(amg_poly_lottes_)
+  case(amg_cheb_4_)
    ! do nothing
-  case(amg_poly_lottes_beta_)
+  case(amg_cheb_4_opt_)
    if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
      call psb_realloc(sm%pdegree,sm%poly_beta,info)
      sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
@@ -87,7 +87,7 @@ subroutine amg_d_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
           & a_err='invalid sm%degree for poly_beta')
      goto 9999
    end if
-  case(amg_poly_new_)
+  case(amg_cheb_1_opt_)

    if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
      !Ok
@@ -58,11 +58,11 @@ subroutine amg_d_poly_smoother_cseti(sm,what,val,info,idx)
    sm%pdegree = val
  case('POLY_VARIANT')
    select case(val)
-    case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_)
+    case(amg_cheb_4_,amg_cheb_4_opt_,amg_cheb_1_opt_)
      sm%variant = val
    case default
-      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val
-      sm%variant = amg_poly_lottes_
+      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_cheb_4_',val
+      sm%variant = amg_cheb_4_
    end select
  case('POLY_RHO_ESTIMATE')
    select case(val)
@@ -77,17 +77,17 @@ subroutine amg_d_poly_smoother_descr(sm,info,iout,coarse,prefix)

  write(iout_,*) trim(prefix_), '  Polynomial smoother   '
  select case(sm%variant)
-  case(amg_poly_lottes_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES'
+  case(amg_cheb_4_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4'
    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
-  case(amg_poly_lottes_beta_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES_BETA'
+  case(amg_cheb_4_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4_OPT'
    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
    if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), '   Coefficients:        ',sm%poly_beta(1:sm%pdegree)
-  case(amg_poly_new_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_NEW'
+  case(amg_cheb_1_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_1_OPT'
    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
    write(iout_,*) trim(prefix_), '    Coefficient:        ',sm%cf_a
@@ -40,7 +40,7 @@ subroutine amg_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_s_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_s_jac_smoother, amg_protect_name => amg_s_jac_smoother_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_s_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_s_poly_smoother, amg_protect_name => amg_s_poly_smoother_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -140,7 +140,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
    call tz%zero()

    select case(sm%variant)
-    case(amg_poly_lottes_)
+    case(amg_cheb_4_)
      if (do_timings) call psb_tic(poly_1)
      block
        real(psb_spk_)      :: cz, cr
@@ -155,7 +155,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
          cz = (2*i*sone-3)/(2*i*sone+sone)
          cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
          if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
+          call psb_upd_xyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info) ! zk = cz * zk-1 + cr * rk-1 
          if (do_timings) call psb_toc(poly_vect)
          if (do_timings) call psb_tic(poly_mv)
          call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
@@ -167,12 +167,12 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
        cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
        cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
        if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,sone,sone,ty,tz,tx,desc_data,info)
        if (do_timings) call psb_toc(poly_vect)
      end block
      if (do_timings) call psb_toc(poly_1)

-    case(amg_poly_lottes_beta_)
+    case(amg_cheb_4_opt_)
      if (do_timings) call psb_tic(poly_2)
      block
        real(psb_spk_)      :: cz, cr
@@ -195,7 +195,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
          cz = (2*i*sone-3)/(2*i*sone+sone)
          cr = (8*i*sone-4)/((2*i*sone+sone)*sm%rho_ba)
          if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
+          call psb_upd_xyz(cr,cz,sm%poly_beta(i),sone,ty,tz,tx,desc_data,info)
          if (do_timings) call psb_toc(poly_vect)
          if (do_timings) call psb_tic(poly_mv)
          call psb_spmm(-sone,sm%pa,tz,sone,r,desc_data,info,work=aux,trans=trans_)
@@ -205,11 +205,11 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
        cz = (2*sm%pdegree*sone-3)/(2*sm%pdegree*sone+sone)
        cr = (8*sm%pdegree*sone-4)/((2*sm%pdegree*sone+sone)*sm%rho_ba)
        if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz(cr,cz,sm%poly_beta(sm%pdegree),sone,ty,tz,tx,desc_data,info)
+        call psb_upd_xyz(cr,cz,sm%poly_beta(sm%pdegree),sone,ty,tz,tx,desc_data,info)
        if (do_timings) call psb_toc(poly_vect)
      end block
      if (do_timings) call psb_toc(poly_2)
-    case(amg_poly_new_)
+    case(amg_cheb_1_opt_)
      if (do_timings) call psb_tic(poly_3)
      block
        real(psb_spk_)      :: sigma, theta, delta, rho_old, rho
@@ -226,7 +226,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
        if (do_timings) call psb_toc(poly_sv)
        call psb_geaxpby((sone/sm%rho_ba),ty,szero,r,desc_data,info)
        if (do_timings) call psb_tic(poly_vect)
-        call psb_abgdxyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
+        call psb_upd_xyz((sone/theta),szero,sone,sone,r,tz,tx,desc_data,info)
        if (do_timings) call psb_toc(poly_vect)

        ! tz == d
@@ -244,7 +244,7 @@ subroutine amg_s_poly_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
          ! d_{k+1} = (rho rho_old) d_k + 2(rho/delta) r_{k+1}
          rho = sone/(2*sigma - rho_old)
          if (do_timings) call psb_tic(poly_vect)
-          call psb_abgdxyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
+          call psb_upd_xyz((2*rho/delta),(rho*rho_old),sone,sone,r,tz,tx,desc_data,info)
          if (do_timings) call psb_toc(poly_vect)
          rho_old = rho
        end do
@@ -75,9 +75,9 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
  nrow_a = a%get_nrows()
  nztota = a%get_nzeros()
  select case(sm%variant)
-  case(amg_poly_lottes_)
+  case(amg_cheb_4_)
    ! do nothing
-  case(amg_poly_lottes_beta_)
+  case(amg_cheb_4_opt_)
    if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
      call psb_realloc(sm%pdegree,sm%poly_beta,info)
      sm%poly_beta(1:sm%pdegree) = amg_d_poly_beta_mat(1:sm%pdegree,sm%pdegree)
@@ -87,7 +87,7 @@ subroutine amg_s_poly_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
           & a_err='invalid sm%degree for poly_beta')
      goto 9999
    end if
-  case(amg_poly_new_)
+  case(amg_cheb_1_opt_)

    if ((1<=sm%pdegree).and.(sm%pdegree<=30)) then
      !Ok
@@ -58,11 +58,11 @@ subroutine amg_s_poly_smoother_cseti(sm,what,val,info,idx)
    sm%pdegree = val
  case('POLY_VARIANT')
    select case(val)
-    case(amg_poly_lottes_,amg_poly_lottes_beta_,amg_poly_new_)
+    case(amg_cheb_4_,amg_cheb_4_opt_,amg_cheb_1_opt_)
      sm%variant = val
    case default
-      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_poly_lottes_',val
-      sm%variant = amg_poly_lottes_
+      write(0,*) 'Invalid choice for POLY_VARIANT, defaulting to amg_cheb_4_',val
+      sm%variant = amg_cheb_4_
    end select
  case('POLY_RHO_ESTIMATE')
    select case(val)
@@ -77,17 +77,17 @@ subroutine amg_s_poly_smoother_descr(sm,info,iout,coarse,prefix)

  write(iout_,*) trim(prefix_), '  Polynomial smoother   '
  select case(sm%variant)
-  case(amg_poly_lottes_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES'
+  case(amg_cheb_4_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4'
    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
-  case(amg_poly_lottes_beta_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_LOTTES_BETA'
+  case(amg_cheb_4_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_4_OPT'
    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
    if (allocated(sm%poly_beta)) write(iout_,*) trim(prefix_), '   Coefficients:        ',sm%poly_beta(1:sm%pdegree)
-  case(amg_poly_new_)
-    write(iout_,*) trim(prefix_), '        variant:        ','POLY_NEW'
+  case(amg_cheb_1_opt_)
+    write(iout_,*) trim(prefix_), '        variant:        ','CHEB_1_OPT'
    write(iout_,*) trim(prefix_), '         Degree:        ',sm%pdegree
    write(iout_,*) trim(prefix_), '         rho_ba:        ',sm%rho_ba
    write(iout_,*) trim(prefix_), '    Coefficient:        ',sm%cf_a
@@ -40,7 +40,7 @@ subroutine amg_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_z_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_z_jac_smoother, amg_protect_name => amg_z_jac_smoother_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -53,7 +53,6 @@ subroutine amg_c_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
  class(psb_i_base_vect_type), intent(in), optional   :: imold
  ! Local variables
  integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
  character(len=20)   :: name='c_ilu_solver_bld', ch_err
@@ -40,7 +40,7 @@ subroutine amg_c_jac_solver_apply(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_c_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_c_jac_solver, amg_protect_name => amg_c_jac_solver_apply
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_c_jac_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_c_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_c_jac_solver, amg_protect_name => amg_c_jac_solver_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -169,7 +169,7 @@ subroutine amg_c_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
     & trans,work,wv,info,init,initu)

  use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use amg_c_krm_solver, amg_protect_name => amg_c_krm_solver_apply_vect

  Implicit None
@@ -53,7 +53,6 @@ subroutine amg_d_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
  class(psb_i_base_vect_type), intent(in), optional   :: imold
  ! Local variables
  integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
  character(len=20)   :: name='d_ilu_solver_bld', ch_err
@@ -40,7 +40,7 @@ subroutine amg_d_jac_solver_apply(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_d_jac_solver, amg_protect_name => amg_d_jac_solver_apply
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_d_jac_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_d_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_d_jac_solver, amg_protect_name => amg_d_jac_solver_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -169,7 +169,7 @@ subroutine amg_d_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
     & trans,work,wv,info,init,initu)

  use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use amg_d_krm_solver, amg_protect_name => amg_d_krm_solver_apply_vect

  Implicit None
@@ -53,7 +53,6 @@ subroutine amg_s_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
  class(psb_i_base_vect_type), intent(in), optional   :: imold
  ! Local variables
  integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
  character(len=20)   :: name='s_ilu_solver_bld', ch_err
@@ -40,7 +40,7 @@ subroutine amg_s_jac_solver_apply(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_s_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_s_jac_solver, amg_protect_name => amg_s_jac_solver_apply
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_s_jac_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_s_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_s_jac_solver, amg_protect_name => amg_s_jac_solver_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -169,7 +169,7 @@ subroutine amg_s_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
     & trans,work,wv,info,init,initu)

  use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use amg_s_krm_solver, amg_protect_name => amg_s_krm_solver_apply_vect

  Implicit None
@@ -53,7 +53,6 @@ subroutine amg_z_ilu_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
  class(psb_i_base_vect_type), intent(in), optional   :: imold
  ! Local variables
  integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota, psb_fctype
-!!$    complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_)   :: np, me, i, err_act, debug_unit, debug_level
  character(len=20)   :: name='z_ilu_solver_bld', ch_err
@@ -40,7 +40,7 @@ subroutine amg_z_jac_solver_apply(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_z_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_z_jac_solver, amg_protect_name => amg_z_jac_solver_apply
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -40,7 +40,7 @@ subroutine amg_z_jac_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,&

  use psb_base_mod
  use amg_z_diag_solver
-  use psb_base_krylov_conv_mod, only : log_conv
+  use psb_base_linsolve_conv_mod, only : log_conv
  use amg_z_jac_solver, amg_protect_name => amg_z_jac_solver_apply_vect
  implicit none
  type(psb_desc_type), intent(in)                 :: desc_data
@@ -169,7 +169,7 @@ subroutine amg_z_krm_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
     & trans,work,wv,info,init,initu)

  use psb_base_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use amg_z_krm_solver, amg_protect_name => amg_z_krm_solver_apply_vect

  Implicit None
@@ -264,7 +264,7 @@ contains
       & ah,ph,bh,xh,cdh,options) bind(c) result(res)
    use psb_base_mod
    use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
    use psb_prec_cbind_mod
    use psb_dkrylov_cbind_mod
    implicit none
@@ -285,7 +285,7 @@ contains
       & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
    use psb_base_mod
    use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
    use psb_objhandle_mod
    use psb_prec_cbind_mod
    use psb_base_string_cbind_mod
@@ -264,7 +264,7 @@ contains
       & ah,ph,bh,xh,cdh,options) bind(c) result(res)
    use psb_base_mod
    use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
    use psb_prec_cbind_mod
    use psb_zkrylov_cbind_mod
    implicit none
@@ -285,7 +285,7 @@ contains
       & ah,ph,bh,xh,eps,cdh,itmax,iter,err,itrace,irst,istop) bind(c) result(res)
    use psb_base_mod
    use psb_prec_mod
-    use psb_krylov_mod
+    use psb_linsolve_mod
    use psb_objhandle_mod
    use psb_prec_cbind_mod
    use psb_base_string_cbind_mod
@@ -9,8 +9,8 @@ HERE=.
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(LIBDIR) $(FMFLAG)$(PSBLAS_INCDIR)
 #PSBLAS_LIBS= -L$(PSBLAS_LIBDIR) -L$(LIBDIR)  $(CPSBLAS_LIB) $(PSBLAS_LIB)
 # -lpsb_krylov_cbind -lpsb_prec_cbind -lpsb_base_cbind 
-PSBC_LIBS= -L$(PSBLAS_LIBDIR) -lpsb_cbind  -lpsb_krylov -lpsb_prec 
-MLDC_LIBS=-L$(LIBDIR) -lmld_cbind  -lmld_prec 
+PSBC_LIBS= -L$(PSBLAS_LIBDIR) -lpsb_cbind  -lpsb_linsolve -lpsb_prec 
+AMGC_LIBS=-L$(LIBDIR) -lamg_cbind  -lamg_prec 
 #
 # Compilers and such
 #
@@ -23,15 +23,15 @@ EXEDIR=./runs
 #UMFLIBS=-lumfpack -lamd -lcholmod -lcolamd -lcamd -lccolamd -L/usr/include/suitesparse
 #UMFFLAGS=-DHave_UMF_ -I/usr/include/suitesparse

- all:  mldec
+ all:  amgec

-mldec: mldec.o
-	$(MPFC) mldec.o  -o mldec $(MLDC_LIBS)  $(PSBC_LIBS) $(PSBCLDLIBS) $(PSBLAS_LIBS) \
+amgec: amgec.o
+	$(MPFC) amgec.o  -o amgec $(AMGC_LIBS)  $(PSBC_LIBS) $(PSBCLDLIBS) $(PSBLAS_LIBS) \
 	  $(UMFLIBS) $(PSBLDLIBS) $(LDLIBS) -lm -lgfortran
 # \
 #	-lifcore -lifcoremt -lguide -limf -lirc -lintlc -lcxaguard -L/opt/intel/fc/10.0.023/lib/ -lm

-	/bin/mv mldec $(EXEDIR)
+	/bin/mv amgec $(EXEDIR)

 .f90.o:
 	$(MPFC) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<
@@ -40,13 +40,13 @@ mldec: mldec.o


 clean: 
-	/bin/rm -f mldec.o  $(EXEDIR)/mldec
+	/bin/rm -f amgec.o  $(EXEDIR)/amgec
 verycleanlib: 
 	(cd ../..; make veryclean)
 lib:
 	(cd ../../; make library)

 tests: all
-	cd runs ; ./mldec < mlde.inp
+	cd runs ; ./amgec < amge.inp


@@ -77,7 +77,7 @@
 #include <math.h>

 #include "psb_base_cbind.h"
-#include "mld_cbind.h"
+#include "amg_cbind.h"


 double  a1(double x, double y, double  z)
@@ -123,7 +123,7 @@ double g(double x, double y, double z)

 #define NBMAX       20

-psb_i_t matgen(psb_i_t ictxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
+psb_i_t matgen(psb_c_ctxt cctxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
 	       psb_c_dspmat *ah,psb_c_descriptor *cdh,
 	       psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
 {
@@ -135,7 +135,7 @@ psb_i_t matgen(psb_i_t ictxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
  psb_l_t irow[10*NBMAX], icol[10*NBMAX];

  info = 0;
-  psb_c_info(ictxt,&iam,&np);
+  psb_c_info(cctxt,&iam,&np);
  deltah = (double) 1.0/(idim+1);
  sqdeltah = deltah*deltah;
  deltah2  = 2.0* deltah;
@@ -253,11 +253,12 @@ void get_hparm(FILE *fp, char *val)

 int main(int argc, char *argv[])
 {
-  psb_i_t ictxt, iam, np;
+  psb_c_ctxt *cctxt;
+  psb_i_t iam, np;
  char methd[40], ptype[40], afmt[8];
  psb_i_t nparms;
  psb_i_t idim,info,istop,itmax,itrace,irst,iter,ret;
-  mld_c_dprec *ph;
+  amg_c_dprec *ph;
  psb_c_dspmat *ah;
  psb_c_dvector *bh, *xh, *rh; 
  psb_i_t nb,nlr, nl;
@@ -269,12 +270,13 @@ int main(int argc, char *argv[])
  psb_c_descriptor *cdh;
  FILE *vectfile;
  
-  ictxt = psb_c_init();
-  psb_c_info(ictxt,&iam,&np);
+  cctxt = psb_c_new_ctxt();
+  psb_c_init(cctxt);
+  psb_c_info(*cctxt,&iam,&np);
  fprintf(stdout,"Initialization: am %d of %d\n",iam,np);

  fflush(stdout);
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
  if (iam == 0) {
    get_iparm(stdin,&nparms);
    get_hparm(stdin,methd);
@@ -287,17 +289,17 @@ int main(int argc, char *argv[])
    get_iparm(stdin,&irst);
  } 
  /* Now broadcast the values, and check they're OK */
-  psb_c_ibcast(ictxt,1,&nparms,0);
-  psb_c_hbcast(ictxt,methd,0);
-  psb_c_hbcast(ictxt,ptype,0);
-  psb_c_hbcast(ictxt,afmt,0);
-  psb_c_ibcast(ictxt,1,&idim,0);
-  psb_c_ibcast(ictxt,1,&istop,0);
-  psb_c_ibcast(ictxt,1,&itmax,0);
-  psb_c_ibcast(ictxt,1,&itrace,0);
-  psb_c_ibcast(ictxt,1,&irst,0);
+  psb_c_ibcast(*cctxt,1,&nparms,0);
+  psb_c_hbcast(*cctxt,methd,0);
+  psb_c_hbcast(*cctxt,ptype,0);
+  psb_c_hbcast(*cctxt,afmt,0);
+  psb_c_ibcast(*cctxt,1,&idim,0);
+  psb_c_ibcast(*cctxt,1,&istop,0);
+  psb_c_ibcast(*cctxt,1,&itmax,0);
+  psb_c_ibcast(*cctxt,1,&itrace,0);
+  psb_c_ibcast(*cctxt,1,&irst,0);
  
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);

  cdh=psb_c_new_descriptor();
  psb_c_set_index_base(0);
@@ -310,15 +312,15 @@ int main(int argc, char *argv[])
  fprintf(stderr,"%d: Input data %d %ld %d %d\n",iam,idim,ng,nb, nl);
  if ((vl=malloc(nb*sizeof(psb_l_t)))==NULL) {
    fprintf(stderr,"On %d: malloc failure\n",iam);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }
  i = ((psb_l_t)iam) * nb;
  for (k=0; k<nl; k++)
    vl[k] = i+k; 

-  if ((info=psb_c_cdall_vl(nl,vl,ictxt,cdh))!=0) {
+  if ((info=psb_c_cdall_vl(nl,vl,*cctxt,cdh))!=0) {
    fprintf(stderr,"From cdall: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }

  bh  = psb_c_new_dvector();
@@ -337,25 +339,25 @@ int main(int argc, char *argv[])

  
  /* Matrix generation  */
-  if (matgen(ictxt,nl,idim,vl,ah,cdh,xh,bh,rh) != 0) {
+  if (matgen(*cctxt,nl,idim,vl,ah,cdh,xh,bh,rh) != 0) {
    fprintf(stderr,"Error during matrix build loop\n");
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }    
-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
  /* Set up the preconditioner */ 
-  ph  = mld_c_dprec_new();
-  mld_c_dprecinit(ictxt,ph,ptype);
-  mld_c_dprecseti(ph,"SMOOTHER_SWEEPS",2);
-  mld_c_dprecseti(ph,"SUB_FILLIN",1);
-  mld_c_dprecsetc(ph,"COARSE_SOLVE","BJAC");
-  mld_c_dprecsetc(ph,"COARSE_SUBSOLVE","ILU");
-  mld_c_dprecseti(ph,"COARSE_FILLIN",0);
-  if ((ret=mld_c_dhierarchy_build(ah,cdh,ph))!=0)
+  ph  = amg_c_dprec_new();
+  amg_c_dprecinit(*cctxt,ph,ptype);
+  amg_c_dprecseti(ph,"SMOOTHER_SWEEPS",2);
+  amg_c_dprecseti(ph,"SUB_FILLIN",1);
+  amg_c_dprecsetc(ph,"COARSE_SOLVE","BJAC");
+  amg_c_dprecsetc(ph,"COARSE_SUBSOLVE","ILU");
+  amg_c_dprecseti(ph,"COARSE_FILLIN",0);
+  if ((ret=amg_c_dhierarchy_build(ah,cdh,ph))!=0)
    fprintf(stderr,"From hierarchy_build: %d\n",ret);
-  if ((ret=mld_c_dsmoothers_build(ah,cdh,ph))!=0)
+  if ((ret=amg_c_dsmoothers_build(ah,cdh,ph))!=0)
    fprintf(stderr,"From smoothers_build: %d\n",ret);

-  psb_c_barrier(ictxt);
+  psb_c_barrier(*cctxt);
  /* Set up the solver options */ 
  psb_c_DefaultSolverOptions(&options);
  options.eps    = 1.e-6;
@@ -365,7 +367,7 @@ int main(int argc, char *argv[])
  options.istop  = istop;
  psb_c_seterraction_ret();
  t1=psb_c_wtime();
-  ret=mld_c_dkrylov(methd,ah,ph,bh,xh,cdh,&options); 
+  ret=amg_c_dkrylov(methd,ah,ph,bh,xh,cdh,&options); 
  t2=psb_c_wtime();
  iter = options.iter;
  err  = options.err;
@@ -413,20 +415,20 @@ int main(int argc, char *argv[])
  /* Clean up memory */ 
  if ((info=psb_c_dgefree(xh,cdh))!=0) {
    fprintf(stderr,"From dgefree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }
  if ((info=psb_c_dgefree(bh,cdh))!=0) {
    fprintf(stderr,"From dgefree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }
  if ((info=psb_c_dgefree(rh,cdh))!=0) {
    fprintf(stderr,"From dgefree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }

  if ((info=psb_c_cdfree(cdh))!=0) {
    fprintf(stderr,"From cdfree: %d\nBailing out\n",info);
-    psb_c_abort(ictxt);
+    psb_c_abort(*cctxt);
  }
  //fprintf(stderr,"pointer  from cdfree: %p\n",cdh->descriptor);
 
@@ -440,6 +442,7 @@ int main(int argc, char *argv[])

  if (iam == 0) fprintf(stderr,"program completed successfully\n");

-  psb_c_barrier(ictxt);
-  psb_c_exit(ictxt);
+  psb_c_barrier(*cctxt);
+  psb_c_exit(*cctxt);
+  free(cctxt);
 }
@@ -409,7 +409,7 @@ save_LDFLAGS=$LDFLAGS;
 ## dnl AC_MSG_NOTICE([psblas dir $pac_cv_psblas_dir])
 ##  PSBLAS_LIBS="-L$pac_cv_psblas_dir/lib"
 ## fi
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"

 dnl ac_compile='${MPIFC-$FC} -c -o conftest${ac_objext} $FMFLAG$PSBLAS_DIR/include $FMFLAG$PSBLAS_DIR/lib conftest.$ac_ext  1>&5'
@@ -484,7 +484,7 @@ dnl AC_MSG_NOTICE([psblas dir $pac_cv_psblas_dir])
 PSBLAS_INCLUDES="$FMFLAG$pac_cv_psblas_dir/modules $PSBLAS_INCLUDES"
 fi
 FCFLAGS=" $PSBLAS_INCLUDES $save_FCFLAGS"
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"

 dnl ac_compile='${MPIFC-$FC} -c -o conftest${ac_objext} $FMFLAG$PSBLAS_DIR/include $FMFLAG$PSBLAS_DIR/lib conftest.$ac_ext  1>&5'
@@ -7417,7 +7417,7 @@ then
   FIFLAG="-I"
   BASEMODNAME=PSB_BASE_MOD
   PRECMODNAME=PSB_PREC_MOD
-   METHDMODNAME=PSB_KRYLOV_MOD
+   METHDMODNAME=PSB_LINSOLVE_MOD
   UTILMODNAME=PSB_UTIL_MOD

 else
@@ -7550,7 +7550,7 @@ printf "%s\n" "$ax_cv_f90_modflag" >&6; }
   FIFLAG=-I
   BASEMODNAME=psb_base_mod
   PRECMODNAME=psb_prec_mod
-   METHDMODNAME=psb_krylov_mod
+   METHDMODNAME=psb_linsolve_mod
   UTILMODNAME=psb_util_mod
 fi

@@ -7683,7 +7683,7 @@ save_LDFLAGS=$LDFLAGS;
 ## dnl AC_MSG_NOTICE([psblas dir $pac_cv_psblas_dir])
 ##  PSBLAS_LIBS="-L$pac_cv_psblas_dir/lib"
 ## fi
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base -L$PSBLAS_LIBDIR"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"

 ac_link='${MPIFC-$FC} -o conftest${ac_exeext} $FCFLAGS conftest.$ac_ext $LDFLAGS $LIBS  1>&5'
@@ -7787,7 +7787,7 @@ elif test "x$pac_cv_psblas_dir" != "x"; then
 PSBLAS_INCLUDES="$FMFLAG$pac_cv_psblas_dir/modules $PSBLAS_INCLUDES"
 fi
 FCFLAGS=" $PSBLAS_INCLUDES $save_FCFLAGS"
-PSBLAS_LIBS="-lpsb_krylov -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
+PSBLAS_LIBS="-lpsb_linsolve -lpsb_prec -lpsb_util -lpsb_base $PSBLAS_LIBS"
 LDFLAGS=" $PSBLAS_LIBS $save_LDFLAGS"


@@ -525,7 +525,7 @@ then
   FIFLAG="-I"
   BASEMODNAME=PSB_BASE_MOD
   PRECMODNAME=PSB_PREC_MOD
-   METHDMODNAME=PSB_KRYLOV_MOD
+   METHDMODNAME=PSB_LINSOLVE_MOD
   UTILMODNAME=PSB_UTIL_MOD

 else 
@@ -536,7 +536,7 @@ else
   FIFLAG=-I
   BASEMODNAME=psb_base_mod
   PRECMODNAME=psb_prec_mod
-   METHDMODNAME=psb_krylov_mod
+   METHDMODNAME=psb_linsolve_mod
   UTILMODNAME=psb_util_mod
 fi

@@ -3,7 +3,7 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_krylov -lamg_prec -lpsb_prec
+AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).

 DFSOBJS=amg_df_sample.o data_input.o
@@ -38,7 +38,7 @@
 program amg_cf_sample
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -38,7 +38,7 @@
 program amg_df_sample
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -38,7 +38,7 @@
 program amg_sf_sample
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -38,7 +38,7 @@
 program amg_zf_sample
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -3,37 +3,42 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_krylov -lamg_prec -lpsb_prec 
+AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec 
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).

 LINKOPT=
 EXEDIR=./runs
-DGEN2D=amg_d_pde2d_base_mod.o amg_d_pde2d_exp_mod.o amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
-DGEN3D=amg_d_pde3d_base_mod.o amg_d_pde3d_exp_mod.o amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o 	
-SGEN2D=amg_s_pde2d_base_mod.o amg_s_pde2d_exp_mod.o amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
-SGEN3D=amg_s_pde3d_base_mod.o amg_s_pde3d_exp_mod.o amg_s_pde3d_gauss_mod.o amg_s_pde3d_box_mod.o 	
+DGEN2D=amg_d_pde2d_poisson_mod.o amg_d_pde2d_exp_mod.o \
+  amg_d_pde2d_gauss_mod.o amg_d_pde2d_box_mod.o
+DGEN3D=amg_d_pde3d_poisson_mod.o amg_d_pde3d_exp_mod.o \
+  amg_d_pde3d_gauss_mod.o amg_d_pde3d_box_mod.o 	
+SGEN2D=amg_s_pde2d_poisson_mod.o amg_s_pde2d_exp_mod.o \
+  amg_s_pde2d_gauss_mod.o amg_s_pde2d_box_mod.o
+SGEN3D=amg_s_pde3d_poisson_mod.o amg_s_pde3d_exp_mod.o \
+  amg_s_pde3d_gauss_mod.o amg_s_pde3d_box_mod.o 	

 all: amg_s_pde3d amg_d_pde3d amg_s_pde2d amg_d_pde2d

 amg_d_pde3d: amg_d_pde3d.o amg_d_genpde_mod.o $(DGEN3D) data_input.o
-	$(FLINK) $(LINKOPT) amg_d_pde3d.o  amg_d_genpde_mod.o $(DGEN3D) data_input.o -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_d_pde3d.o  amg_d_genpde_mod.o $(DGEN3D) data_input.o \
+         -o amg_d_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_d_pde3d $(EXEDIR)

 amg_s_pde3d: amg_s_pde3d.o amg_s_genpde_mod.o $(SGEN3D) data_input.o
-	$(FLINK) $(LINKOPT) amg_s_pde3d.o  amg_s_genpde_mod.o  $(SGEN3D) data_input.o -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_s_pde3d.o  amg_s_genpde_mod.o  $(SGEN3D) data_input.o \
+         -o amg_s_pde3d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_s_pde3d $(EXEDIR)

 amg_d_pde2d: amg_d_pde2d.o amg_d_genpde_mod.o  $(DGEN2D) data_input.o
-	$(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_d_pde2d.o amg_d_genpde_mod.o $(DGEN2D) data_input.o \
+         -o amg_d_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_d_pde2d $(EXEDIR)

 amg_s_pde2d: amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o
-	$(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
+	$(FLINK) $(LINKOPT) amg_s_pde2d.o amg_s_genpde_mod.o $(SGEN2D) data_input.o \
+          -o amg_s_pde2d $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
 	/bin/mv amg_s_pde2d $(EXEDIR)

-amg_d_pde3d_rebld: amg_d_pde3d_rebld.o  data_input.o
-	$(FLINK) $(LINKOPT) amg_d_pde3d_rebld.o  data_input.o -o amg_d_pde3d_rebld $(AMG_LIBS) $(PSBLAS_LIBS) $(LDLIBS)
-	/bin/mv amg_d_pde3d_rebld $(EXEDIR)

 amg_d_pde3d.o amg_s_pde3d.o amg_d_pde2d.o amg_s_pde2d.o: data_input.o

@@ -42,6 +47,11 @@ amg_s_pde3d.o: amg_s_genpde_mod.o $(SGEN3D)
 amg_d_pde2d.o: amg_d_genpde_mod.o $(DGEN2D)
 amg_s_pde2d.o: amg_s_genpde_mod.o $(SGEN2D)

+amg_d_genpde_mod.o: $(DGEN3D)
+amg_s_genpde_mod.o: $(SGEN3D)
+amg_d_genpde_mod.o: $(DGEN2D)
+amg_s_genpde_mod.o: $(SGEN2D)
+
 check: all
 	cd runs && ./amg_d_pde2d <amg_pde2d.inp && ./amg_s_pde2d<amg_pde2d.inp

@@ -66,10 +66,10 @@
 program amg_d_pde2d
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
-  use amg_d_pde2d_base_mod
+  use amg_d_pde2d_poisson_mod
  use amg_d_pde2d_exp_mod
  use amg_d_pde2d_box_mod
  use amg_d_pde2d_gauss_mod
@@ -81,7 +81,8 @@ program amg_d_pde2d

  ! input parameters
  character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
  integer(psb_ipk_) :: idim
  integer(psb_epk_) :: system_size

@@ -146,6 +147,8 @@ program amg_d_pde2d
    integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
    integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
    character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue    ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr         ! number of overlap layers
    character(len=32)  :: restr        ! restriction over application of AS
    character(len=32)  :: prol         ! prolongation over application of AS
@@ -162,6 +165,8 @@ program amg_d_pde2d
    integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
    integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
    character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue2   ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr2        ! number of overlap layers
    character(len=32)  :: restr2       ! restriction  over application of AS
    character(len=32)  :: prol2        ! prolongation over application of AS
@@ -244,18 +249,22 @@ program amg_d_pde2d
  call psb_barrier(ctxt)
  t1 = psb_wtime()
  select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,b1_base,b2_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,&
+         & b1_poisson,b2_poisson,c_poisson,g_poisson,info)
  case("EXP")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info)
+         & a1_exp,a2_exp,&
+         & b1_exp,b2_exp,c_exp,g_exp,info)
  case("BOX")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_box,a2_box,b1_box,b2_box,c_box,g_box,info)
+         & a1_box,a2_box,&
+         & b1_box,b2_box,c_box,g_box,info)
  case("GAUSS")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_gauss,a2_gauss,b1_gauss,b2_gauss,c_gauss,g_gauss,info)
+         & a1_gauss,a2_gauss,&
+         & b1_gauss,b2_gauss,c_gauss,g_gauss,info)
  case default
    info=psb_err_from_subroutine_
    ch_err='amg_gen_pdecoeff'
@@ -344,7 +353,12 @@ program amg_d_pde2d
    call prec%set('smoother_sweeps', p_choice%jsweeps,    info)
    call prec%set('poly_degree',     p_choice%degree,    info)
    call prec%set('poly_variant',    p_choice%pvariant,  info)
-    
+    if (p_choice%prhovalue > dzero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
    select case (psb_toupper(p_choice%smther))
    case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
      ! do nothing
@@ -376,6 +390,12 @@ program amg_d_pde2d
      call prec%set('smoother_sweeps', p_choice%jsweeps2,  info,pos='post')
      call prec%set('poly_degree',     p_choice%degree2,   info,pos='post')
      call prec%set('poly_variant',    p_choice%pvariant2, info,pos='post')
+      if (p_choice%prhovalue > dzero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
      select case (psb_toupper(p_choice%smther2))
      case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
        ! do nothing
@@ -452,9 +472,16 @@ program amg_d_pde2d
  !
  call psb_barrier(ctxt)
  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,&
-       & istop=s_choice%istopc,irst=s_choice%irst)
+  if (psb_toupper(trim(s_choice%kmethd)) == 'RICHARDSON') then
+    call psb_richardson(a,prec,b,x,s_choice%eps,&
+         & desc_a,info,itmax=s_choice%itmax,iter=iter,&
+         & err=err,itrace=s_choice%itrace,&
+         & istop=s_choice%istopc)
+  else
+    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,&
+         & istop=s_choice%istopc,irst=s_choice%irst)
+  end if
  call psb_barrier(ctxt)
  tslv = psb_wtime() - t1

@@ -592,7 +619,9 @@ contains
      call read_data(prec%smther,inp_unit)     ! smoother type
      call read_data(prec%jsweeps,inp_unit)    ! (pre-)smoother / 1-lev prec sweeps
      call read_data(prec%degree,inp_unit)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr,inp_unit)       ! number of overlap layers
      call read_data(prec%restr,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol,inp_unit)       ! prolongation over application of AS
@@ -607,6 +636,8 @@ contains
      call read_data(prec%jsweeps2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%degree2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%pvariant2,inp_unit)   !
+      call read_data(prec%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr2,inp_unit)       ! number of overlap layers
      call read_data(prec%restr2,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol2,inp_unit)       ! prolongation over application of AS
@@ -679,6 +710,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps)
    call psb_bcast(ctxt,prec%degree)
    call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
    call psb_bcast(ctxt,prec%novr)
    call psb_bcast(ctxt,prec%restr)
    call psb_bcast(ctxt,prec%prol)
@@ -693,6 +726,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps2)
    call psb_bcast(ctxt,prec%degree2)
    call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
    call psb_bcast(ctxt,prec%novr2)
    call psb_bcast(ctxt,prec%restr2)
    call psb_bcast(ctxt,prec%prol2)
@@ -34,56 +34,56 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_d_pde2d_base_mod
+module amg_d_pde2d_poisson_mod
  use psb_base_mod, only : psb_dpk_, dzero, done
  real(psb_dpk_), save, private :: epsilon=done/80
 contains
-  subroutine pde_set_parm2d_base(dat)
+  subroutine pde_set_parm2d_poisson(dat)
    real(psb_dpk_), intent(in) :: dat
    epsilon = dat
-  end subroutine pde_set_parm2d_base
+  end subroutine pde_set_parm2d_poisson
  !
  ! functions parametrizing the differential equation
  !
-  function b1_base(x,y)
+  function b1_poisson(x,y)
    implicit none 
-    real(psb_dpk_) :: b1_base
+    real(psb_dpk_) :: b1_poisson
    real(psb_dpk_), intent(in) :: x,y
-    b1_base = dzero/1.414_psb_dpk_
-  end function b1_base
-  function b2_base(x,y)
+    b1_poisson = dzero
+  end function b1_poisson
+  function b2_poisson(x,y)
    implicit none 
-    real(psb_dpk_) ::  b2_base
+    real(psb_dpk_) ::  b2_poisson
    real(psb_dpk_), intent(in) :: x,y
-    b2_base = dzero/1.414_psb_dpk_
-  end function b2_base
-  function c_base(x,y)
+    b2_poisson = dzero
+  end function b2_poisson
+  function c_poisson(x,y)
    implicit none 
-    real(psb_dpk_) ::  c_base
+    real(psb_dpk_) ::  c_poisson
    real(psb_dpk_), intent(in) :: x,y
-    c_base = dzero
-  end function c_base
-  function a1_base(x,y)
+    c_poisson = dzero
+  end function c_poisson
+  function a1_poisson(x,y)
    implicit none 
-    real(psb_dpk_) ::  a1_base
+    real(psb_dpk_) ::  a1_poisson
    real(psb_dpk_), intent(in) :: x,y
-    a1_base=done*epsilon
-  end function a1_base
-  function a2_base(x,y)
+    a1_poisson=done*epsilon
+  end function a1_poisson
+  function a2_poisson(x,y)
    implicit none 
-    real(psb_dpk_) ::  a2_base
+    real(psb_dpk_) ::  a2_poisson
    real(psb_dpk_), intent(in) :: x,y
-    a2_base=done*epsilon
-  end function a2_base
-  function g_base(x,y)
+    a2_poisson=done*epsilon
+  end function a2_poisson
+  function g_poisson(x,y)
    implicit none 
-    real(psb_dpk_) ::  g_base
+    real(psb_dpk_) ::  g_poisson
    real(psb_dpk_), intent(in) :: x,y
-    g_base = dzero
+    g_poisson = dzero
    if (x == done) then
-      g_base = done
+      g_poisson = done
    else if (x == dzero) then
-      g_base = done
+      g_poisson = done
    end if
-  end function g_base
-end module amg_d_pde2d_base_mod
+  end function g_poisson
+end module amg_d_pde2d_poisson_mod
@@ -67,10 +67,10 @@
 program amg_d_pde3d
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
-  use amg_d_pde3d_base_mod
+  use amg_d_pde3d_poisson_mod
  use amg_d_pde3d_exp_mod
  use amg_d_pde3d_box_mod
  use amg_d_pde3d_gauss_mod
@@ -82,7 +82,8 @@ program amg_d_pde3d

  ! input parameters
  character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
  integer(psb_ipk_) :: idim
  integer(psb_epk_) :: system_size

@@ -147,6 +148,8 @@ program amg_d_pde3d
    integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
    integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
    character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue    ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr         ! number of overlap layers
    character(len=32)  :: restr        ! restriction over application of AS
    character(len=32)  :: prol         ! prolongation over application of AS
@@ -163,6 +166,8 @@ program amg_d_pde3d
    integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
    integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
    character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_dpk_)     :: prhovalue2   ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr2        ! number of overlap layers
    character(len=32)  :: restr2       ! restriction  over application of AS
    character(len=32)  :: prol2        ! prolongation over application of AS
@@ -246,18 +251,22 @@ program amg_d_pde3d
  call psb_barrier(ctxt)
  t1 = psb_wtime()
  select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,a3_poisson,&
+         & b1_poisson,b2_poisson,b3_poisson,c_poisson,g_poisson,info)
  case("EXP")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info)
+         & a1_exp,a2_exp,a3_exp,&
+         & b1_exp,b2_exp,b3_exp,c_exp,g_exp,info)
  case("BOX")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_box,a2_box,a3_box,b1_box,b2_box,b3_box,c_box,g_box,info)
+         & a1_box,a2_box,a3_box,&
+         & b1_box,b2_box,b3_box,c_box,g_box,info)
  case("GAUSS")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
+         & a1_gauss,a2_gauss,a3_gauss,&
+         & b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
  case default
    info=psb_err_from_subroutine_
    ch_err='amg_gen_pdecoeff'
@@ -348,7 +357,12 @@ program amg_d_pde3d
    call prec%set('smoother_sweeps', p_choice%jsweeps,    info)
    call prec%set('poly_degree',     p_choice%degree,    info)
    call prec%set('poly_variant',    p_choice%pvariant,  info)
-    
+    if (p_choice%prhovalue > dzero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
    select case (psb_toupper(p_choice%smther))
    case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
      ! do nothing
@@ -380,6 +394,12 @@ program amg_d_pde3d
      call prec%set('smoother_sweeps', p_choice%jsweeps2,  info,pos='post')
      call prec%set('poly_degree',     p_choice%degree2,   info,pos='post')
      call prec%set('poly_variant',    p_choice%pvariant2, info,pos='post')
+      if (p_choice%prhovalue > dzero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
      select case (psb_toupper(p_choice%smther2))
      case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
        ! do nothing
@@ -456,9 +476,16 @@ program amg_d_pde3d
  !
  call psb_barrier(ctxt)
  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,&
-       & istop=s_choice%istopc,irst=s_choice%irst)
+  if (psb_toupper(trim(s_choice%kmethd)) == 'RICHARDSON') then
+    call psb_richardson(a,prec,b,x,s_choice%eps,&
+         & desc_a,info,itmax=s_choice%itmax,iter=iter,&
+         & err=err,itrace=s_choice%itrace,&
+         & istop=s_choice%istopc)
+  else
+    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,&
+         & istop=s_choice%istopc,irst=s_choice%irst)
+  end if
  call psb_barrier(ctxt)
  tslv = psb_wtime() - t1

@@ -596,7 +623,9 @@ contains
      call read_data(prec%smther,inp_unit)     ! smoother type
      call read_data(prec%jsweeps,inp_unit)    ! (pre-)smoother / 1-lev prec sweeps
      call read_data(prec%degree,inp_unit)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr,inp_unit)       ! number of overlap layers
      call read_data(prec%restr,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol,inp_unit)       ! prolongation over application of AS
@@ -611,6 +640,8 @@ contains
      call read_data(prec%jsweeps2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%degree2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%pvariant2,inp_unit)   !
+      call read_data(prec%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr2,inp_unit)       ! number of overlap layers
      call read_data(prec%restr2,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol2,inp_unit)       ! prolongation over application of AS
@@ -683,6 +714,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps)
    call psb_bcast(ctxt,prec%degree)
    call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
    call psb_bcast(ctxt,prec%novr)
    call psb_bcast(ctxt,prec%restr)
    call psb_bcast(ctxt,prec%prol)
@@ -697,6 +730,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps2)
    call psb_bcast(ctxt,prec%degree2)
    call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
    call psb_bcast(ctxt,prec%novr2)
    call psb_bcast(ctxt,prec%restr2)
    call psb_bcast(ctxt,prec%prol2)
@@ -34,68 +34,68 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_d_pde3d_base_mod
+module amg_d_pde3d_poisson_mod
  use psb_base_mod, only : psb_dpk_, done, dzero
  real(psb_dpk_), save, private :: epsilon=done/80
 contains
-  subroutine pde_set_parm3d_base(dat)
+  subroutine pde_set_parm3d_poisson(dat)
    real(psb_dpk_), intent(in) :: dat
    epsilon = dat
-  end subroutine pde_set_parm3d_base
+  end subroutine pde_set_parm3d_poisson
  !
  ! functions parametrizing the differential equation
  !
-  function b1_base(x,y,z)
+  function b1_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) :: b1_base
+    real(psb_dpk_) :: b1_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    b1_base=dzero/sqrt(3.0_psb_dpk_)
-  end function b1_base
-  function b2_base(x,y,z)
+    b1_poisson=dzero
+  end function b1_poisson
+  function b2_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  b2_base
+    real(psb_dpk_) ::  b2_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    b2_base=dzero/sqrt(3.0_psb_dpk_)
-  end function b2_base
-  function b3_base(x,y,z)
+    b2_poisson=dzero
+  end function b2_poisson
+  function b3_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  b3_base
+    real(psb_dpk_) ::  b3_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    b3_base=dzero/sqrt(3.0_psb_dpk_)
-  end function b3_base
-  function c_base(x,y,z)
+    b3_poisson=dzero
+  end function b3_poisson
+  function c_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  c_base
+    real(psb_dpk_) ::  c_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    c_base=dzero
-  end function c_base
-  function a1_base(x,y,z)
+    c_poisson=dzero
+  end function c_poisson
+  function a1_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  a1_base
+    real(psb_dpk_) ::  a1_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    a1_base=epsilon
-  end function a1_base
-  function a2_base(x,y,z)
+    a1_poisson=epsilon
+  end function a1_poisson
+  function a2_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  a2_base
+    real(psb_dpk_) ::  a2_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    a2_base=epsilon
-  end function a2_base
-  function a3_base(x,y,z)
+    a2_poisson=epsilon
+  end function a2_poisson
+  function a3_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  a3_base
+    real(psb_dpk_) ::  a3_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    a3_base=epsilon
-  end function a3_base
-  function g_base(x,y,z)
+    a3_poisson=epsilon
+  end function a3_poisson
+  function g_poisson(x,y,z)
    implicit none 
-    real(psb_dpk_) ::  g_base
+    real(psb_dpk_) ::  g_poisson
    real(psb_dpk_), intent(in) :: x,y,z
-    g_base = dzero
+    g_poisson = dzero
    if (x == done) then
-      g_base = done
+      g_poisson = done
    else if (x == dzero) then
-      g_base = done
+      g_poisson = done
    end if
-  end function g_base
-end module amg_d_pde3d_base_mod
+  end function g_poisson
+end module amg_d_pde3d_poisson_mod
@@ -66,10 +66,10 @@
 program amg_s_pde2d
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
-  use amg_s_pde2d_base_mod
+  use amg_s_pde2d_poisson_mod
  use amg_s_pde2d_exp_mod
  use amg_s_pde2d_box_mod
  use amg_s_pde2d_gauss_mod
@@ -81,7 +81,8 @@ program amg_s_pde2d

  ! input parameters
  character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
  integer(psb_ipk_) :: idim
  integer(psb_epk_) :: system_size

@@ -146,6 +147,8 @@ program amg_s_pde2d
    integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
    integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
    character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue    ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr         ! number of overlap layers
    character(len=32)  :: restr        ! restriction over application of AS
    character(len=32)  :: prol         ! prolongation over application of AS
@@ -162,6 +165,8 @@ program amg_s_pde2d
    integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
    integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
    character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue2   ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr2        ! number of overlap layers
    character(len=32)  :: restr2       ! restriction  over application of AS
    character(len=32)  :: prol2        ! prolongation over application of AS
@@ -244,18 +249,22 @@ program amg_s_pde2d
  call psb_barrier(ctxt)
  t1 = psb_wtime()
  select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,b1_base,b2_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,&
+         & b1_poisson,b2_poisson,c_poisson,g_poisson,info)
  case("EXP")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_exp,a2_exp,b1_exp,b2_exp,c_exp,g_exp,info)
+         & a1_exp,a2_exp,&
+         & b1_exp,b2_exp,c_exp,g_exp,info)
  case("BOX")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_box,a2_box,b1_box,b2_box,c_box,g_box,info)
+         & a1_box,a2_box,&
+         & b1_box,b2_box,c_box,g_box,info)
  case("GAUSS")
    call amg_gen_pde2d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_gauss,a2_gauss,b1_gauss,b2_gauss,c_gauss,g_gauss,info)
+         & a1_gauss,a2_gauss,&
+         & b1_gauss,b2_gauss,c_gauss,g_gauss,info)
  case default
    info=psb_err_from_subroutine_
    ch_err='amg_gen_pdecoeff'
@@ -344,7 +353,12 @@ program amg_s_pde2d
    call prec%set('smoother_sweeps', p_choice%jsweeps,    info)
    call prec%set('poly_degree',     p_choice%degree,    info)
    call prec%set('poly_variant',    p_choice%pvariant,  info)
-    
+    if (p_choice%prhovalue > szero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
    select case (psb_toupper(p_choice%smther))
    case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
      ! do nothing
@@ -376,6 +390,12 @@ program amg_s_pde2d
      call prec%set('smoother_sweeps', p_choice%jsweeps2,  info,pos='post')
      call prec%set('poly_degree',     p_choice%degree2,   info,pos='post')
      call prec%set('poly_variant',    p_choice%pvariant2, info,pos='post')
+      if (p_choice%prhovalue > szero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
      select case (psb_toupper(p_choice%smther2))
      case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
        ! do nothing
@@ -452,9 +472,16 @@ program amg_s_pde2d
  !
  call psb_barrier(ctxt)
  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,&
-       & istop=s_choice%istopc,irst=s_choice%irst)
+  if (psb_toupper(trim(s_choice%kmethd)) == 'RICHARDSON') then
+    call psb_richardson(a,prec,b,x,s_choice%eps,&
+         & desc_a,info,itmax=s_choice%itmax,iter=iter,&
+         & err=err,itrace=s_choice%itrace,&
+         & istop=s_choice%istopc)
+  else
+    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,&
+         & istop=s_choice%istopc,irst=s_choice%irst)
+  end if
  call psb_barrier(ctxt)
  tslv = psb_wtime() - t1

@@ -592,7 +619,9 @@ contains
      call read_data(prec%smther,inp_unit)     ! smoother type
      call read_data(prec%jsweeps,inp_unit)    ! (pre-)smoother / 1-lev prec sweeps
      call read_data(prec%degree,inp_unit)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr,inp_unit)       ! number of overlap layers
      call read_data(prec%restr,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol,inp_unit)       ! prolongation over application of AS
@@ -607,6 +636,8 @@ contains
      call read_data(prec%jsweeps2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%degree2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%pvariant2,inp_unit)   !
+      call read_data(prec%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr2,inp_unit)       ! number of overlap layers
      call read_data(prec%restr2,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol2,inp_unit)       ! prolongation over application of AS
@@ -679,6 +710,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps)
    call psb_bcast(ctxt,prec%degree)
    call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
    call psb_bcast(ctxt,prec%novr)
    call psb_bcast(ctxt,prec%restr)
    call psb_bcast(ctxt,prec%prol)
@@ -693,6 +726,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps2)
    call psb_bcast(ctxt,prec%degree2)
    call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
    call psb_bcast(ctxt,prec%novr2)
    call psb_bcast(ctxt,prec%restr2)
    call psb_bcast(ctxt,prec%prol2)
@@ -34,56 +34,56 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_s_pde2d_base_mod
+module amg_s_pde2d_poisson_mod
  use psb_base_mod, only : psb_spk_, szero, sone
  real(psb_spk_), save, private :: epsilon=sone/80
 contains
-  subroutine pde_set_parm2d_base(dat)
+  subroutine pde_set_parm2d_poisson(dat)
    real(psb_spk_), intent(in) :: dat
    epsilon = dat
-  end subroutine pde_set_parm2d_base
+  end subroutine pde_set_parm2d_poisson
  !
  ! functions parametrizing the differential equation
  !
-  function b1_base(x,y)
+  function b1_poisson(x,y)
    implicit none 
-    real(psb_spk_) :: b1_base
+    real(psb_spk_) :: b1_poisson
    real(psb_spk_), intent(in) :: x,y
-    b1_base = szero/1.414_psb_spk_
-  end function b1_base
-  function b2_base(x,y)
+    b1_poisson = szero
+  end function b1_poisson
+  function b2_poisson(x,y)
    implicit none 
-    real(psb_spk_) ::  b2_base
+    real(psb_spk_) ::  b2_poisson
    real(psb_spk_), intent(in) :: x,y
-    b2_base = szero/1.414_psb_spk_
-  end function b2_base
-  function c_base(x,y)
+    b2_poisson = szero
+  end function b2_poisson
+  function c_poisson(x,y)
    implicit none 
-    real(psb_spk_) ::  c_base
+    real(psb_spk_) ::  c_poisson
    real(psb_spk_), intent(in) :: x,y
-    c_base = szero
-  end function c_base
-  function a1_base(x,y)
+    c_poisson = szero
+  end function c_poisson
+  function a1_poisson(x,y)
    implicit none 
-    real(psb_spk_) ::  a1_base
+    real(psb_spk_) ::  a1_poisson
    real(psb_spk_), intent(in) :: x,y
-    a1_base=sone*epsilon
-  end function a1_base
-  function a2_base(x,y)
+    a1_poisson=sone*epsilon
+  end function a1_poisson
+  function a2_poisson(x,y)
    implicit none 
-    real(psb_spk_) ::  a2_base
+    real(psb_spk_) ::  a2_poisson
    real(psb_spk_), intent(in) :: x,y
-    a2_base=sone*epsilon
-  end function a2_base
-  function g_base(x,y)
+    a2_poisson=sone*epsilon
+  end function a2_poisson
+  function g_poisson(x,y)
    implicit none 
-    real(psb_spk_) ::  g_base
+    real(psb_spk_) ::  g_poisson
    real(psb_spk_), intent(in) :: x,y
-    g_base = szero
+    g_poisson = szero
    if (x == sone) then
-      g_base = sone
+      g_poisson = sone
    else if (x == szero) then
-      g_base = sone
+      g_poisson = sone
    end if
-  end function g_base
-end module amg_s_pde2d_base_mod
+  end function g_poisson
+end module amg_s_pde2d_poisson_mod
@@ -67,10 +67,10 @@
 program amg_s_pde3d
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
-  use amg_s_pde3d_base_mod
+  use amg_s_pde3d_poisson_mod
  use amg_s_pde3d_exp_mod
  use amg_s_pde3d_box_mod
  use amg_s_pde3d_gauss_mod
@@ -82,7 +82,8 @@ program amg_s_pde3d

  ! input parameters
  character(len=20) :: kmethd, ptype
-  character(len=5)  :: afmt, pdecoeff
+  character(len=5)  :: afmt
+  character(len=32)  :: pdecoeff
  integer(psb_ipk_) :: idim
  integer(psb_epk_) :: system_size

@@ -147,6 +148,8 @@ program amg_s_pde3d
    integer(psb_ipk_)  :: jsweeps      ! (pre-)smoother / 1-lev prec. sweeps
    integer(psb_ipk_)  :: degree       ! degree for polynomial smoother
    character(len=32)  :: pvariant     ! polynomial  variant
+    character(len=32)  :: prhovariant  ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue    ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr         ! number of overlap layers
    character(len=32)  :: restr        ! restriction over application of AS
    character(len=32)  :: prol         ! prolongation over application of AS
@@ -163,6 +166,8 @@ program amg_s_pde3d
    integer(psb_ipk_)  :: jsweeps2     ! post-smoother sweeps
    integer(psb_ipk_)  :: degree2      ! degree for polynomial smoother
    character(len=32)  :: pvariant2    ! polynomial  variant
+    character(len=32)  :: prhovariant2 ! how to estimate rho(M^{-1}A)
+    real(psb_spk_)     :: prhovalue2   ! if previous is set value, we set it from this one
    integer(psb_ipk_)  :: novr2        ! number of overlap layers
    character(len=32)  :: restr2       ! restriction  over application of AS
    character(len=32)  :: prol2        ! prolongation over application of AS
@@ -246,18 +251,22 @@ program amg_s_pde3d
  call psb_barrier(ctxt)
  t1 = psb_wtime()
  select case(psb_toupper(trim(pdecoeff)))
-  case("CONST")
+  case("POISSON")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_base,a2_base,a3_base,b1_base,b2_base,b3_base,c_base,g_base,info)
+         & a1_poisson,a2_poisson,a3_poisson,&
+         & b1_poisson,b2_poisson,b3_poisson,c_poisson,g_poisson,info)
  case("EXP")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_exp,a2_exp,a3_exp,b1_exp,b2_exp,b3_exp,c_exp,g_exp,info)
+         & a1_exp,a2_exp,a3_exp,&
+         & b1_exp,b2_exp,b3_exp,c_exp,g_exp,info)
  case("BOX")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_box,a2_box,a3_box,b1_box,b2_box,b3_box,c_box,g_box,info)
+         & a1_box,a2_box,a3_box,&
+         & b1_box,b2_box,b3_box,c_box,g_box,info)
  case("GAUSS")
    call amg_gen_pde3d(ctxt,idim,a,b,x,desc_a,afmt,&
-        & a1_gauss,a2_gauss,a3_gauss,b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
+         & a1_gauss,a2_gauss,a3_gauss,&
+         & b1_gauss,b2_gauss,b3_gauss,c_gauss,g_gauss,info)
  case default
    info=psb_err_from_subroutine_
    ch_err='amg_gen_pdecoeff'
@@ -348,7 +357,12 @@ program amg_s_pde3d
    call prec%set('smoother_sweeps', p_choice%jsweeps,    info)
    call prec%set('poly_degree',     p_choice%degree,    info)
    call prec%set('poly_variant',    p_choice%pvariant,  info)
-    
+    if (p_choice%prhovalue > szero ) then
+            call prec%set('poly_rho_ba', p_choice%prhovalue, info)
+    else
+            call prec%set('poly_rho_estimate', p_choice%prhovariant, info)
+    end if
+
    select case (psb_toupper(p_choice%smther))
    case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
      ! do nothing
@@ -380,6 +394,12 @@ program amg_s_pde3d
      call prec%set('smoother_sweeps', p_choice%jsweeps2,  info,pos='post')
      call prec%set('poly_degree',     p_choice%degree2,   info,pos='post')
      call prec%set('poly_variant',    p_choice%pvariant2, info,pos='post')
+      if (p_choice%prhovalue > szero ) then
+         call prec%set('poly_rho_ba', p_choice%prhovalue2, info,pos='post')
+      else
+         call prec%set('poly_rho_estimate', p_choice%prhovariant2, info,pos='post')
+      end if
+
      select case (psb_toupper(p_choice%smther2))
      case ('GS','BWGS','FBGS','JACOBI','L1-JACOBI','L1-FBGS')
        ! do nothing
@@ -456,9 +476,16 @@ program amg_s_pde3d
  !
  call psb_barrier(ctxt)
  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,&
-       & istop=s_choice%istopc,irst=s_choice%irst)
+  if (psb_toupper(trim(s_choice%kmethd)) == 'RICHARDSON') then
+    call psb_richardson(a,prec,b,x,s_choice%eps,&
+         & desc_a,info,itmax=s_choice%itmax,iter=iter,&
+         & err=err,itrace=s_choice%itrace,&
+         & istop=s_choice%istopc)
+  else
+    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,&
+         & istop=s_choice%istopc,irst=s_choice%irst)
+  end if
  call psb_barrier(ctxt)
  tslv = psb_wtime() - t1

@@ -596,7 +623,9 @@ contains
      call read_data(prec%smther,inp_unit)     ! smoother type
      call read_data(prec%jsweeps,inp_unit)    ! (pre-)smoother / 1-lev prec sweeps
      call read_data(prec%degree,inp_unit)     ! (pre-)smoother / 1-lev prec sweeps
-      call read_data(prec%pvariant,inp_unit)   ! 
+      call read_data(prec%pvariant,inp_unit)   !
+      call read_data(prec%prhovariant,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr,inp_unit)       ! number of overlap layers
      call read_data(prec%restr,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol,inp_unit)       ! prolongation over application of AS
@@ -611,6 +640,8 @@ contains
      call read_data(prec%jsweeps2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%degree2,inp_unit)    ! (post-)smoother sweeps
      call read_data(prec%pvariant2,inp_unit)   !
+      call read_data(prec%prhovariant2,inp_unit)! how to estimate rho(M^{-1}A)
+      call read_data(prec%prhovalue2,inp_unit)  ! if previous is set value, we set it from this one
      call read_data(prec%novr2,inp_unit)       ! number of overlap layers
      call read_data(prec%restr2,inp_unit)      ! restriction  over application of AS
      call read_data(prec%prol2,inp_unit)       ! prolongation over application of AS
@@ -683,6 +714,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps)
    call psb_bcast(ctxt,prec%degree)
    call psb_bcast(ctxt,prec%pvariant)
+    call psb_bcast(ctxt,prec%prhovariant)
+    call psb_bcast(ctxt,prec%prhovalue)  
    call psb_bcast(ctxt,prec%novr)
    call psb_bcast(ctxt,prec%restr)
    call psb_bcast(ctxt,prec%prol)
@@ -697,6 +730,8 @@ contains
    call psb_bcast(ctxt,prec%jsweeps2)
    call psb_bcast(ctxt,prec%degree2)
    call psb_bcast(ctxt,prec%pvariant2)
+    call psb_bcast(ctxt,prec%prhovariant2)
+    call psb_bcast(ctxt,prec%prhovalue2)  
    call psb_bcast(ctxt,prec%novr2)
    call psb_bcast(ctxt,prec%restr2)
    call psb_bcast(ctxt,prec%prol2)
@@ -34,68 +34,68 @@
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
 !
-module amg_s_pde3d_base_mod
+module amg_s_pde3d_poisson_mod
  use psb_base_mod, only : psb_spk_, sone, szero
  real(psb_spk_), save, private :: epsilon=sone/80
 contains
-  subroutine pde_set_parm3d_base(dat)
+  subroutine pde_set_parm3d_poisson(dat)
    real(psb_spk_), intent(in) :: dat
    epsilon = dat
-  end subroutine pde_set_parm3d_base
+  end subroutine pde_set_parm3d_poisson
  !
  ! functions parametrizing the differential equation
  !
-  function b1_base(x,y,z)
+  function b1_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) :: b1_base
+    real(psb_spk_) :: b1_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    b1_base=szero/sqrt(3.0_psb_spk_)
-  end function b1_base
-  function b2_base(x,y,z)
+    b1_poisson=szero
+  end function b1_poisson
+  function b2_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  b2_base
+    real(psb_spk_) ::  b2_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    b2_base=szero/sqrt(3.0_psb_spk_)
-  end function b2_base
-  function b3_base(x,y,z)
+    b2_poisson=szero
+  end function b2_poisson
+  function b3_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  b3_base
+    real(psb_spk_) ::  b3_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    b3_base=szero/sqrt(3.0_psb_spk_)
-  end function b3_base
-  function c_base(x,y,z)
+    b3_poisson=szero
+  end function b3_poisson
+  function c_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  c_base
+    real(psb_spk_) ::  c_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    c_base=szero
-  end function c_base
-  function a1_base(x,y,z)
+    c_poisson=szero
+  end function c_poisson
+  function a1_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  a1_base
+    real(psb_spk_) ::  a1_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    a1_base=epsilon
-  end function a1_base
-  function a2_base(x,y,z)
+    a1_poisson=epsilon
+  end function a1_poisson
+  function a2_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  a2_base
+    real(psb_spk_) ::  a2_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    a2_base=epsilon
-  end function a2_base
-  function a3_base(x,y,z)
+    a2_poisson=epsilon
+  end function a2_poisson
+  function a3_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  a3_base
+    real(psb_spk_) ::  a3_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    a3_base=epsilon
-  end function a3_base
-  function g_base(x,y,z)
+    a3_poisson=epsilon
+  end function a3_poisson
+  function g_poisson(x,y,z)
    implicit none 
-    real(psb_spk_) ::  g_base
+    real(psb_spk_) ::  g_poisson
    real(psb_spk_), intent(in) :: x,y,z
-    g_base = szero
+    g_poisson = szero
    if (x == sone) then
-      g_base = sone
+      g_poisson = sone
    else if (x == szero) then
-      g_base = sone
+      g_poisson = sone
    end if
-  end function g_base
-end module amg_s_pde3d_base_mod
+  end function g_poisson
+end module amg_s_pde3d_poisson_mod
@@ -1,8 +1,10 @@
 %%%%%%%%%%%  General  arguments % Lines starting with % are ignored.
 CSR                         ! Storage format CSR COO JAD
-0200                        ! IDIM; domain size. Linear system size is IDIM**2
-CONST                       ! PDECOEFF: CONST, EXP, BOX, GAUSS Coefficients of the PDE
-CG                          ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
+0150                        ! IDIM; domain size. Linear system size is IDIM**2
+POISSON                     ! PDECOEFF: POISSON, EXP, BOX, GAUSS
+%                           ! Coefficients of the PDE
+CG                          ! Iterative method:
+%                           !    BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
 00500                       ! ITMAX
 1                           ! ITRACE
@@ -16,11 +18,12 @@ ML                        ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML
 FBGS                        ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY  r 1-level, repeats previous.
 6                           ! Number of sweeps for smoother
 1                           ! degree for polynomial smoother
-POLY_LOTTES_BETA            ! Polynomial variant
+CHEB_4_OPT            ! Polynomial variant
 %  Fields to be added for POLY
-%    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+POLY_RHO_EST_POWER     %    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
 %    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
-%    POLY_RHO_BA                   set to value
+%    POLY_RHO_BA            set to value
+1.0
 %
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
@@ -35,7 +38,11 @@ LLK                         ! AINV variant
 NONE                        ! Second (post) smoother, ignored if NONE
 6                           ! Number of sweeps for (post) smoother
 1                           ! degree for polynomial smoother
-POLY_LOTTES_BETA                 ! Polynomial variant
+CHEB_4_OPT                 ! Polynomial variant
+POLY_RHO_EST_POWER     %    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+%    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
+%    POLY_RHO_BA            set to value
+1.0
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
 NONE                        ! AS prolongation operator: NONE SUM AVG
@@ -1,8 +1,10 @@
 %%%%%%%%%%%  General  arguments % Lines starting with % are ignored.
 CSR                         ! Storage format CSR COO JAD
 0150                        ! IDIM; domain size. Linear system size is IDIM**3
-CONST                       ! PDECOEFF: CONST, EXP, BOX, GAUSS Coefficients of the PDE
-CG                          ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
+POISSON                     ! PDECOEFF: POISSON, EXP, BOX, GAUSS
+%                           ! Coefficients of the PDE
+CG                          ! Iterative method:
+%                           !    BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
 00500                       ! ITMAX
 1                           ! ITRACE
@@ -13,14 +15,15 @@ ML-VBM-VCYCLE-FBGS-D-BJAC       ! Longer descriptive name for preconditioner (up
 ML                        ! Preconditioner type: NONE JACOBI GS FBGS BJAC AS ML POLY
 %
 %%%%%%%%%%%  First smoother (for all levels but coarsest) %%%%%%%%%%%%%%%%
-FBGS                        ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY  r 1-level, repeats previous.
+POLY                        ! Smoother type JACOBI FBGS GS BWGS BJAC AS POLY  r 1-level, repeats previous.
 6                           ! Number of sweeps for smoother
-1                           ! degree for polynomial smoother
-POLY_LOTTES_BETA            ! Polynomial variant
+6                           ! degree for polynomial smoother
+CHEB_4_OPT            ! Polynomial variant
 %  Fields to be added for POLY
-%    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+POLY_RHO_EST_POWER     %    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
 %    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
-%    POLY_RHO_BA                   set to value
+%    POLY_RHO_BA            set to value
+1.0
 %
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
@@ -35,7 +38,11 @@ LLK                         ! AINV variant
 NONE                        ! Second (post) smoother, ignored if NONE
 6                           ! Number of sweeps for (post) smoother
 1                           ! degree for polynomial smoother
-POLY_LOTTES_BETA                 ! Polynomial variant
+CHEB_4_OPT                 ! Polynomial variant
+POLY_RHO_EST_POWER     %    POLY_RHO_ESTIMATE             Currently only POLY_RHO_EST_POWER
+%    POLY_RHO_ESTIMATE_ITERATIONS  default = 20
+%    POLY_RHO_BA            set to value
+1.0
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO
 NONE                        ! AS prolongation operator: NONE SUM AVG
@@ -3,7 +3,7 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_krylov -lamg_prec -lpsb_prec
+AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
 LINKOPT=

@@ -47,7 +47,7 @@
 program amg_cexample_1lev
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -62,7 +62,7 @@
 program amg_cexample_ml
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input

@@ -47,7 +47,7 @@
 program amg_dexample_1lev
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -62,7 +62,7 @@
 program amg_dexample_ml
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input

@@ -47,7 +47,7 @@
 program amg_sexample_1lev
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -62,7 +62,7 @@
 program amg_sexample_ml
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input

@@ -47,7 +47,7 @@
 program amg_zexample_1lev
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  implicit none
@@ -62,7 +62,7 @@
 program amg_zexample_ml
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input

@@ -3,7 +3,7 @@ AMGINCDIR=$(AMGDIR)/include
 include $(AMGINCDIR)/Make.inc.amg4psblas
 AMGMODDIR=$(AMGDIR)/modules
 AMGLIBDIR=$(AMGDIR)/lib
-AMG_LIBS=-L$(AMGLIBDIR) -lpsb_krylov -lamg_prec -lpsb_prec
+AMG_LIBS=-L$(AMGLIBDIR) -lpsb_linsolve -lamg_prec -lpsb_prec
 FINCLUDES=$(FMFLAG). $(FMFLAG)$(AMGMODDIR) $(FMFLAG)$(AMGINCDIR) $(PSBLAS_INCLUDES) $(FIFLAG).
 LINKOPT=

@@ -60,7 +60,7 @@
 program amg_dexample_1lev
  use psb_base_mod
  use amg_prec_mod
-  use psb_krylov_mod
+  use psb_linsolve_mod
  use psb_util_mod
  use data_input
  use amg_d_pde_mod
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
sfilippone	032543d625	Adjusted CBIND test program.	2024-11-11 13:06:48 +01:00
sfilippone	07149a02ad	Fixes for krylov->linsolve. Test program to be completed.	2024-11-11 09:11:07 +01:00
sfilippone	08a0c744b1	Switch from KRYLOV to LINSOLVE	2024-11-10 14:00:51 +01:00
sfilippone	60324084d8	New Richardson solver.	2024-11-08 18:25:13 +01:00
Salvatore Filippone	aec5a52c7f	Missing CBIND clean target	2024-10-16 10:56:53 +02:00
sfilippone	dfc261cf34	Modify error message in smth_bld	2024-09-12 18:45:20 +02:00
sfilippone	cab98295e2	Improve handling of pointers in hierarchy_bld	2024-08-06 10:28:02 +02:00
sfilippone	5a83c63810	Unused lines in ilu_solver_bld	2024-08-02 15:00:44 +02:00
sfilippone	2e43f55455	Update sample/advanced/pdegeb	2024-08-02 15:00:28 +02:00
sfilippone	244fcda207	Update test program input	2024-08-02 13:48:42 +02:00
sfilippone	ca6fce0765	Fix for test program	2024-08-02 13:46:41 +02:00
sfilippone	b6f92354d3	Reworked sample data generation	2024-08-02 12:39:09 +02:00
sfilippone	1b7fe6a9a7	Fix input data in samples pdegen	2024-07-30 09:49:46 +02:00
sfilippone	14ea4d9c15	Fix sample programs for polynomial smoothers	2024-07-29 17:35:01 +02:00
sfilippone	3ee333baac	Change name abgdxyz into upd_xyz	2024-07-29 16:59:59 +02:00
sfilippone	ecb41dfbbf	Update pde3d.inp	2024-07-16 13:41:51 +02:00
sfilippone	33ac3f786b	Merge latest changes from polysmooth	2024-07-16 13:12:10 +02:00
sfilippone	474c6a3634	Merge branch 'PolySmooth' into development	2024-07-05 16:50:26 +02:00
sfilippone	c9605d1b29	Mods for OpenMP	2024-06-04 13:25:12 +02:00
sfilippone	767b606bb2	Do away with -DOMP	2024-05-30 16:14:32 +02:00
sfilippone	8492c07521	Merge PSBCXXDEFINES into CXXDEFINES	2024-05-30 16:14:12 +02:00
sfilippone	17698c2725	Changes to OpenMP MathcBox version needs -DOMP	2024-05-30 14:43:26 +02:00
Salvatore Filippone	2ef4459b18	Added "COARSE_INVFILL"	2024-02-21 11:51:18 +01:00
sfilippone	c2fd0ac66d	Disable MATCHBOX with SERIAL_MPI and add error message	2023-11-26 11:43:10 +01:00
sfilippone	5387e206b1	Fixed sample programs	2023-11-24 16:07:22 +01:00
sfilippone	fc34385341	Fix matrix generation for samples	2023-11-16 17:55:42 +01:00
sfilippone	5fbdfb1436	Fix free for jac_solver	2023-11-16 17:55:24 +01:00