Fixes for L1-JAC vs Gauss-Seidel

mlprec/impl/smoother/mld_c_l1_jac_smoother_bld.f90

@@ -108,21 +108,9 @@ subroutine mld_c_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
       sm%nd_nnz_tot = sm%nd%get_nzeros()
       call psb_sum(ictxt,sm%nd_nnz_tot)
       arwsum = sm%nd%arwsum(info)
+      if (info == 0) call sm%sv%set_xtra_d(arwsum)
       call a%csclip(tmpa,info,&
            & jmax=nrow_a,rscale=.false.,cscale=.false.)      
-      call tmpa%mv_to(tmpcoo)
-      call tmpcoo%set_dupl(psb_dupl_add_)
-      nz = tmpcoo%get_nzeros()
-      call tmpcoo%reallocate(nz+n_row)
-      do i=1, n_row
-        tmpcoo%ia(nz+i)  = i
-        tmpcoo%ja(nz+i)  = i
-        tmpcoo%val(nz+i) = arwsum(i)
-      end do
-      call tmpcoo%set_nzeros(nz+n_row)
-      call tmpcoo%fix(info)
-      call tmpcoo%mv_to_fmt(tmpcsr,info)
-      call tmpa%mv_from(tmpcsr)
       call sm%sv%build(tmpa,desc_a,info,amold=amold,vmold=vmold)
     end if
   end select

mlprec/impl/smoother/mld_d_l1_jac_smoother_bld.f90

@@ -108,21 +108,9 @@ subroutine mld_d_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
       sm%nd_nnz_tot = sm%nd%get_nzeros()
       call psb_sum(ictxt,sm%nd_nnz_tot)
       arwsum = sm%nd%arwsum(info)
+      if (info == 0) call sm%sv%set_xtra_d(arwsum)
       call a%csclip(tmpa,info,&
            & jmax=nrow_a,rscale=.false.,cscale=.false.)      
-      call tmpa%mv_to(tmpcoo)
-      call tmpcoo%set_dupl(psb_dupl_add_)
-      nz = tmpcoo%get_nzeros()
-      call tmpcoo%reallocate(nz+n_row)
-      do i=1, n_row
-        tmpcoo%ia(nz+i)  = i
-        tmpcoo%ja(nz+i)  = i
-        tmpcoo%val(nz+i) = arwsum(i)
-      end do
-      call tmpcoo%set_nzeros(nz+n_row)
-      call tmpcoo%fix(info)
-      call tmpcoo%mv_to_fmt(tmpcsr,info)
-      call tmpa%mv_from(tmpcsr)
       call sm%sv%build(tmpa,desc_a,info,amold=amold,vmold=vmold)
     end if
   end select

mlprec/impl/smoother/mld_s_l1_jac_smoother_bld.f90

@@ -108,21 +108,9 @@ subroutine mld_s_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
       sm%nd_nnz_tot = sm%nd%get_nzeros()
       call psb_sum(ictxt,sm%nd_nnz_tot)
       arwsum = sm%nd%arwsum(info)
+      if (info == 0) call sm%sv%set_xtra_d(arwsum)
       call a%csclip(tmpa,info,&
            & jmax=nrow_a,rscale=.false.,cscale=.false.)      
-      call tmpa%mv_to(tmpcoo)
-      call tmpcoo%set_dupl(psb_dupl_add_)
-      nz = tmpcoo%get_nzeros()
-      call tmpcoo%reallocate(nz+n_row)
-      do i=1, n_row
-        tmpcoo%ia(nz+i)  = i
-        tmpcoo%ja(nz+i)  = i
-        tmpcoo%val(nz+i) = arwsum(i)
-      end do
-      call tmpcoo%set_nzeros(nz+n_row)
-      call tmpcoo%fix(info)
-      call tmpcoo%mv_to_fmt(tmpcsr,info)
-      call tmpa%mv_from(tmpcsr)
       call sm%sv%build(tmpa,desc_a,info,amold=amold,vmold=vmold)
     end if
   end select

mlprec/impl/smoother/mld_z_l1_jac_smoother_bld.f90

@@ -108,21 +108,9 @@ subroutine mld_z_l1_jac_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
       sm%nd_nnz_tot = sm%nd%get_nzeros()
       call psb_sum(ictxt,sm%nd_nnz_tot)
       arwsum = sm%nd%arwsum(info)
+      if (info == 0) call sm%sv%set_xtra_d(arwsum)
       call a%csclip(tmpa,info,&
            & jmax=nrow_a,rscale=.false.,cscale=.false.)      
-      call tmpa%mv_to(tmpcoo)
-      call tmpcoo%set_dupl(psb_dupl_add_)
-      nz = tmpcoo%get_nzeros()
-      call tmpcoo%reallocate(nz+n_row)
-      do i=1, n_row
-        tmpcoo%ia(nz+i)  = i
-        tmpcoo%ja(nz+i)  = i
-        tmpcoo%val(nz+i) = arwsum(i)
-      end do
-      call tmpcoo%set_nzeros(nz+n_row)
-      call tmpcoo%fix(info)
-      call tmpcoo%mv_to_fmt(tmpcsr,info)
-      call tmpa%mv_from(tmpcsr)
       call sm%sv%build(tmpa,desc_a,info,amold=amold,vmold=vmold)
     end if
   end select

mlprec/impl/solver/mld_c_bwgs_solver_bld.f90

@@ -79,7 +79,45 @@ subroutine mld_c_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     !
     call a%tril(sv%l,info,diag=-ione,jmax=nrow_a,u=sv%u)     
-
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm
+        type(psb_c_coo_sparse_mat) :: tcoo
+        !
+        ! For BWGS:  LX = L - D, UX = U + D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
   else
 
     info = psb_err_missing_override_method_

mlprec/impl/solver/mld_c_gs_solver_bld.f90

@@ -77,6 +77,45 @@ subroutine mld_c_gs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     ! 
     call a%tril(sv%l,info,diag=izero,jmax=nrow_a,u=sv%u)
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm, dp
+        type(psb_c_coo_sparse_mat) :: tcoo
+        !
+        ! For GS:  LX = L + D, UX = U - D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
 
   else
 

mlprec/impl/solver/mld_d_bwgs_solver_bld.f90

@@ -79,7 +79,45 @@ subroutine mld_d_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     !
     call a%tril(sv%l,info,diag=-ione,jmax=nrow_a,u=sv%u)     
-
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm
+        type(psb_d_coo_sparse_mat) :: tcoo
+        !
+        ! For BWGS:  LX = L - D, UX = U + D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
   else
 
     info = psb_err_missing_override_method_

mlprec/impl/solver/mld_d_gs_solver_bld.f90

@@ -77,6 +77,45 @@ subroutine mld_d_gs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     ! 
     call a%tril(sv%l,info,diag=izero,jmax=nrow_a,u=sv%u)
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm, dp
+        type(psb_d_coo_sparse_mat) :: tcoo
+        !
+        ! For GS:  LX = L + D, UX = U - D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
 
   else
 

mlprec/impl/solver/mld_s_bwgs_solver_bld.f90

@@ -79,7 +79,45 @@ subroutine mld_s_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     !
     call a%tril(sv%l,info,diag=-ione,jmax=nrow_a,u=sv%u)     
-
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm
+        type(psb_s_coo_sparse_mat) :: tcoo
+        !
+        ! For BWGS:  LX = L - D, UX = U + D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
   else
 
     info = psb_err_missing_override_method_

mlprec/impl/solver/mld_s_gs_solver_bld.f90

@@ -77,6 +77,45 @@ subroutine mld_s_gs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     ! 
     call a%tril(sv%l,info,diag=izero,jmax=nrow_a,u=sv%u)
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm, dp
+        type(psb_s_coo_sparse_mat) :: tcoo
+        !
+        ! For GS:  LX = L + D, UX = U - D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
 
   else
 

mlprec/impl/solver/mld_z_bwgs_solver_bld.f90

@@ -79,7 +79,45 @@ subroutine mld_z_bwgs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     !
     call a%tril(sv%l,info,diag=-ione,jmax=nrow_a,u=sv%u)     
-
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm
+        type(psb_z_coo_sparse_mat) :: tcoo
+        !
+        ! For BWGS:  LX = L - D, UX = U + D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
   else
 
     info = psb_err_missing_override_method_

mlprec/impl/solver/mld_z_gs_solver_bld.f90

@@ -77,6 +77,45 @@ subroutine mld_z_gs_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
     ! be handled by the outer Jacobi smoother.
     ! 
     call a%tril(sv%l,info,diag=izero,jmax=nrow_a,u=sv%u)
+    !
+    ! Is this an L1-GS solver? 
+    !
+    if (allocated(sv%xtra)) then
+      block
+        integer(psb_ipk_) :: k, nz, nrm, dp
+        type(psb_z_coo_sparse_mat) :: tcoo
+        !
+        ! For GS:  LX = L + D, UX = U - D
+        !
+        call sv%l%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%l%mv_from(tcoo)
+
+        call sv%u%mv_to(tcoo)
+        nrm = min(psb_size(sv%xtra),tcoo%get_nrows(),tcoo%get_ncols())
+        nz = tcoo%get_nzeros()
+        call tcoo%ensure_size(nz+nrm)
+        call tcoo%set_dupl(psb_dupl_add_)
+        do k=1,nrm
+          tcoo%ia(nz+k) = k
+          tcoo%ja(nz+k) = k
+          tcoo%val(nz+k) = -sv%xtra(k)
+        end do
+        call tcoo%set_nzeros(nz+nrm)
+        call tcoo%fix(info)
+        call sv%u%mv_from(tcoo)        
+      end block
+    end if
 
   else
 

mlprec/mld_c_base_solver_mod.f90

@@ -112,12 +112,14 @@ module mld_c_base_solver_mod
     procedure, nopass   :: get_id    => c_base_solver_get_id
     procedure, nopass   :: is_iterative => c_base_solver_is_iterative
     procedure, pass(sv) :: is_global => c_base_solver_is_global
+    procedure, pass(sv) :: set_xtra_d => c_base_solver_set_xtra_d
   end type mld_c_base_solver_type
 
   private :: c_base_solver_sizeof, c_base_solver_default,&
        &  c_base_solver_get_nzeros, c_base_solver_get_fmt, &
        &  c_base_solver_is_iterative, c_base_solver_get_id, &
-       &  c_base_solver_get_wrksize, c_base_solver_is_global
+       &  c_base_solver_get_wrksize, c_base_solver_is_global, &
+       &  c_base_solver_set_xtra_d
 
 
   interface  
@@ -418,4 +420,14 @@ contains
     val = 0
   end function c_base_solver_get_wrksize
 
+  subroutine c_base_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_c_base_solver_type), intent(inout) :: sv
+    real(psb_spk_), intent(in)  :: d(:)
+    ! Do nothing for base version
+
+    return
+  end subroutine c_base_solver_set_xtra_d
+
 end module mld_c_base_solver_mod

mlprec/mld_c_gs_solver.f90

@@ -59,6 +59,7 @@ module mld_c_gs_solver
     type(psb_cspmat_type)      :: l, u
     integer(psb_ipk_)          :: sweeps
     real(psb_spk_)             :: eps
+    real(psb_spk_), allocatable  :: xtra(:)
   contains
     procedure, pass(sv) :: dump    => mld_c_gs_solver_dmp
     procedure, pass(sv) :: check   => c_gs_solver_check
@@ -77,6 +78,7 @@ module mld_c_gs_solver
     procedure, pass(sv) :: default => c_gs_solver_default
     procedure, pass(sv) :: sizeof  => c_gs_solver_sizeof
     procedure, pass(sv) :: get_nzeros => c_gs_solver_get_nzeros
+    procedure, pass(sv) :: set_xtra_d => c_gs_solver_set_xtra_d    
     procedure, nopass   :: get_wrksz => c_gs_solver_get_wrksize
     procedure, nopass   :: get_fmt    => c_gs_solver_get_fmt
     procedure, nopass   :: get_id    => c_gs_solver_get_id
@@ -585,4 +587,16 @@ contains
     val = 2
   end function c_gs_solver_get_wrksize
 
+
+  subroutine c_gs_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_c_gs_solver_type), intent(inout) :: sv
+    real(psb_spk_), intent(in)  :: d(:)
+
+    sv%xtra = d
+    
+    return
+  end subroutine c_gs_solver_set_xtra_d
+
 end module mld_c_gs_solver

mlprec/mld_d_base_solver_mod.f90

@@ -112,12 +112,14 @@ module mld_d_base_solver_mod
     procedure, nopass   :: get_id    => d_base_solver_get_id
     procedure, nopass   :: is_iterative => d_base_solver_is_iterative
     procedure, pass(sv) :: is_global => d_base_solver_is_global
+    procedure, pass(sv) :: set_xtra_d => d_base_solver_set_xtra_d
   end type mld_d_base_solver_type
 
   private :: d_base_solver_sizeof, d_base_solver_default,&
        &  d_base_solver_get_nzeros, d_base_solver_get_fmt, &
        &  d_base_solver_is_iterative, d_base_solver_get_id, &
-       &  d_base_solver_get_wrksize, d_base_solver_is_global
+       &  d_base_solver_get_wrksize, d_base_solver_is_global, &
+       &  d_base_solver_set_xtra_d
 
 
   interface  
@@ -418,4 +420,14 @@ contains
     val = 0
   end function d_base_solver_get_wrksize
 
+  subroutine d_base_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_d_base_solver_type), intent(inout) :: sv
+    real(psb_dpk_), intent(in)  :: d(:)
+    ! Do nothing for base version
+
+    return
+  end subroutine d_base_solver_set_xtra_d
+
 end module mld_d_base_solver_mod

mlprec/mld_d_gs_solver.f90

@@ -59,6 +59,7 @@ module mld_d_gs_solver
     type(psb_dspmat_type)      :: l, u
     integer(psb_ipk_)          :: sweeps
     real(psb_dpk_)             :: eps
+    real(psb_dpk_), allocatable  :: xtra(:)
   contains
     procedure, pass(sv) :: dump    => mld_d_gs_solver_dmp
     procedure, pass(sv) :: check   => d_gs_solver_check
@@ -77,6 +78,7 @@ module mld_d_gs_solver
     procedure, pass(sv) :: default => d_gs_solver_default
     procedure, pass(sv) :: sizeof  => d_gs_solver_sizeof
     procedure, pass(sv) :: get_nzeros => d_gs_solver_get_nzeros
+    procedure, pass(sv) :: set_xtra_d => d_gs_solver_set_xtra_d    
     procedure, nopass   :: get_wrksz => d_gs_solver_get_wrksize
     procedure, nopass   :: get_fmt    => d_gs_solver_get_fmt
     procedure, nopass   :: get_id    => d_gs_solver_get_id
@@ -585,4 +587,16 @@ contains
     val = 2
   end function d_gs_solver_get_wrksize
 
+
+  subroutine d_gs_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_d_gs_solver_type), intent(inout) :: sv
+    real(psb_dpk_), intent(in)  :: d(:)
+
+    sv%xtra = d
+    
+    return
+  end subroutine d_gs_solver_set_xtra_d
+
 end module mld_d_gs_solver

mlprec/mld_s_base_solver_mod.f90

@@ -112,12 +112,14 @@ module mld_s_base_solver_mod
     procedure, nopass   :: get_id    => s_base_solver_get_id
     procedure, nopass   :: is_iterative => s_base_solver_is_iterative
     procedure, pass(sv) :: is_global => s_base_solver_is_global
+    procedure, pass(sv) :: set_xtra_d => s_base_solver_set_xtra_d
   end type mld_s_base_solver_type
 
   private :: s_base_solver_sizeof, s_base_solver_default,&
        &  s_base_solver_get_nzeros, s_base_solver_get_fmt, &
        &  s_base_solver_is_iterative, s_base_solver_get_id, &
-       &  s_base_solver_get_wrksize, s_base_solver_is_global
+       &  s_base_solver_get_wrksize, s_base_solver_is_global, &
+       &  s_base_solver_set_xtra_d
 
 
   interface  
@@ -418,4 +420,14 @@ contains
     val = 0
   end function s_base_solver_get_wrksize
 
+  subroutine s_base_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_s_base_solver_type), intent(inout) :: sv
+    real(psb_spk_), intent(in)  :: d(:)
+    ! Do nothing for base version
+
+    return
+  end subroutine s_base_solver_set_xtra_d
+
 end module mld_s_base_solver_mod

mlprec/mld_s_gs_solver.f90

@@ -59,6 +59,7 @@ module mld_s_gs_solver
     type(psb_sspmat_type)      :: l, u
     integer(psb_ipk_)          :: sweeps
     real(psb_spk_)             :: eps
+    real(psb_spk_), allocatable  :: xtra(:)
   contains
     procedure, pass(sv) :: dump    => mld_s_gs_solver_dmp
     procedure, pass(sv) :: check   => s_gs_solver_check
@@ -77,6 +78,7 @@ module mld_s_gs_solver
     procedure, pass(sv) :: default => s_gs_solver_default
     procedure, pass(sv) :: sizeof  => s_gs_solver_sizeof
     procedure, pass(sv) :: get_nzeros => s_gs_solver_get_nzeros
+    procedure, pass(sv) :: set_xtra_d => s_gs_solver_set_xtra_d    
     procedure, nopass   :: get_wrksz => s_gs_solver_get_wrksize
     procedure, nopass   :: get_fmt    => s_gs_solver_get_fmt
     procedure, nopass   :: get_id    => s_gs_solver_get_id
@@ -585,4 +587,16 @@ contains
     val = 2
   end function s_gs_solver_get_wrksize
 
+
+  subroutine s_gs_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_s_gs_solver_type), intent(inout) :: sv
+    real(psb_spk_), intent(in)  :: d(:)
+
+    sv%xtra = d
+    
+    return
+  end subroutine s_gs_solver_set_xtra_d
+
 end module mld_s_gs_solver

mlprec/mld_z_base_solver_mod.f90

@@ -112,12 +112,14 @@ module mld_z_base_solver_mod
     procedure, nopass   :: get_id    => z_base_solver_get_id
     procedure, nopass   :: is_iterative => z_base_solver_is_iterative
     procedure, pass(sv) :: is_global => z_base_solver_is_global
+    procedure, pass(sv) :: set_xtra_d => z_base_solver_set_xtra_d
   end type mld_z_base_solver_type
 
   private :: z_base_solver_sizeof, z_base_solver_default,&
        &  z_base_solver_get_nzeros, z_base_solver_get_fmt, &
        &  z_base_solver_is_iterative, z_base_solver_get_id, &
-       &  z_base_solver_get_wrksize, z_base_solver_is_global
+       &  z_base_solver_get_wrksize, z_base_solver_is_global, &
+       &  z_base_solver_set_xtra_d
 
 
   interface  
@@ -418,4 +420,14 @@ contains
     val = 0
   end function z_base_solver_get_wrksize
 
+  subroutine z_base_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_z_base_solver_type), intent(inout) :: sv
+    real(psb_dpk_), intent(in)  :: d(:)
+    ! Do nothing for base version
+
+    return
+  end subroutine z_base_solver_set_xtra_d
+
 end module mld_z_base_solver_mod

mlprec/mld_z_gs_solver.f90

@@ -59,6 +59,7 @@ module mld_z_gs_solver
     type(psb_zspmat_type)      :: l, u
     integer(psb_ipk_)          :: sweeps
     real(psb_dpk_)             :: eps
+    real(psb_dpk_), allocatable  :: xtra(:)
   contains
     procedure, pass(sv) :: dump    => mld_z_gs_solver_dmp
     procedure, pass(sv) :: check   => z_gs_solver_check
@@ -77,6 +78,7 @@ module mld_z_gs_solver
     procedure, pass(sv) :: default => z_gs_solver_default
     procedure, pass(sv) :: sizeof  => z_gs_solver_sizeof
     procedure, pass(sv) :: get_nzeros => z_gs_solver_get_nzeros
+    procedure, pass(sv) :: set_xtra_d => z_gs_solver_set_xtra_d    
     procedure, nopass   :: get_wrksz => z_gs_solver_get_wrksize
     procedure, nopass   :: get_fmt    => z_gs_solver_get_fmt
     procedure, nopass   :: get_id    => z_gs_solver_get_id
@@ -585,4 +587,16 @@ contains
     val = 2
   end function z_gs_solver_get_wrksize
 
+
+  subroutine z_gs_solver_set_xtra_d(sv,d) 
+    implicit none 
+    ! Arguments
+    class(mld_z_gs_solver_type), intent(inout) :: sv
+    real(psb_dpk_), intent(in)  :: d(:)
+
+    sv%xtra = d
+    
+    return
+  end subroutine z_gs_solver_set_xtra_d
+
 end module mld_z_gs_solver