mld2p4-2:

mlprec/impl/mld_cmlprec_aply.f90
 mlprec/impl/mld_dmlprec_aply.f90
 mlprec/impl/mld_smlprec_aply.f90
 mlprec/impl/mld_zmlprec_aply.f90
 mlprec/impl/smoother/mld_c_as_smoother_apply.f90
 mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_c_base_smoother_apply.f90
 mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_c_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_d_as_smoother_apply.f90
 mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_d_as_smoother_bld.f90
 mlprec/impl/smoother/mld_d_as_smoother_cnv.f90
 mlprec/impl/smoother/mld_d_base_smoother_apply.f90
 mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_d_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_s_as_smoother_apply.f90
 mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_s_base_smoother_apply.f90
 mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_s_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_z_as_smoother_apply.f90
 mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_z_base_smoother_apply.f90
 mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90
 mlprec/impl/smoother/mld_z_jac_smoother_apply.f90
 mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90
 mlprec/impl/solver/mld_c_base_solver_apply.f90
 mlprec/impl/solver/mld_c_base_solver_apply_vect.f90
 mlprec/impl/solver/mld_c_bwgs_solver_apply.f90
 mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90
 mlprec/impl/solver/mld_c_diag_solver_apply.f90
 mlprec/impl/solver/mld_c_diag_solver_apply_vect.f90
 mlprec/impl/solver/mld_c_gs_solver_apply.f90
 mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90
 mlprec/impl/solver/mld_c_id_solver_apply.f90
 mlprec/impl/solver/mld_c_id_solver_apply_vect.f90
 mlprec/impl/solver/mld_c_ilu_solver_apply.f90
 mlprec/impl/solver/mld_c_ilu_solver_apply_vect.f90
 mlprec/impl/solver/mld_c_mumps_solver_apply.F90
 mlprec/impl/solver/mld_c_mumps_solver_apply_vect.F90
 mlprec/impl/solver/mld_d_base_solver_apply.f90
 mlprec/impl/solver/mld_d_base_solver_apply_vect.f90
 mlprec/impl/solver/mld_d_bwgs_solver_apply.f90
 mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90
 mlprec/impl/solver/mld_d_diag_solver_apply.f90
 mlprec/impl/solver/mld_d_diag_solver_apply_vect.f90
 mlprec/impl/solver/mld_d_gs_solver_apply.f90
 mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90
 mlprec/impl/solver/mld_d_id_solver_apply.f90
 mlprec/impl/solver/mld_d_id_solver_apply_vect.f90
 mlprec/impl/solver/mld_d_ilu_solver_apply.f90
 mlprec/impl/solver/mld_d_ilu_solver_apply_vect.f90
 mlprec/impl/solver/mld_d_mumps_solver_apply.F90
 mlprec/impl/solver/mld_d_mumps_solver_apply_vect.F90
 mlprec/impl/solver/mld_s_base_solver_apply.f90
 mlprec/impl/solver/mld_s_base_solver_apply_vect.f90
 mlprec/impl/solver/mld_s_bwgs_solver_apply.f90
 mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90
 mlprec/impl/solver/mld_s_diag_solver_apply.f90
 mlprec/impl/solver/mld_s_diag_solver_apply_vect.f90
 mlprec/impl/solver/mld_s_gs_solver_apply.f90
 mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90
 mlprec/impl/solver/mld_s_id_solver_apply.f90
 mlprec/impl/solver/mld_s_id_solver_apply_vect.f90
 mlprec/impl/solver/mld_s_ilu_solver_apply.f90
 mlprec/impl/solver/mld_s_ilu_solver_apply_vect.f90
 mlprec/impl/solver/mld_s_mumps_solver_apply.F90
 mlprec/impl/solver/mld_s_mumps_solver_apply_vect.F90
 mlprec/impl/solver/mld_z_base_solver_apply.f90
 mlprec/impl/solver/mld_z_base_solver_apply_vect.f90
 mlprec/impl/solver/mld_z_bwgs_solver_apply.f90
 mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90
 mlprec/impl/solver/mld_z_diag_solver_apply.f90
 mlprec/impl/solver/mld_z_diag_solver_apply_vect.f90
 mlprec/impl/solver/mld_z_gs_solver_apply.f90
 mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90
 mlprec/impl/solver/mld_z_id_solver_apply.f90
 mlprec/impl/solver/mld_z_id_solver_apply_vect.f90
 mlprec/impl/solver/mld_z_ilu_solver_apply.f90
 mlprec/impl/solver/mld_z_ilu_solver_apply_vect.f90
 mlprec/impl/solver/mld_z_mumps_solver_apply.F90
 mlprec/impl/solver/mld_z_mumps_solver_apply_vect.F90
 mlprec/mld_base_prec_type.F90
 mlprec/mld_c_as_smoother.f90
 mlprec/mld_c_base_smoother_mod.f90
 mlprec/mld_c_base_solver_mod.f90
 mlprec/mld_c_diag_solver.f90
 mlprec/mld_c_gs_solver.f90
 mlprec/mld_c_id_solver.f90
 mlprec/mld_c_ilu_solver.f90
 mlprec/mld_c_jac_smoother.f90
 mlprec/mld_c_mumps_solver.F90
 mlprec/mld_c_slu_solver.F90
 mlprec/mld_c_sludist_solver.F90
 mlprec/mld_c_umf_solver.F90
 mlprec/mld_d_as_smoother.f90
 mlprec/mld_d_base_smoother_mod.f90
 mlprec/mld_d_base_solver_mod.f90
 mlprec/mld_d_diag_solver.f90
 mlprec/mld_d_gs_solver.f90
 mlprec/mld_d_id_solver.f90
 mlprec/mld_d_ilu_solver.f90
 mlprec/mld_d_jac_smoother.f90
 mlprec/mld_d_mumps_solver.F90
 mlprec/mld_d_slu_solver.F90
 mlprec/mld_d_sludist_solver.F90
 mlprec/mld_d_umf_solver.F90
 mlprec/mld_s_as_smoother.f90
 mlprec/mld_s_base_smoother_mod.f90
 mlprec/mld_s_base_solver_mod.f90
 mlprec/mld_s_diag_solver.f90
 mlprec/mld_s_gs_solver.f90
 mlprec/mld_s_id_solver.f90
 mlprec/mld_s_ilu_solver.f90
 mlprec/mld_s_jac_smoother.f90
 mlprec/mld_s_mumps_solver.F90
 mlprec/mld_s_slu_solver.F90
 mlprec/mld_s_sludist_solver.F90
 mlprec/mld_s_umf_solver.F90
 mlprec/mld_z_as_smoother.f90
 mlprec/mld_z_base_smoother_mod.f90
 mlprec/mld_z_base_solver_mod.f90
 mlprec/mld_z_diag_solver.f90
 mlprec/mld_z_gs_solver.f90
 mlprec/mld_z_id_solver.f90
 mlprec/mld_z_ilu_solver.f90
 mlprec/mld_z_jac_smoother.f90
 mlprec/mld_z_mumps_solver.F90
 mlprec/mld_z_slu_solver.F90
 mlprec/mld_z_sludist_solver.F90
 mlprec/mld_z_umf_solver.F90

Fixed parms print for V/Wcycle. 
Reworked MLPREC_APLY; now the multilevel code is more flexible and
readable.
Introduced option for initial vector in smoothers and solvers;
non-iterative solvers (currently all except GS) ignore it.

mlprec/impl/smoother/mld_c_as_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
-  
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_apply
   implicit none 
@@ -50,11 +50,13 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   integer(psb_ipk_), intent(in)         :: sweeps
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, nrow_d, i
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='c_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -63,6 +65,12 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -84,6 +92,7 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   n_col = sm%desc_data%get_local_cols()
   nrow_d = desc_data%get_local_rows()
   isz    = max(n_row,N_COL)
+
   if ((6*isz) <= size(work)) then 
     ww => work(1:isz)
     tx => work(isz+1:2*isz)
@@ -129,10 +138,10 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else if ((sm%novr == 0).and.(sweeps == 1)) then 
+  else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -226,7 +235,7 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
         goto 9999
       end select
 
-      call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(cone,tx,czero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -299,7 +308,24 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
       ! to compute an approximate solution of a linear system.
       !
       !
+      select case (init_)
+      case('Z') 
         ty = czero
+      case('Y')
+        call psb_geaxpby(cone,y,czero,ty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(cone,initu,czero,ty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
+      
       do i=1, sweeps
         select case(trans_)
         case('N')
@@ -386,7 +412,7 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_apply_vect
   implicit none 
@@ -50,13 +50,15 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                  :: sweeps
   complex(psb_spk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col, nrow_d, i
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: vx(:) 
   type(psb_c_vect_type) :: vtx, vty, vww
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='c_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -65,6 +67,12 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -135,7 +143,7 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -240,7 +248,7 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         goto 9999
       end select
 
-      call sm%sv%apply(cone,vtx,czero,vty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(cone,vtx,czero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -313,7 +321,23 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
       ! to compute an approximate solution of a linear system.
       !
       !
-      call vty%set(czero)
+      select case (init_)
+      case('Z') 
+        call vty%zero()
+      case('Y')
+        call psb_geaxpby(cone,y,czero,vty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(cone,initu,czero,vty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
       
       do i=1, sweeps
         select case(trans_)
@@ -401,7 +425,7 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(cone,vww,czero,vty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(cone,vww,czero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_c_base_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_c_base_smoother_mod, mld_protect_name =>  mld_c_base_smoother_apply
   implicit none 
@@ -50,6 +50,8 @@ subroutine mld_c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
   integer(psb_ipk_), intent(in)                :: sweeps
   complex(psb_spk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='c_base_smoother_apply'
@@ -66,7 +68,7 @@ subroutine mld_c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
 
   else
     if (allocated(sm%sv)) then 
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90

@@ -37,8 +37,7 @@
 !!$ 
 !!$
 subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-  
+     &  trans,sweeps,work,info,init,initu)
-     &  trans,sweeps,work,info)
   use psb_base_mod
   use mld_c_base_smoother_mod, mld_protect_name =>  mld_c_base_smoother_apply_vect
   implicit none 
@@ -51,6 +50,8 @@ subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
   integer(psb_ipk_), intent(in)                    :: sweeps
   complex(psb_spk_),target, intent(inout)            :: work(:)
   integer(psb_ipk_), intent(out)                   :: info
+  character, intent(in), optional                  :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='c_base_smoother_apply'
@@ -67,7 +68,7 @@ subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
 
   else
     if (allocated(sm%sv)) then
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_c_jac_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply
   implicit none 
@@ -50,17 +50,27 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
   integer(psb_ipk_), intent(in)         :: sweeps
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)    :: n_row,n_col
-  complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
+  complex(psb_spk_), allocatable  :: tx(:),ty(:)
+  complex(psb_spk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='c_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -112,8 +122,8 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else  if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then 
+  else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -122,38 +132,47 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
       goto 9999
     endif
 
-  else if (sweeps  > 1) then 
+  else if (sweeps >= 1) then 
-
     !
     !
     ! Apply multiple sweeps of a block-Jacobi solver
     ! to compute an approximate solution of a linear system.
     !
     !
-    allocate(tx(n_col),ty(n_col),stat=info)
+    
-    if (info /= psb_success_) then 
+    call psb_geasb(tx,desc_data,info)
-      info=psb_err_alloc_request_
+    call psb_geasb(ty,desc_data,info)
-      call psb_errpush(info,name,& 
+
-           & i_err=(/2*n_col,izero,izero,izero,izero/),&
+    select case (init_)
-           & a_err='complex(psb_spk_)')
+    case('Z') 
+      tx(:) = czero
+    case('Y')
+      call psb_geaxpby(cone,y,czero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
         goto 9999
       end if
+      call psb_geaxpby(cone,initu,czero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
-    tx = czero
-    ty = czero
     do i=1, sweeps
       !
       ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
       ! block diagonal part and the remaining part of the local matrix
       ! and Y(j) is the approximate solution at sweep j.
       !
-      ty(1:n_row) = x(1:n_row)
+      call psb_geaxpby(cone,x,czero,ty,desc_data,info)
-      call psb_spmm(-cone,sm%nd,tx,cone,ty,desc_data,info,&
+      call psb_spmm(-cone,sm%nd,tx,cone,ty,desc_data,info,work=aux,trans=trans_)
-           & work=aux,trans=trans_)
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(cone,ty,czero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(cone,ty,czero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply_vect
@@ -51,23 +51,31 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                   :: sweeps
   complex(psb_spk_),target, intent(inout)           :: work(:)
   integer(psb_ipk_), intent(out)                  :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col
   type(psb_c_vect_type)  :: tx, ty
   complex(psb_spk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='c_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
-  case('T')
+  case('T','C')
-  case('C')
   case default
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
@@ -107,8 +115,6 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     end if
   endif
 
-!!$  write(0,*) 'Jacobi   smoother with ',sweeps
-  
   if (sweeps == 0) then
     
     !
@@ -118,7 +124,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
   else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -134,9 +140,27 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     ! to compute an approximate solution of a linear system.
     !
     !
-    call tx%bld(x%get_nrows(),mold=x%v)
+    
-    call tx%set(czero)
+    call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
-    call ty%bld(x%get_nrows(),mold=x%v)
+    call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
+
+    select case (init_)
+    case('Z') 
+      call tx%zero()
+    case('Y')
+      call psb_geaxpby(cone,y,czero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
+        goto 9999
+      end if
+      call psb_geaxpby(cone,initu,czero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
     do i=1, sweeps
       !
@@ -149,7 +173,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(cone,ty,czero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(cone,ty,czero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_d_as_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
-  
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_apply
   implicit none 
@@ -50,11 +50,13 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   integer(psb_ipk_), intent(in)         :: sweeps
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, nrow_d, i
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -63,6 +65,12 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -84,6 +92,7 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   n_col = sm%desc_data%get_local_cols()
   nrow_d = desc_data%get_local_rows()
   isz    = max(n_row,N_COL)
+
   if ((6*isz) <= size(work)) then 
     ww => work(1:isz)
     tx => work(isz+1:2*isz)
@@ -129,10 +138,10 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else if ((sm%novr == 0).and.(sweeps == 1)) then 
+  else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -226,7 +235,7 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
         goto 9999
       end select
 
-      call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -299,7 +308,24 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
       ! to compute an approximate solution of a linear system.
       !
       !
+      select case (init_)
+      case('Z') 
         ty = dzero
+      case('Y')
+        call psb_geaxpby(done,y,dzero,ty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(done,initu,dzero,ty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
+      
       do i=1, sweeps
         select case(trans_)
         case('N')
@@ -386,7 +412,7 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_apply_vect
   implicit none 
@@ -50,13 +50,15 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                  :: sweeps
   real(psb_dpk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col, nrow_d, i
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: vx(:) 
   type(psb_d_vect_type) :: vtx, vty, vww
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -65,6 +67,12 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -135,7 +143,7 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -240,7 +248,7 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         goto 9999
       end select
 
-      call sm%sv%apply(done,vtx,dzero,vty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(done,vtx,dzero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -313,7 +321,23 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
       ! to compute an approximate solution of a linear system.
       !
       !
-      call vty%set(dzero)
+      select case (init_)
+      case('Z') 
+        call vty%zero()
+      case('Y')
+        call psb_geaxpby(done,y,dzero,vty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(done,initu,dzero,vty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
       
       do i=1, sweeps
         select case(trans_)
@@ -401,7 +425,7 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(done,vww,dzero,vty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(done,vww,dzero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_d_as_smoother_bld.f90

@@ -142,12 +142,14 @@ subroutine mld_d_as_smoother_bld(a,desc_a,sm,upd,info,amold,vmold,imold)
          & blck%get_nrows(), blck%get_nzeros()
 
   End if
-  if (info == psb_success_) call sm%sv%build(a,sm%desc_data,upd,info,&
+  if (info == psb_success_) &
+       & call sm%sv%build(a,sm%desc_data,upd,info,&
        &  blck,amold=amold,vmold=vmold)
 
   nrow_a = a%get_nrows()
   n_row  = sm%desc_data%get_local_rows()
   n_col  = sm%desc_data%get_local_cols()
+
   if (info == psb_success_) call a%csclip(sm%nd,info,&
        & jmin=nrow_a+1,rscale=.false.,cscale=.false.)
   if (info == psb_success_) call blck%csclip(atmp,info,&

mlprec/impl/smoother/mld_d_base_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_d_base_smoother_mod, mld_protect_name =>  mld_d_base_smoother_apply
   implicit none 
@@ -50,6 +50,8 @@ subroutine mld_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
   integer(psb_ipk_), intent(in)                :: sweeps
   real(psb_dpk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='d_base_smoother_apply'
@@ -66,7 +68,7 @@ subroutine mld_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
 
   else
     if (allocated(sm%sv)) then 
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90

@@ -37,8 +37,7 @@
 !!$ 
 !!$
 subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-  
+     &  trans,sweeps,work,info,init,initu)
-     &  trans,sweeps,work,info)
   use psb_base_mod
   use mld_d_base_smoother_mod, mld_protect_name =>  mld_d_base_smoother_apply_vect
   implicit none 
@@ -51,6 +50,8 @@ subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
   integer(psb_ipk_), intent(in)                    :: sweeps
   real(psb_dpk_),target, intent(inout)            :: work(:)
   integer(psb_ipk_), intent(out)                   :: info
+  character, intent(in), optional                  :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='d_base_smoother_apply'
@@ -67,7 +68,7 @@ subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
 
   else
     if (allocated(sm%sv)) then
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_d_jac_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply
   implicit none 
@@ -50,17 +50,27 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
   integer(psb_ipk_), intent(in)         :: sweeps
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)    :: n_row,n_col
-  real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
+  real(psb_dpk_), allocatable  :: tx(:),ty(:)
+  real(psb_dpk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -112,8 +122,8 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else  if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then 
+  else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -122,38 +132,47 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
       goto 9999
     endif
 
-  else if (sweeps  > 1) then 
+  else if (sweeps >= 1) then 
-
     !
     !
     ! Apply multiple sweeps of a block-Jacobi solver
     ! to compute an approximate solution of a linear system.
     !
     !
-    allocate(tx(n_col),ty(n_col),stat=info)
+    
-    if (info /= psb_success_) then 
+    call psb_geasb(tx,desc_data,info)
-      info=psb_err_alloc_request_
+    call psb_geasb(ty,desc_data,info)
-      call psb_errpush(info,name,& 
+
-           & i_err=(/2*n_col,izero,izero,izero,izero/),&
+    select case (init_)
-           & a_err='real(psb_dpk_)')
+    case('Z') 
+      tx(:) = dzero
+    case('Y')
+      call psb_geaxpby(done,y,dzero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
         goto 9999
       end if
+      call psb_geaxpby(done,initu,dzero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
-    tx = dzero
-    ty = dzero
     do i=1, sweeps
       !
       ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
       ! block diagonal part and the remaining part of the local matrix
       ! and Y(j) is the approximate solution at sweep j.
       !
-      ty(1:n_row) = x(1:n_row)
+      call psb_geaxpby(done,x,dzero,ty,desc_data,info)
-      call psb_spmm(-done,sm%nd,tx,done,ty,desc_data,info,&
+      call psb_spmm(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_)
-           & work=aux,trans=trans_)
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(done,ty,dzero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(done,ty,dzero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply_vect
@@ -51,23 +51,31 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                   :: sweeps
   real(psb_dpk_),target, intent(inout)           :: work(:)
   integer(psb_ipk_), intent(out)                  :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col
   type(psb_d_vect_type)  :: tx, ty
   real(psb_dpk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
-  case('T')
+  case('T','C')
-  case('C')
   case default
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
@@ -107,8 +115,6 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     end if
   endif
 
-!!$  write(0,*) 'Jacobi   smoother with ',sweeps
-  
   if (sweeps == 0) then
     
     !
@@ -118,7 +124,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
   else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -134,9 +140,27 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     ! to compute an approximate solution of a linear system.
     !
     !
-    call tx%bld(x%get_nrows(),mold=x%v)
+    
-    call tx%set(dzero)
+    call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
-    call ty%bld(x%get_nrows(),mold=x%v)
+    call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
+
+    select case (init_)
+    case('Z') 
+      call tx%zero()
+    case('Y')
+      call psb_geaxpby(done,y,dzero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
+        goto 9999
+      end if
+      call psb_geaxpby(done,initu,dzero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
     do i=1, sweeps
       !
@@ -149,7 +173,7 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(done,ty,dzero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(done,ty,dzero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_s_as_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
-  
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_apply
   implicit none 
@@ -50,11 +50,13 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   integer(psb_ipk_), intent(in)         :: sweeps
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, nrow_d, i
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='s_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -63,6 +65,12 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -84,6 +92,7 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   n_col = sm%desc_data%get_local_cols()
   nrow_d = desc_data%get_local_rows()
   isz    = max(n_row,N_COL)
+
   if ((6*isz) <= size(work)) then 
     ww => work(1:isz)
     tx => work(isz+1:2*isz)
@@ -129,10 +138,10 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else if ((sm%novr == 0).and.(sweeps == 1)) then 
+  else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -226,7 +235,7 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
         goto 9999
       end select
 
-      call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(sone,tx,szero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -299,7 +308,24 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
       ! to compute an approximate solution of a linear system.
       !
       !
+      select case (init_)
+      case('Z') 
         ty = szero
+      case('Y')
+        call psb_geaxpby(sone,y,szero,ty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(sone,initu,szero,ty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
+      
       do i=1, sweeps
         select case(trans_)
         case('N')
@@ -386,7 +412,7 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_apply_vect
   implicit none 
@@ -50,13 +50,15 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                  :: sweeps
   real(psb_spk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col, nrow_d, i
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: vx(:) 
   type(psb_s_vect_type) :: vtx, vty, vww
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='s_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -65,6 +67,12 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -135,7 +143,7 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -240,7 +248,7 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         goto 9999
       end select
 
-      call sm%sv%apply(sone,vtx,szero,vty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(sone,vtx,szero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -313,7 +321,23 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
       ! to compute an approximate solution of a linear system.
       !
       !
-      call vty%set(szero)
+      select case (init_)
+      case('Z') 
+        call vty%zero()
+      case('Y')
+        call psb_geaxpby(sone,y,szero,vty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(sone,initu,szero,vty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
       
       do i=1, sweeps
         select case(trans_)
@@ -401,7 +425,7 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(sone,vww,szero,vty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(sone,vww,szero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_s_base_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_s_base_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_s_base_smoother_mod, mld_protect_name =>  mld_s_base_smoother_apply
   implicit none 
@@ -50,6 +50,8 @@ subroutine mld_s_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
   integer(psb_ipk_), intent(in)                :: sweeps
   real(psb_spk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='s_base_smoother_apply'
@@ -66,7 +68,7 @@ subroutine mld_s_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
 
   else
     if (allocated(sm%sv)) then 
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90

@@ -37,8 +37,7 @@
 !!$ 
 !!$
 subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-  
+     &  trans,sweeps,work,info,init,initu)
-     &  trans,sweeps,work,info)
   use psb_base_mod
   use mld_s_base_smoother_mod, mld_protect_name =>  mld_s_base_smoother_apply_vect
   implicit none 
@@ -51,6 +50,8 @@ subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
   integer(psb_ipk_), intent(in)                    :: sweeps
   real(psb_spk_),target, intent(inout)            :: work(:)
   integer(psb_ipk_), intent(out)                   :: info
+  character, intent(in), optional                  :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='s_base_smoother_apply'
@@ -67,7 +68,7 @@ subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
 
   else
     if (allocated(sm%sv)) then
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_s_jac_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply
   implicit none 
@@ -50,17 +50,27 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
   integer(psb_ipk_), intent(in)         :: sweeps
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)    :: n_row,n_col
-  real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
+  real(psb_spk_), allocatable  :: tx(:),ty(:)
+  real(psb_spk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='s_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -112,8 +122,8 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else  if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then 
+  else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -122,38 +132,47 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
       goto 9999
     endif
 
-  else if (sweeps  > 1) then 
+  else if (sweeps >= 1) then 
-
     !
     !
     ! Apply multiple sweeps of a block-Jacobi solver
     ! to compute an approximate solution of a linear system.
     !
     !
-    allocate(tx(n_col),ty(n_col),stat=info)
+    
-    if (info /= psb_success_) then 
+    call psb_geasb(tx,desc_data,info)
-      info=psb_err_alloc_request_
+    call psb_geasb(ty,desc_data,info)
-      call psb_errpush(info,name,& 
+
-           & i_err=(/2*n_col,izero,izero,izero,izero/),&
+    select case (init_)
-           & a_err='real(psb_spk_)')
+    case('Z') 
+      tx(:) = szero
+    case('Y')
+      call psb_geaxpby(sone,y,szero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
         goto 9999
       end if
+      call psb_geaxpby(sone,initu,szero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
-    tx = szero
-    ty = szero
     do i=1, sweeps
       !
       ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
       ! block diagonal part and the remaining part of the local matrix
       ! and Y(j) is the approximate solution at sweep j.
       !
-      ty(1:n_row) = x(1:n_row)
+      call psb_geaxpby(sone,x,szero,ty,desc_data,info)
-      call psb_spmm(-sone,sm%nd,tx,sone,ty,desc_data,info,&
+      call psb_spmm(-sone,sm%nd,tx,sone,ty,desc_data,info,work=aux,trans=trans_)
-           & work=aux,trans=trans_)
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(sone,ty,szero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(sone,ty,szero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply_vect
@@ -51,23 +51,31 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                   :: sweeps
   real(psb_spk_),target, intent(inout)           :: work(:)
   integer(psb_ipk_), intent(out)                  :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col
   type(psb_s_vect_type)  :: tx, ty
   real(psb_spk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='s_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
-  case('T')
+  case('T','C')
-  case('C')
   case default
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
@@ -107,8 +115,6 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     end if
   endif
 
-!!$  write(0,*) 'Jacobi   smoother with ',sweeps
-  
   if (sweeps == 0) then
     
     !
@@ -118,7 +124,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
   else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -134,9 +140,27 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     ! to compute an approximate solution of a linear system.
     !
     !
-    call tx%bld(x%get_nrows(),mold=x%v)
+    
-    call tx%set(szero)
+    call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
-    call ty%bld(x%get_nrows(),mold=x%v)
+    call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
+
+    select case (init_)
+    case('Z') 
+      call tx%zero()
+    case('Y')
+      call psb_geaxpby(sone,y,szero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
+        goto 9999
+      end if
+      call psb_geaxpby(sone,initu,szero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
     do i=1, sweeps
       !
@@ -149,7 +173,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(sone,ty,szero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(sone,ty,szero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_z_as_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,&
-  
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_apply
   implicit none 
@@ -50,11 +50,13 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   integer(psb_ipk_), intent(in)         :: sweeps
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, nrow_d, i
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='z_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -63,6 +65,12 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -84,6 +92,7 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
   n_col = sm%desc_data%get_local_cols()
   nrow_d = desc_data%get_local_rows()
   isz    = max(n_row,N_COL)
+
   if ((6*isz) <= size(work)) then 
     ww => work(1:isz)
     tx => work(isz+1:2*isz)
@@ -129,10 +138,10 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else if ((sm%novr == 0).and.(sweeps == 1)) then 
+  else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -226,7 +235,7 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
         goto 9999
       end select
 
-      call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(zone,tx,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -299,7 +308,24 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
       ! to compute an approximate solution of a linear system.
       !
       !
+      select case (init_)
+      case('Z') 
         ty = zzero
+      case('Y')
+        call psb_geaxpby(zone,y,zzero,ty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(zone,initu,zzero,ty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
+      
       do i=1, sweeps
         select case(trans_)
         case('N')
@@ -386,7 +412,7 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_apply_vect
   implicit none 
@@ -50,13 +50,15 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                  :: sweeps
   complex(psb_dpk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col, nrow_d, i
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: vx(:) 
   type(psb_z_vect_type) :: vtx, vty, vww
   integer(psb_ipk_)  :: ictxt,np,me, err_act,isz,int_err(5)
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='z_as_smoother_apply', ch_err
 
   call psb_erractionsave(err_act)
@@ -65,6 +67,12 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   ictxt = desc_data%get_context()
   call psb_info(ictxt,me,np)
 
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -135,7 +143,7 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   else if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then 
     !
     ! Shortcut: in this case it's just the same
-    ! as Block Jacobi.
+    ! as Block Jacobi. Moreover, if .not.sv%is_iterative, there's no need to pass init
     !
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
@@ -240,7 +248,7 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
         goto 9999
       end select
 
-      call sm%sv%apply(zone,vtx,zzero,vty,sm%desc_data,trans_,aux,info) 
+      call sm%sv%apply(zone,vtx,zzero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) then
         call psb_errpush(psb_err_internal_error_,name,&
@@ -313,7 +321,23 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
       ! to compute an approximate solution of a linear system.
       !
       !
-      call vty%set(zzero)
+      select case (init_)
+      case('Z') 
+        call vty%zero()
+      case('Y')
+        call psb_geaxpby(zone,y,zzero,vty,sm%desc_data,info)
+      case('U')
+        if (.not.present(initu)) then
+          call psb_errpush(psb_err_internal_error_,name,&
+               & a_err='missing initu to smoother_apply')
+          goto 9999
+        end if
+        call psb_geaxpby(zone,initu,zzero,vty,sm%desc_data,info)
+      case default
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='wrong  init to smoother_apply')
+        goto 9999
+      end select
       
       do i=1, sweeps
         select case(trans_)
@@ -401,7 +425,7 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
         if (info /= psb_success_) exit
 
-        call sm%sv%apply(zone,vww,zzero,vty,sm%desc_data,trans_,aux,info) 
+        call sm%sv%apply(zone,vww,zzero,vty,sm%desc_data,trans_,aux,info,init='Y') 
 
         if (info /= psb_success_) exit
 

mlprec/impl/smoother/mld_z_base_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_z_base_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_z_base_smoother_mod, mld_protect_name =>  mld_z_base_smoother_apply
   implicit none 
@@ -50,6 +50,8 @@ subroutine mld_z_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
   integer(psb_ipk_), intent(in)                :: sweeps
   complex(psb_dpk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='z_base_smoother_apply'
@@ -66,7 +68,7 @@ subroutine mld_z_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo
 
   else
     if (allocated(sm%sv)) then 
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90

@@ -37,8 +37,7 @@
 !!$ 
 !!$
 subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-  
+     &  trans,sweeps,work,info,init,initu)
-     &  trans,sweeps,work,info)
   use psb_base_mod
   use mld_z_base_smoother_mod, mld_protect_name =>  mld_z_base_smoother_apply_vect
   implicit none 
@@ -51,6 +50,8 @@ subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
   integer(psb_ipk_), intent(in)                    :: sweeps
   complex(psb_dpk_),target, intent(inout)            :: work(:)
   integer(psb_ipk_), intent(out)                   :: info
+  character, intent(in), optional                  :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='z_base_smoother_apply'
@@ -67,7 +68,7 @@ subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
 
   else
     if (allocated(sm%sv)) then
-      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info)
+      call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info,init=init, initu=initu)
     else
       info = 1121
     endif

mlprec/impl/smoother/mld_z_jac_smoother_apply.f90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info)
+subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,&
-  
+     & trans,sweeps,work,info,init,initu)
   use psb_base_mod
   use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply
   implicit none 
@@ -50,17 +50,27 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
   integer(psb_ipk_), intent(in)         :: sweeps
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)    :: n_row,n_col
-  complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
+  complex(psb_dpk_), allocatable  :: tx(:),ty(:)
+  complex(psb_dpk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='z_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
@@ -112,8 +122,8 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
     !
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
-  else  if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then 
+  else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -122,38 +132,47 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor
       goto 9999
     endif
 
-  else if (sweeps  > 1) then 
+  else if (sweeps >= 1) then 
-
     !
     !
     ! Apply multiple sweeps of a block-Jacobi solver
     ! to compute an approximate solution of a linear system.
     !
     !
-    allocate(tx(n_col),ty(n_col),stat=info)
+    
-    if (info /= psb_success_) then 
+    call psb_geasb(tx,desc_data,info)
-      info=psb_err_alloc_request_
+    call psb_geasb(ty,desc_data,info)
-      call psb_errpush(info,name,& 
+
-           & i_err=(/2*n_col,izero,izero,izero,izero/),&
+    select case (init_)
-           & a_err='complex(psb_dpk_)')
+    case('Z') 
+      tx(:) = zzero
+    case('Y')
+      call psb_geaxpby(zone,y,zzero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
         goto 9999
       end if
+      call psb_geaxpby(zone,initu,zzero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
-    tx = zzero
-    ty = zzero
     do i=1, sweeps
       !
       ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
       ! block diagonal part and the remaining part of the local matrix
       ! and Y(j) is the approximate solution at sweep j.
       !
-      ty(1:n_row) = x(1:n_row)
+      call psb_geaxpby(zone,x,zzero,ty,desc_data,info)
-      call psb_spmm(-zone,sm%nd,tx,zone,ty,desc_data,info,&
+      call psb_spmm(-zone,sm%nd,tx,zone,ty,desc_data,info,work=aux,trans=trans_)
-           & work=aux,trans=trans_)
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(zone,ty,zzero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(zone,ty,zzero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90

@@ -37,7 +37,7 @@
 !!$ 
 !!$
 subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 
-     & sweeps,work,info)
+     & sweeps,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply_vect
@@ -51,23 +51,31 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
   integer(psb_ipk_), intent(in)                   :: sweeps
   complex(psb_dpk_),target, intent(inout)           :: work(:)
   integer(psb_ipk_), intent(out)                  :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: n_row,n_col
   type(psb_z_vect_type)  :: tx, ty
   complex(psb_dpk_), pointer :: ww(:), aux(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='z_jac_smoother_apply'
 
   call psb_erractionsave(err_act)
 
   info = psb_success_
 
+  
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
   trans_ = psb_toupper(trans)
   select case(trans_)
   case('N')
-  case('T')
+  case('T','C')
-  case('C')
   case default
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
@@ -107,8 +115,6 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     end if
   endif
 
-!!$  write(0,*) 'Jacobi   smoother with ',sweeps
-  
   if (sweeps == 0) then
     
     !
@@ -118,7 +124,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     call psb_geaxpby(alpha,x,beta,y,desc_data,info) 
 
   else   if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then 
-
+    !  if .not.sv%is_iterative, there's no need to pass init
     call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) 
 
     if (info /= psb_success_) then
@@ -134,9 +140,27 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
     ! to compute an approximate solution of a linear system.
     !
     !
-    call tx%bld(x%get_nrows(),mold=x%v)
+    
-    call tx%set(zzero)
+    call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
-    call ty%bld(x%get_nrows(),mold=x%v)
+    call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
+
+    select case (init_)
+    case('Z') 
+      call tx%zero()
+    case('Y')
+      call psb_geaxpby(zone,y,zzero,tx,desc_data,info)
+    case('U')
+      if (.not.present(initu)) then
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='missing initu to smoother_apply')
+        goto 9999
+      end if
+      call psb_geaxpby(zone,initu,zzero,tx,desc_data,info)
+    case default
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='wrong  init to smoother_apply')
+      goto 9999
+    end select
 
     do i=1, sweeps
       !
@@ -149,7 +173,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
 
       if (info /= psb_success_) exit
 
-      call sm%sv%apply(zone,ty,zzero,tx,desc_data,trans_,aux,info) 
+      call sm%sv%apply(zone,ty,zzero,tx,desc_data,trans_,aux,info,init='Y') 
 
       if (info /= psb_success_) exit
     end do

mlprec/impl/solver/mld_c_base_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_base_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_base_solver_mod, mld_protect_name =>  mld_c_base_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_c_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   complex(psb_spk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='c_base_solver_apply'

mlprec/impl/solver/mld_c_base_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_base_solver_mod, mld_protect_name =>  mld_c_base_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_c_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   complex(psb_spk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='c_base_solver_apply_vect'

mlprec/impl/solver/mld_c_bwgs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_gs_solver, mld_protect_name => mld_c_bwgs_solver_apply
@@ -49,13 +50,15 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
-  character(len=20)  :: name='d_bwgs_solver_apply'
+  character(len=20)  :: name='c_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
   ictxt = desc_data%get_ctxt()
@@ -72,6 +75,13 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = czero
+  case('Y')
+    call psb_geaxpby(cone,y,czero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(cone,initu,czero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(cone,y,czero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(cone,x,czero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-cone,sv%l,xit,cone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(cone,sv%u,wv,czero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_c_bwgs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_gs_solver, mld_protect_name => mld_c_bwgs_solver_apply_vect
@@ -49,14 +50,16 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                   :: trans
   complex(psb_spk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_c_vect_type)  :: wv, xit
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
-  character(len=20)  :: name='d_bwgs_solver_apply'
+  character(len=20)  :: name='c_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
   ictxt = desc_data%get_ctxt()
@@ -73,6 +76,13 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(cone,y,czero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(cone,initu,czero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(cone,y,czero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(cone,x,czero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-cone,sv%l,xit,cone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(cone,sv%u,wv,czero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_c_diag_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_diag_solver, mld_protect_name => mld_c_diag_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   complex(psb_spk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)   :: n_row,n_col
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
   
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_c_diag_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_diag_solver, mld_protect_name => mld_c_diag_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   complex(psb_spk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_c_gs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_apply
@@ -49,12 +50,14 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='c_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -72,6 +75,13 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = czero
+  case('Y')
+    call psb_geaxpby(cone,y,czero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(cone,initu,czero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(cone,y,czero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(cone,x,czero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-cone,sv%u,xit,cone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(cone,sv%l,wv,czero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_apply_vect
@@ -49,13 +50,15 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                   :: trans
   complex(psb_spk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_c_vect_type)  :: wv, xit
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='c_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -73,6 +76,13 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(cone,y,czero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(cone,initu,czero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(cone,y,czero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(cone,x,czero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-cone,sv%u,xit,cone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(cone,sv%l,wv,czero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_c_id_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_c_id_solver, mld_protect_name => mld_c_id_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_c_id_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_c_id_solver, mld_protect_name => mld_c_id_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                  :: trans
   complex(psb_spk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_c_ilu_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_ilu_solver, mld_protect_name => mld_c_ilu_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_c_ilu_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_c_ilu_solver, mld_protect_name => mld_c_ilu_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
   character(len=1),intent(in)                   :: trans
   complex(psb_spk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   type(psb_c_vect_type)  :: wv, wv1
@@ -70,6 +73,9 @@ subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_c_mumps_solver_apply.F90

@@ -36,9 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
-
+     & trans,work,info,init,initu)
-subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   use psb_base_mod
   use mld_c_mumps_solver
   implicit none 
@@ -50,6 +49,8 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)          :: trans
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional       :: init
+  complex(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row, n_col, nglob
   complex(psb_spk_), allocatable     :: ww(:)
@@ -70,6 +71,9 @@ subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   nglob = desc_data%get_global_rows()
   n_row = desc_data%get_local_rows()

mlprec/impl/solver/mld_c_mumps_solver_apply_vect.F90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine c_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
-  subroutine c_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+     & trans,work,info,init,initu)
   use psb_base_mod
   use mld_c_mumps_solver
   implicit none 
@@ -49,6 +49,8 @@
   character(len=1),intent(in)          :: trans
   complex(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional                :: init
+  type(psb_c_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: err_act
   character(len=20)  :: name='c_mumps_solver_apply_vect'
@@ -58,6 +60,9 @@
   call psb_erractionsave(err_act)
 
   info = psb_success_
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call x%v%sync()
   call y%v%sync()

mlprec/impl/solver/mld_d_base_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_base_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_base_solver_mod, mld_protect_name =>  mld_d_base_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_d_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   real(psb_dpk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='d_base_solver_apply'

mlprec/impl/solver/mld_d_base_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_base_solver_mod, mld_protect_name =>  mld_d_base_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_d_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   real(psb_dpk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='d_base_solver_apply_vect'

mlprec/impl/solver/mld_d_bwgs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_gs_solver, mld_protect_name => mld_d_bwgs_solver_apply
@@ -49,12 +50,14 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -72,6 +75,13 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = dzero
+  case('Y')
+    call psb_geaxpby(done,y,dzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(done,initu,dzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_d_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(done,y,dzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(done,x,dzero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-done,sv%l,xit,done,wv,desc_data,info,doswap=.false.)
         call psb_spsm(done,sv%u,wv,dzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_d_bwgs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_gs_solver, mld_protect_name => mld_d_bwgs_solver_apply_vect
@@ -49,13 +50,15 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                   :: trans
   real(psb_dpk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_d_vect_type)  :: wv, xit
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -73,6 +76,13 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(done,y,dzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(done,initu,dzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_d_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(done,y,dzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(done,x,dzero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-done,sv%l,xit,done,wv,desc_data,info,doswap=.false.)
         call psb_spsm(done,sv%u,wv,dzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_d_diag_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_diag_solver, mld_protect_name => mld_d_diag_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   real(psb_dpk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)   :: n_row,n_col
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
   
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_d_diag_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_diag_solver, mld_protect_name => mld_d_diag_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   real(psb_dpk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_d_gs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_gs_solver, mld_protect_name => mld_d_gs_solver_apply
@@ -49,12 +50,14 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -72,6 +75,13 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = dzero
+  case('Y')
+    call psb_geaxpby(done,y,dzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(done,initu,dzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_d_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(done,y,dzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(done,x,dzero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-done,sv%u,xit,done,wv,desc_data,info,doswap=.false.)
         call psb_spsm(done,sv%l,wv,dzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_d_gs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_gs_solver, mld_protect_name => mld_d_gs_solver_apply_vect
@@ -49,13 +50,15 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                   :: trans
   real(psb_dpk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_d_vect_type)  :: wv, xit
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='d_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -73,6 +76,13 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(done,y,dzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(done,initu,dzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_d_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(done,y,dzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(done,x,dzero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-done,sv%u,xit,done,wv,desc_data,info,doswap=.false.)
         call psb_spsm(done,sv%l,wv,dzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_d_id_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_d_id_solver, mld_protect_name => mld_d_id_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_d_id_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_d_id_solver, mld_protect_name => mld_d_id_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                  :: trans
   real(psb_dpk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_d_ilu_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_ilu_solver, mld_protect_name => mld_d_ilu_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_d_ilu_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_d_ilu_solver, mld_protect_name => mld_d_ilu_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
   character(len=1),intent(in)                   :: trans
   real(psb_dpk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   type(psb_d_vect_type)  :: wv, wv1
@@ -70,6 +73,9 @@ subroutine mld_d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_d_mumps_solver_apply.F90

@@ -36,9 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
-
+     & trans,work,info,init,initu)
-subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   use psb_base_mod
   use mld_d_mumps_solver
   implicit none 
@@ -50,6 +49,8 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)          :: trans
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional       :: init
+  real(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row, n_col, nglob
   real(psb_dpk_), allocatable     :: ww(:)
@@ -70,6 +71,9 @@ subroutine d_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   nglob = desc_data%get_global_rows()
   n_row = desc_data%get_local_rows()

mlprec/impl/solver/mld_d_mumps_solver_apply_vect.F90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine d_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
-  subroutine d_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+     & trans,work,info,init,initu)
   use psb_base_mod
   use mld_d_mumps_solver
   implicit none 
@@ -49,6 +49,8 @@
   character(len=1),intent(in)          :: trans
   real(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional                :: init
+  type(psb_d_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: err_act
   character(len=20)  :: name='d_mumps_solver_apply_vect'
@@ -58,6 +60,9 @@
   call psb_erractionsave(err_act)
 
   info = psb_success_
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call x%v%sync()
   call y%v%sync()

mlprec/impl/solver/mld_s_base_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_base_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_base_solver_mod, mld_protect_name =>  mld_s_base_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_s_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   real(psb_spk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='s_base_solver_apply'

mlprec/impl/solver/mld_s_base_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_base_solver_mod, mld_protect_name =>  mld_s_base_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_s_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   real(psb_spk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='s_base_solver_apply_vect'

mlprec/impl/solver/mld_s_bwgs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_gs_solver, mld_protect_name => mld_s_bwgs_solver_apply
@@ -49,13 +50,15 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
-  character(len=20)  :: name='d_bwgs_solver_apply'
+  character(len=20)  :: name='s_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
   ictxt = desc_data%get_ctxt()
@@ -72,6 +75,13 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = szero
+  case('Y')
+    call psb_geaxpby(sone,y,szero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(sone,initu,szero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_s_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(sone,y,szero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(sone,x,szero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-sone,sv%l,xit,sone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(sone,sv%u,wv,szero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_s_bwgs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_gs_solver, mld_protect_name => mld_s_bwgs_solver_apply_vect
@@ -49,14 +50,16 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                   :: trans
   real(psb_spk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_s_vect_type)  :: wv, xit
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
-  character(len=20)  :: name='d_bwgs_solver_apply'
+  character(len=20)  :: name='s_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
   ictxt = desc_data%get_ctxt()
@@ -73,6 +76,13 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(sone,y,szero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(sone,initu,szero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_s_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(sone,y,szero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(sone,x,szero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-sone,sv%l,xit,sone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(sone,sv%u,wv,szero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_s_diag_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_diag_solver, mld_protect_name => mld_s_diag_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   real(psb_spk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)   :: n_row,n_col
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
   
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_s_diag_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_diag_solver, mld_protect_name => mld_s_diag_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   real(psb_spk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_s_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_s_gs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_apply
@@ -49,12 +50,14 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='s_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -72,6 +75,13 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = szero
+  case('Y')
+    call psb_geaxpby(sone,y,szero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(sone,initu,szero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(sone,y,szero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(sone,x,szero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-sone,sv%u,xit,sone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(sone,sv%l,wv,szero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_apply_vect
@@ -49,13 +50,15 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                   :: trans
   real(psb_spk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_s_vect_type)  :: wv, xit
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   real(psb_spk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='s_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -73,6 +76,13 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(sone,y,szero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(sone,initu,szero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(sone,y,szero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(sone,x,szero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-sone,sv%u,xit,sone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(sone,sv%l,wv,szero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_s_id_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_s_id_solver, mld_protect_name => mld_s_id_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_s_id_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_s_id_solver, mld_protect_name => mld_s_id_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                  :: trans
   real(psb_spk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_s_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_s_ilu_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_ilu_solver, mld_protect_name => mld_s_ilu_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_s_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_s_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_s_ilu_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_s_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_s_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_s_ilu_solver, mld_protect_name => mld_s_ilu_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_s_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
   character(len=1),intent(in)                   :: trans
   real(psb_spk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   type(psb_s_vect_type)  :: wv, wv1
@@ -70,6 +73,9 @@ subroutine mld_s_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_s_mumps_solver_apply.F90

@@ -36,9 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
-
+     & trans,work,info,init,initu)
-subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   use psb_base_mod
   use mld_s_mumps_solver
   implicit none 
@@ -50,6 +49,8 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)          :: trans
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional       :: init
+  real(psb_spk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row, n_col, nglob
   real(psb_spk_), allocatable     :: ww(:)
@@ -70,6 +71,9 @@ subroutine s_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   nglob = desc_data%get_global_rows()
   n_row = desc_data%get_local_rows()

mlprec/impl/solver/mld_s_mumps_solver_apply_vect.F90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine s_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
-  subroutine s_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+     & trans,work,info,init,initu)
   use psb_base_mod
   use mld_s_mumps_solver
   implicit none 
@@ -49,6 +49,8 @@
   character(len=1),intent(in)          :: trans
   real(psb_spk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional                :: init
+  type(psb_s_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: err_act
   character(len=20)  :: name='s_mumps_solver_apply_vect'
@@ -58,6 +60,9 @@
   call psb_erractionsave(err_act)
 
   info = psb_success_
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call x%v%sync()
   call y%v%sync()

mlprec/impl/solver/mld_z_base_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_base_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_base_solver_mod, mld_protect_name =>  mld_z_base_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_z_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   complex(psb_dpk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='z_base_solver_apply'

mlprec/impl/solver/mld_z_base_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_base_solver_mod, mld_protect_name =>  mld_z_base_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_z_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   complex(psb_dpk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_) :: err_act
   character(len=20) :: name='z_base_solver_apply_vect'

mlprec/impl/solver/mld_z_bwgs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_gs_solver, mld_protect_name => mld_z_bwgs_solver_apply
@@ -49,13 +50,15 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
-  character(len=20)  :: name='d_bwgs_solver_apply'
+  character(len=20)  :: name='z_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
   ictxt = desc_data%get_ctxt()
@@ -72,6 +75,13 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = zzero
+  case('Y')
+    call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(zone,initu,zzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_z_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-zone,sv%l,xit,zone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(zone,sv%u,wv,zzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_z_bwgs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_gs_solver, mld_protect_name => mld_z_bwgs_solver_apply_vect
@@ -49,14 +50,16 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                   :: trans
   complex(psb_dpk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_z_vect_type)  :: wv, xit
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
-  character(len=20)  :: name='d_bwgs_solver_apply'
+  character(len=20)  :: name='z_bwgs_solver_apply'
 
   call psb_erractionsave(err_act)
   ictxt = desc_data%get_ctxt()
@@ -73,6 +76,13 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(zone,initu,zzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_z_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
         ! Update with L. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-zone,sv%l,xit,zone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(zone,sv%u,wv,zzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_z_diag_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_diag_solver, mld_protect_name => mld_z_diag_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)                :: trans
   complex(psb_dpk_),target, intent(inout)      :: work(:)
   integer(psb_ipk_), intent(out)             :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)   :: n_row,n_col
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
   
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_z_diag_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_diag_solver, mld_protect_name => mld_z_diag_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
   character(len=1),intent(in)                    :: trans
   complex(psb_dpk_),target, intent(inout)          :: work(:)
   integer(psb_ipk_), intent(out)                 :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_z_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,i
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_z_gs_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     &trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_apply
@@ -49,12 +50,14 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col, itx
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:),wv(:),xit(:)
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='z_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -72,6 +75,13 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -110,6 +120,23 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
 
   call psb_geasb(wv,desc_data,info) 
   call psb_geasb(xit,desc_data,info) 
+  select case (init_)
+  case('Z') 
+    xit(:) = zzero
+  case('Y')
+    call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(zone,initu,zzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -118,21 +145,12 @@ subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-zone,sv%u,xit,zone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(zone,sv%l,wv,zzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_apply_vect
@@ -49,13 +50,15 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                   :: trans
   complex(psb_dpk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col, itx
   type(psb_z_vect_type)  :: wv, xit
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
   complex(psb_dpk_), allocatable :: temp(:)
   integer(psb_ipk_)   :: ictxt,np,me,i, err_act
-  character          :: trans_
+  character          :: trans_, init_
   character(len=20)  :: name='z_gs_solver_apply'
 
   call psb_erractionsave(err_act)
@@ -73,6 +76,13 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     goto 9999
   end select
   
+  if (present(init)) then
+    init_ = psb_toupper(init)
+  else
+    init_='Z'
+  end if
+
+
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()
 
@@ -113,6 +123,23 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
 
   call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
   call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) 
+  select case (init_)
+  case('Z') 
+    call xit%zero()
+  case('Y')
+    call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
+  case('U')
+    if (.not.present(initu)) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='missing initu to smoother_apply')
+      goto 9999
+    end if
+    call psb_geaxpby(zone,initu,zzero,xit,desc_data,info)
+  case default
+    call psb_errpush(psb_err_internal_error_,name,&
+         & a_err='wrong  init to smoother_apply')
+    goto 9999
+  end select
   
   select case(trans_)
   case('N')
@@ -121,21 +148,12 @@ subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
       ! Fixed number of iterations
       !
       !
-      !  WARNING: this is not completely satisfactory. We are assuming here Y
-      !  as the initial guess, but this is only working if we are called from the
-      !  current JAC smoother loop. A good solution would be to have a separate
-      !  input argument as the initial guess
-      !  
-!!$      write(0,*) 'GS Iteration with ',sv%sweeps
-      call psb_geaxpby(zone,y,zzero,xit,desc_data,info)
       do itx=1,sv%sweeps
         call psb_geaxpby(zone,x,zzero,wv,desc_data,info)
         ! Update with U. The off-diagonal block is taken care
         ! from the Jacobi smoother, hence this is purely local. 
         call psb_spmm(-zone,sv%u,xit,zone,wv,desc_data,info,doswap=.false.)
         call psb_spsm(zone,sv%l,wv,zzero,xit,desc_data,info)
-!!$        temp = xit%get_vect()
-!!$        write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row)
       end do
       
       call psb_geaxpby(alpha,xit,beta,y,desc_data,info)

mlprec/impl/solver/mld_z_id_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_z_id_solver, mld_protect_name => mld_z_id_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_z_id_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
 
   use psb_base_mod
   use mld_z_id_solver, mld_protect_name => mld_z_id_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
   character(len=1),intent(in)                  :: trans
   complex(psb_dpk_),target, intent(inout)        :: work(:)
   integer(psb_ipk_), intent(out)               :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)  :: n_row,n_col
   integer(psb_ipk_)  :: ictxt,np,me,i, err_act
@@ -68,6 +71,9 @@ subroutine mld_z_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,inf
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call psb_geaxpby(alpha,x,beta,y,desc_data,info)    
 

mlprec/impl/solver/mld_z_ilu_solver_apply.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_ilu_solver, mld_protect_name => mld_z_ilu_solver_apply
@@ -49,6 +50,8 @@ subroutine mld_z_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)           :: trans
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)        :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row,n_col
   complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
@@ -68,6 +71,9 @@ subroutine mld_z_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_z_ilu_solver_apply_vect.f90

@@ -36,7 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-subroutine mld_z_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+subroutine mld_z_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+     & trans,work,info,init,initu)
   
   use psb_base_mod
   use mld_z_ilu_solver, mld_protect_name => mld_z_ilu_solver_apply_vect
@@ -49,6 +50,8 @@ subroutine mld_z_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
   character(len=1),intent(in)                   :: trans
   complex(psb_dpk_),target, intent(inout)         :: work(:)
   integer(psb_ipk_), intent(out)                :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)   :: n_row,n_col
   type(psb_z_vect_type)  :: wv, wv1
@@ -70,6 +73,9 @@ subroutine mld_z_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,in
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   n_row = desc_data%get_local_rows()
   n_col = desc_data%get_local_cols()

mlprec/impl/solver/mld_z_mumps_solver_apply.F90

@@ -36,9 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,&
-
+     & trans,work,info,init,initu)
-subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   use psb_base_mod
   use mld_z_mumps_solver
   implicit none 
@@ -50,6 +49,8 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
   character(len=1),intent(in)          :: trans
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional       :: init
+  complex(psb_dpk_),intent(inout), optional :: initu(:)
 
   integer(psb_ipk_)  :: n_row, n_col, nglob
   complex(psb_dpk_), allocatable     :: ww(:)
@@ -70,6 +71,9 @@ subroutine z_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
     call psb_errpush(psb_err_iarg_invalid_i_,name)
     goto 9999
   end select
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   nglob = desc_data%get_global_rows()
   n_row = desc_data%get_local_rows()

mlprec/impl/solver/mld_z_mumps_solver_apply_vect.F90

@@ -36,8 +36,8 @@
 !!$  POSSIBILITY OF SUCH DAMAGE.
 !!$ 
 !!$
-
+subroutine z_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
-  subroutine z_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+     & trans,work,info,init,initu)
   use psb_base_mod
   use mld_z_mumps_solver
   implicit none 
@@ -49,6 +49,8 @@
   character(len=1),intent(in)          :: trans
   complex(psb_dpk_),target, intent(inout) :: work(:)
   integer(psb_ipk_), intent(out)       :: info
+  character, intent(in), optional                :: init
+  type(psb_z_vect_type),intent(inout), optional   :: initu
 
   integer(psb_ipk_)    :: err_act
   character(len=20)  :: name='z_mumps_solver_apply_vect'
@@ -58,6 +60,9 @@
   call psb_erractionsave(err_act)
 
   info = psb_success_
+  !
+  ! For non-iterative solvers, init and initu are ignored.
+  !
 
   call x%v%sync()
   call y%v%sync()

mlprec/mld_base_prec_type.F90

@@ -545,16 +545,17 @@ contains
     integer(psb_ipk_), intent(out)            :: info
 
     info = psb_success_
-    if (pm%ml_type>mld_no_ml_) then
+    if ((pm%ml_type>=mld_no_ml_).and.(pm%ml_type<=mld_max_ml_type_)) then
 
       write(iout,*) '  Multilevel type: ',&
            &   ml_names(pm%ml_type)
-      write(iout,*) '  Smoother position: ',&
+      select case (pm%ml_type)
-           & smooth_pos_names(pm%smoother_pos)
+      case (mld_add_ml_)
-      if (pm%ml_type == mld_add_ml_) then
         write(iout,*) '  Number of smoother sweeps : ',&
              & pm%sweeps
-      else 
+      case (mld_mult_ml_) 
+        write(iout,*) '  Smoother position: ',&
+             & smooth_pos_names(pm%smoother_pos)
         select case (pm%smoother_pos)
         case (mld_pre_smooth_)
           write(iout,*) '  Number of smoother sweeps : ',&
@@ -568,7 +569,13 @@ contains
                &  '  post: ',&
                &  pm%sweeps_post
         end select
-      end if
+      case (mld_vcycle_ml_, mld_wcycle_ml_, mld_kcycle_ml_, mld_kcyclesym_ml_)
+        write(iout,*) '  Number of smoother sweeps : pre: ',&
+             &  pm%sweeps_pre ,&
+             &  '  post: ',&
+             &  pm%sweeps_post
+      end select
+      
       write(iout,*) '  Aggregation: ', &
            &   aggr_names(pm%aggr_alg)
       write(iout,*) '               with initial ordering: ',&
@@ -586,6 +593,9 @@ contains
           write(iout,*) '  Damping omega computation: unknown value in iprcparm!!'
         end if
       end if
+    else
+      write(iout,*) '  Multilevel type: Unkonwn value. Something is amis....',&
+           & pm%ml_type           
     end if
     
     return

mlprec/mld_c_as_smoother.f90

@@ -97,7 +97,7 @@ module mld_c_as_smoother
   
   interface 
     subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-      & trans,sweeps,work,info)
+      & trans,sweeps,work,info,init,initu)
       import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
            & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, &
            & psb_desc_type, psb_ipk_
@@ -111,12 +111,14 @@ module mld_c_as_smoother
       integer(psb_ipk_), intent(in)                  :: sweeps
       complex(psb_spk_),target, intent(inout)          :: work(:)
       integer(psb_ipk_), intent(out)                 :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_as_smoother_apply_vect
   end interface
   
   interface
     subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,& 
-         & trans,sweeps,work,info)
+         & trans,sweeps,work,info,init,initu)
       import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
            & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_,&
            & psb_desc_type, psb_ipk_
@@ -130,6 +132,8 @@ module mld_c_as_smoother
       integer(psb_ipk_), intent(in)         :: sweeps
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_as_smoother_apply
   end interface
   

mlprec/mld_c_base_smoother_mod.f90

@@ -125,7 +125,7 @@ module mld_c_base_smoother_mod
 
   interface 
     subroutine mld_c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,& 
-         & trans,sweeps,work,info)
+         & trans,sweeps,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & mld_c_base_smoother_type, psb_ipk_
@@ -138,12 +138,14 @@ module mld_c_base_smoother_mod
       integer(psb_ipk_), intent(in)                :: sweeps
       complex(psb_spk_),target, intent(inout)        :: work(:)
       integer(psb_ipk_), intent(out)               :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_base_smoother_apply
   end interface
   
   interface 
     subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-         &  trans,sweeps,work,info)
+         &  trans,sweeps,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & mld_c_base_smoother_type, psb_ipk_
@@ -156,6 +158,8 @@ module mld_c_base_smoother_mod
       integer(psb_ipk_), intent(in)                    :: sweeps
       complex(psb_spk_),target, intent(inout)            :: work(:)
       integer(psb_ipk_), intent(out)                   :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_base_smoother_apply_vect
   end interface
   

mlprec/mld_c_base_solver_mod.f90

@@ -116,7 +116,8 @@ module mld_c_base_solver_mod
 
 
   interface  
-    subroutine mld_c_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_base_solver_apply(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
        & mld_c_base_solver_type, psb_ipk_
@@ -129,12 +130,15 @@ module mld_c_base_solver_mod
       character(len=1),intent(in)                :: trans
       complex(psb_spk_),target, intent(inout)      :: work(:)
       integer(psb_ipk_), intent(out)             :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_base_solver_apply
   end interface 
   
       
   interface 
-    subroutine mld_c_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & mld_c_base_solver_type, psb_ipk_
@@ -147,6 +151,8 @@ module mld_c_base_solver_mod
       character(len=1),intent(in)                    :: trans
       complex(psb_spk_),target, intent(inout)          :: work(:)
       integer(psb_ipk_), intent(out)                 :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_base_solver_apply_vect
   end interface
   

mlprec/mld_c_diag_solver.f90

@@ -71,7 +71,7 @@ module mld_c_diag_solver
 
   interface 
     subroutine mld_c_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,& 
-         & trans,work,info)
+         & trans,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
        & mld_c_diag_solver_type, psb_ipk_
@@ -83,11 +83,14 @@ module mld_c_diag_solver
       character(len=1),intent(in)                    :: trans
       complex(psb_spk_),target, intent(inout)          :: work(:)
       integer(psb_ipk_), intent(out)                 :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_diag_solver_apply_vect
   end interface
   
   interface 
-    subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
        & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
        & mld_c_diag_solver_type, psb_ipk_
@@ -99,6 +102,8 @@ module mld_c_diag_solver
       character(len=1),intent(in)                :: trans
       complex(psb_spk_),target, intent(inout)      :: work(:)
       integer(psb_ipk_), intent(out)             :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_diag_solver_apply
   end interface
   

mlprec/mld_c_gs_solver.f90

@@ -96,7 +96,8 @@ module mld_c_gs_solver
 
 
   interface 
-    subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_gs_solver_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -108,8 +109,11 @@ module mld_c_gs_solver
       character(len=1),intent(in)                   :: trans
       complex(psb_spk_),target, intent(inout)         :: work(:)
       integer(psb_ipk_), intent(out)                :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_gs_solver_apply_vect
-    subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_bwgs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_bwgs_solver_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -121,11 +125,13 @@ module mld_c_gs_solver
       character(len=1),intent(in)                   :: trans
       complex(psb_spk_),target, intent(inout)         :: work(:)
       integer(psb_ipk_), intent(out)                :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_bwgs_solver_apply_vect
   end interface
 
   interface 
-    subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_gs_solver_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -137,8 +143,11 @@ module mld_c_gs_solver
       character(len=1),intent(in)           :: trans
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_gs_solver_apply
-    subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_bwgs_solver_apply(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_bwgs_solver_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -150,6 +159,8 @@ module mld_c_gs_solver
       character(len=1),intent(in)           :: trans
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_bwgs_solver_apply
   end interface
 

mlprec/mld_c_id_solver.f90

@@ -64,7 +64,8 @@ module mld_c_id_solver
        &  c_id_solver_descr
 
   interface 
-    subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & 
            & mld_c_id_solver_type, psb_ipk_
@@ -76,11 +77,14 @@ module mld_c_id_solver
       character(len=1),intent(in)                  :: trans
       complex(psb_spk_),target, intent(inout)        :: work(:)
       integer(psb_ipk_), intent(out)               :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_id_solver_apply_vect
   end interface
   
   interface 
-    subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, psb_cspmat_type,  psb_c_base_sparse_mat, &
            & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
            & mld_c_id_solver_type, psb_ipk_
@@ -92,6 +96,8 @@ module mld_c_id_solver
       character(len=1),intent(in)          :: trans
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)                 :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_id_solver_apply
   end interface
 

mlprec/mld_c_ilu_solver.f90

@@ -88,7 +88,8 @@ module mld_c_ilu_solver
 
 
   interface 
-    subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_ilu_solver_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -100,11 +101,14 @@ module mld_c_ilu_solver
       character(len=1),intent(in)                   :: trans
       complex(psb_spk_),target, intent(inout)         :: work(:)
       integer(psb_ipk_), intent(out)                :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_ilu_solver_apply_vect
   end interface
 
   interface 
-    subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_ilu_solver_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -116,6 +120,8 @@ module mld_c_ilu_solver
       character(len=1),intent(in)           :: trans
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_ilu_solver_apply
   end interface
 

mlprec/mld_c_jac_smoother.f90

@@ -74,7 +74,7 @@ module mld_c_jac_smoother
 
   interface 
     subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& 
-         & sweeps,work,info)
+         & sweeps,work,info,init,initu)
       import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
            & psb_ipk_
@@ -88,12 +88,14 @@ module mld_c_jac_smoother
       integer(psb_ipk_), intent(in)                   :: sweeps
       complex(psb_spk_),target, intent(inout)           :: work(:)
       integer(psb_ipk_), intent(out)                  :: info
+      character, intent(in), optional                :: init
+      type(psb_c_vect_type),intent(inout), optional   :: initu
     end subroutine mld_c_jac_smoother_apply_vect
   end interface
   
   interface 
     subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,& 
-         & sweeps,work,info)
+         & sweeps,work,info,init,initu)
       import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, &
            & psb_ipk_
@@ -106,6 +108,8 @@ module mld_c_jac_smoother
       integer(psb_ipk_), intent(in)         :: sweeps
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      complex(psb_spk_),intent(inout), optional :: initu(:)
     end subroutine mld_c_jac_smoother_apply
   end interface
   

mlprec/mld_c_mumps_solver.F90

@@ -93,7 +93,8 @@ module mld_c_mumps_solver
 #endif
 
   interface 
-    subroutine c_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine c_mumps_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+         & trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_mumps_solver_type, psb_c_vect_type, psb_dpk_, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -105,14 +106,13 @@ module mld_c_mumps_solver
       character(len=1),intent(in)          :: trans
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer, intent(out)                 :: info
-
+      character, intent(in), optional                :: init
-      integer(psb_ipk_)    :: err_act
+      type(psb_c_vect_type),intent(inout), optional   :: initu
-      character(len=20)  :: name='c_mumps_solver_apply_vect'
     end subroutine c_mumps_solver_apply_vect
   end interface
 
   interface
-    subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+    subroutine c_mumps_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info,init,initu)
       import :: psb_desc_type, mld_c_mumps_solver_type, psb_c_vect_type, psb_dpk_, psb_spk_, &
            & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_
       implicit none 
@@ -124,13 +124,8 @@ module mld_c_mumps_solver
       character(len=1),intent(in)          :: trans
       complex(psb_spk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)                 :: info
-
+      character, intent(in), optional       :: init
-      integer(psb_ipk_)    :: n_row, n_col, nglob
+      complex(psb_spk_),intent(inout), optional :: initu(:)
-      complex(psb_spk_), pointer     :: ww(:)
-      complex(psb_spk_), allocatable, target :: gx(:)
-      integer(psb_ipk_)  :: ictxt,np,me,i, err_act
-      character          :: trans_
-      character(len=20)  :: name='c_mumps_solver_apply'
     end subroutine c_mumps_solver_apply
   end interface
 

mlprec/mld_c_slu_solver.F90

@@ -116,7 +116,8 @@ module mld_c_slu_solver
 
 contains
 
-  subroutine c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+  subroutine c_slu_solver_apply(alpha,sv,x,beta,y,desc_data,&
+       & trans,work,info,init,initu)
     use psb_base_mod
     implicit none 
     type(psb_desc_type), intent(in)      :: desc_data
@@ -127,6 +128,8 @@ contains
     character(len=1),intent(in)          :: trans
     complex(psb_spk_),target, intent(inout) :: work(:)
     integer, intent(out)                 :: info
+    character, intent(in), optional       :: init
+    complex(psb_spk_),intent(inout), optional :: initu(:)
 
     integer    :: n_row,n_col
     complex(psb_spk_), pointer :: ww(:)
@@ -146,6 +149,9 @@ contains
       call psb_errpush(psb_err_iarg_invalid_i_,name)
       goto 9999
     end select
+    !
+    ! For non-iterative solvers, init and initu are ignored.
+    !
 
     n_row = desc_data%get_local_rows()
     n_col = desc_data%get_local_cols()
@@ -198,7 +204,8 @@ contains
 
   end subroutine c_slu_solver_apply
   
-  subroutine c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+  subroutine c_slu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+       & trans,work,info,init,initu)
     use psb_base_mod
     implicit none 
     type(psb_desc_type), intent(in)      :: desc_data
@@ -209,6 +216,8 @@ contains
     character(len=1),intent(in)          :: trans
     complex(psb_spk_),target, intent(inout) :: work(:)
     integer, intent(out)                 :: info
+    character, intent(in), optional                :: init
+    type(psb_c_vect_type),intent(inout), optional   :: initu
 
     integer    :: err_act
     character(len=20)  :: name='c_slu_solver_apply_vect'
@@ -216,6 +225,9 @@ contains
     call psb_erractionsave(err_act)
 
     info = psb_success_
+    !
+    ! For non-iterative solvers, init and initu are ignored.
+    !
 
     call x%v%sync()
     call y%v%sync()

mlprec/mld_c_sludist_solver.F90

@@ -113,7 +113,8 @@ module mld_c_sludist_solver
 
 contains
 
-  subroutine c_sludist_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+  subroutine c_sludist_solver_apply(alpha,sv,x,beta,y,desc_data,&
+       & trans,work,info,init,initu)
     use psb_base_mod
     implicit none 
     type(psb_desc_type), intent(in)      :: desc_data
@@ -124,6 +125,8 @@ contains
     character(len=1),intent(in)          :: trans
     complex(psb_spk_),target, intent(inout) :: work(:)
     integer, intent(out)                 :: info
+    character, intent(in), optional       :: init
+    complex(psb_spk_),intent(inout), optional :: initu(:)
 
     integer    :: n_row,n_col
     complex(psb_spk_), pointer :: ww(:)
@@ -143,6 +146,10 @@ contains
       call psb_errpush(psb_err_iarg_invalid_i_,name)
       goto 9999
     end select
+    !
+    ! For non-iterative solvers, init and initu are ignored.
+    !
+    
 
     n_row = desc_data%get_local_rows()
     n_col = desc_data%get_local_cols()
@@ -197,7 +204,8 @@ contains
 
   end subroutine c_sludist_solver_apply
 
-  subroutine c_sludist_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+  subroutine c_sludist_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+       & trans,work,info,init,initu)
     use psb_base_mod
     implicit none 
     type(psb_desc_type), intent(in)      :: desc_data
@@ -208,6 +216,8 @@ contains
     character(len=1),intent(in)          :: trans
     complex(psb_spk_),target, intent(inout) :: work(:)
     integer, intent(out)                 :: info
+    character, intent(in), optional                :: init
+    type(psb_c_vect_type),intent(inout), optional   :: initu
 
     integer    :: err_act
     character(len=20)  :: name='c_sludist_solver_apply_vect'
@@ -215,6 +225,10 @@ contains
     call psb_erractionsave(err_act)
 
     info = psb_success_
+    !
+    ! For non-iterative solvers, init and initu are ignored.
+    !
+    
 
     call x%v%sync()
     call y%v%sync()

mlprec/mld_c_umf_solver.F90

@@ -116,7 +116,8 @@ module mld_c_umf_solver
 
 contains
 
-  subroutine c_umf_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info)
+  subroutine c_umf_solver_apply(alpha,sv,x,beta,y,desc_data,&
+       & trans,work,info,init,initu)
     use psb_base_mod
     implicit none 
     type(psb_desc_type), intent(in)      :: desc_data
@@ -127,6 +128,8 @@ contains
     character(len=1),intent(in)          :: trans
     complex(psb_spk_),target, intent(inout) :: work(:)
     integer, intent(out)                 :: info
+    character, intent(in), optional       :: init
+    complex(psb_spk_),intent(inout), optional :: initu(:)
 
     integer    :: n_row,n_col
     complex(psb_spk_), pointer :: ww(:)
@@ -146,6 +149,9 @@ contains
       call psb_errpush(psb_err_iarg_invalid_i_,name)
       goto 9999
     end select
+    !
+    ! For non-iterative solvers, init and initu are ignored.
+    !
 
     n_row = desc_data%get_local_rows()
     n_col = desc_data%get_local_cols()
@@ -202,7 +208,8 @@ contains
 
   end subroutine c_umf_solver_apply
 
-  subroutine c_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info)
+  subroutine c_umf_solver_apply_vect(alpha,sv,x,beta,y,desc_data,&
+       & trans,work,info,init,initu)
     use psb_base_mod
     implicit none 
     type(psb_desc_type), intent(in)      :: desc_data
@@ -213,6 +220,8 @@ contains
     character(len=1),intent(in)          :: trans
     complex(psb_spk_),target, intent(inout) :: work(:)
     integer, intent(out)                 :: info
+    character, intent(in), optional                :: init
+    type(psb_c_vect_type),intent(inout), optional   :: initu
 
     integer    :: err_act
     character(len=20)  :: name='c_umf_solver_apply_vect'
@@ -220,6 +229,9 @@ contains
     call psb_erractionsave(err_act)
 
     info = psb_success_
+    !
+    ! For non-iterative solvers, init and initu are ignored.
+    !
 
     call x%v%sync()
     call y%v%sync()

mlprec/mld_d_as_smoother.f90

@@ -97,7 +97,7 @@ module mld_d_as_smoother
   
   interface 
     subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-      & trans,sweeps,work,info)
+      & trans,sweeps,work,info,init,initu)
       import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
            & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, &
            & psb_desc_type, psb_ipk_
@@ -111,12 +111,14 @@ module mld_d_as_smoother
       integer(psb_ipk_), intent(in)                  :: sweeps
       real(psb_dpk_),target, intent(inout)          :: work(:)
       integer(psb_ipk_), intent(out)                 :: info
+      character, intent(in), optional                :: init
+      type(psb_d_vect_type),intent(inout), optional   :: initu
     end subroutine mld_d_as_smoother_apply_vect
   end interface
   
   interface
     subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,& 
-         & trans,sweeps,work,info)
+         & trans,sweeps,work,info,init,initu)
       import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
            & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_,&
            & psb_desc_type, psb_ipk_
@@ -130,6 +132,8 @@ module mld_d_as_smoother
       integer(psb_ipk_), intent(in)         :: sweeps
       real(psb_dpk_),target, intent(inout) :: work(:)
       integer(psb_ipk_), intent(out)        :: info
+      character, intent(in), optional       :: init
+      real(psb_dpk_),intent(inout), optional :: initu(:)
     end subroutine mld_d_as_smoother_apply
   end interface
   

mlprec/mld_d_base_smoother_mod.f90

@@ -125,7 +125,7 @@ module mld_d_base_smoother_mod
 
   interface 
     subroutine mld_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,& 
-         & trans,sweeps,work,info)
+         & trans,sweeps,work,info,init,initu)
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
            & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
            & mld_d_base_smoother_type, psb_ipk_
@@ -138,12 +138,14 @@ module mld_d_base_smoother_mod
       integer(psb_ipk_), intent(in)                :: sweeps
       real(psb_dpk_),target, intent(inout)        :: work(:)
       integer(psb_ipk_), intent(out)               :: info
+      character, intent(in), optional       :: init
+      real(psb_dpk_),intent(inout), optional :: initu(:)
     end subroutine mld_d_base_smoother_apply
   end interface
   
   interface 
     subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
-         &  trans,sweeps,work,info)
+         &  trans,sweeps,work,info,init,initu)
       import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
            & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
            & mld_d_base_smoother_type, psb_ipk_
@@ -156,6 +158,8 @@ module mld_d_base_smoother_mod
       integer(psb_ipk_), intent(in)                    :: sweeps
       real(psb_dpk_),target, intent(inout)            :: work(:)
       integer(psb_ipk_), intent(out)                   :: info
+      character, intent(in), optional                :: init
+      type(psb_d_vect_type),intent(inout), optional   :: initu
     end subroutine mld_d_base_smoother_apply_vect
   end interface