New wrk methods in psblas%prec

krylov/psb_ckrylov.f90

@@ -151,7 +151,7 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   procedure(psb_ckryl_rest_vect) :: psb_crgmres_vect, psb_ccgstabl_vect, psb_cgcr_vect
   procedure(psb_ckryl_cond_vect) :: psb_ccg_vect, psb_cfcg_vect
 
-  
+  logical           :: do_alloc_wrk
   integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_
   character(len=20)             :: name
 
@@ -172,6 +172,9 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
     itrace_ = -1
   end if
 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a)
+
   select case(psb_toupper(method))
   case('CG') 
     call  psb_ccg_vect(a,prec,b,x,eps,desc_a,info,&
@@ -205,6 +208,8 @@ Subroutine psb_ckrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   end select
 
+  if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
+  
   if(info /= psb_success_) then
     call psb_errpush(info,name)
     goto 9999

krylov/psb_dkrylov.f90

@@ -151,7 +151,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   procedure(psb_dkryl_rest_vect) :: psb_drgmres_vect, psb_dcgstabl_vect, psb_dgcr_vect
   procedure(psb_dkryl_cond_vect) :: psb_dcg_vect, psb_dfcg_vect
 
-  
+  logical           :: do_alloc_wrk
   integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_
   character(len=20)             :: name
 
@@ -172,6 +172,9 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
     itrace_ = -1
   end if
 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a)
+
   select case(psb_toupper(method))
   case('CG') 
     call  psb_dcg_vect(a,prec,b,x,eps,desc_a,info,&
@@ -205,6 +208,8 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   end select
 
+  if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
+  
   if(info /= psb_success_) then
     call psb_errpush(info,name)
     goto 9999

krylov/psb_skrylov.f90

@@ -151,7 +151,7 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   procedure(psb_skryl_rest_vect) :: psb_srgmres_vect, psb_scgstabl_vect, psb_sgcr_vect
   procedure(psb_skryl_cond_vect) :: psb_scg_vect, psb_sfcg_vect
 
-  
+  logical           :: do_alloc_wrk
   integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_
   character(len=20)             :: name
 
@@ -172,6 +172,9 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
     itrace_ = -1
   end if
 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a)
+
   select case(psb_toupper(method))
   case('CG') 
     call  psb_scg_vect(a,prec,b,x,eps,desc_a,info,&
@@ -205,6 +208,8 @@ Subroutine psb_skrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   end select
 
+  if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
+  
   if(info /= psb_success_) then
     call psb_errpush(info,name)
     goto 9999

krylov/psb_zkrylov.f90

@@ -151,7 +151,7 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
   procedure(psb_zkryl_rest_vect) :: psb_zrgmres_vect, psb_zcgstabl_vect, psb_zgcr_vect
   procedure(psb_zkryl_cond_vect) :: psb_zcg_vect, psb_zfcg_vect
 
-  
+  logical           :: do_alloc_wrk
   integer(psb_ipk_) :: ictxt,me,np,err_act, itrace_
   character(len=20)             :: name
 
@@ -172,6 +172,9 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
     itrace_ = -1
   end if
 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v,desc=desc_a)
+
   select case(psb_toupper(method))
   case('CG') 
     call  psb_zcg_vect(a,prec,b,x,eps,desc_a,info,&
@@ -205,6 +208,8 @@ Subroutine psb_zkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
          &itmax,iter,err,itrace=itrace_,istop=istop)
   end select
 
+  if ((info==psb_success_).and.do_alloc_wrk) call prec%free_wrk(info)
+  
   if(info /= psb_success_) then
     call psb_errpush(info,name)
     goto 9999

prec/impl/psb_c_bjacprec_impl.f90

@@ -90,6 +90,7 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
   integer(psb_ipk_) :: ictxt,np,me
   integer(psb_ipk_) :: err_act, ierr(5)
   integer(psb_ipk_) :: debug_level, debug_unit
+  logical            :: do_alloc_wrk
   character          :: trans_
   character(len=20)  :: name='c_bjac_prec_apply'
   character(len=20)  :: ch_err
@@ -154,55 +155,57 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
     goto 9999      
   end if
 
-  call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
-  call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) 
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
-
+
-  select case(prec%iprcparm(psb_f_type_))
+  associate (wv => prec%wrk(1), wv1 => prec%wrk(2))
-  case(psb_f_ilu_n_) 
+    
-
+    select case(prec%iprcparm(psb_f_type_))
-    select case(trans_)
+    case(psb_f_ilu_n_) 
-    case('N')
+      
-      call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,&
+      select case(trans_)
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
+      case('N')
-      if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
+        call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('T')
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-      call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,&
+        
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
+      case('T')
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
+        call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('C')
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-
+        
-      call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,&
+      case('C')
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        
-
+        call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,&
-      call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_)
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-  
+        
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
+        call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_)
-           & beta,y,desc_data,info,&
+        
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
-
+             & beta,y,desc_data,info,&
-    end select
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-    if (info /= psb_success_) then 
+        
-      ch_err="psb_spsm"
+      end select
+      if (info /= psb_success_) then 
+        ch_err="psb_spsm"
+        goto 9999
+      end if
+      
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='Invalid factorization')
       goto 9999
-    end if
+    end select
-
-
-  case default
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='Invalid factorization')
-    goto 9999
-  end select
 
+  end associate
+  
   call psb_halo(y,desc_data,info,data=psb_comm_mov_)
 
-  call wv%free(info)
+  if (do_alloc_wrk) call prec%free_wrk(info)
-  call wv1%free(info)
   if (n_col <= size(work)) then 
     if ((4*n_col+n_col) <= size(work)) then 
     else

prec/impl/psb_d_bjacprec_impl.f90

@@ -90,6 +90,7 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
   integer(psb_ipk_) :: ictxt,np,me
   integer(psb_ipk_) :: err_act, ierr(5)
   integer(psb_ipk_) :: debug_level, debug_unit
+  logical            :: do_alloc_wrk
   character          :: trans_
   character(len=20)  :: name='d_bjac_prec_apply'
   character(len=20)  :: ch_err
@@ -154,55 +155,57 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
     goto 9999      
   end if
 
-  call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
-  call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) 
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
-
+
-  select case(prec%iprcparm(psb_f_type_))
+  associate (wv => prec%wrk(1), wv1 => prec%wrk(2))
-  case(psb_f_ilu_n_) 
+    
-
+    select case(prec%iprcparm(psb_f_type_))
-    select case(trans_)
+    case(psb_f_ilu_n_) 
-    case('N')
+      
-      call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,&
+      select case(trans_)
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
+      case('N')
-      if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
+        call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('T')
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-      call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,&
+        
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
+      case('T')
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
+        call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('C')
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-
+        
-      call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,&
+      case('C')
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        
-
+        call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,&
-      call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_)
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-  
+        
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
+        call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_)
-           & beta,y,desc_data,info,&
+        
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
-
+             & beta,y,desc_data,info,&
-    end select
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-    if (info /= psb_success_) then 
+        
-      ch_err="psb_spsm"
+      end select
+      if (info /= psb_success_) then 
+        ch_err="psb_spsm"
+        goto 9999
+      end if
+      
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='Invalid factorization')
       goto 9999
-    end if
+    end select
-
-
-  case default
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='Invalid factorization')
-    goto 9999
-  end select
 
+  end associate
+  
   call psb_halo(y,desc_data,info,data=psb_comm_mov_)
 
-  call wv%free(info)
+  if (do_alloc_wrk) call prec%free_wrk(info)
-  call wv1%free(info)
   if (n_col <= size(work)) then 
     if ((4*n_col+n_col) <= size(work)) then 
     else

prec/impl/psb_s_bjacprec_impl.f90

@@ -90,6 +90,7 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
   integer(psb_ipk_) :: ictxt,np,me
   integer(psb_ipk_) :: err_act, ierr(5)
   integer(psb_ipk_) :: debug_level, debug_unit
+  logical            :: do_alloc_wrk
   character          :: trans_
   character(len=20)  :: name='s_bjac_prec_apply'
   character(len=20)  :: ch_err
@@ -154,55 +155,57 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
     goto 9999      
   end if
 
-  call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
-  call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) 
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
-
+
-  select case(prec%iprcparm(psb_f_type_))
+  associate (wv => prec%wrk(1), wv1 => prec%wrk(2))
-  case(psb_f_ilu_n_) 
+    
-
+    select case(prec%iprcparm(psb_f_type_))
-    select case(trans_)
+    case(psb_f_ilu_n_) 
-    case('N')
+      
-      call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,&
+      select case(trans_)
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
+      case('N')
-      if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
+        call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('T')
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-      call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,&
+        
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
+      case('T')
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
+        call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('C')
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-
+        
-      call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,&
+      case('C')
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        
-
+        call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,&
-      call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_)
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-  
+        
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
+        call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_)
-           & beta,y,desc_data,info,&
+        
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
-
+             & beta,y,desc_data,info,&
-    end select
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-    if (info /= psb_success_) then 
+        
-      ch_err="psb_spsm"
+      end select
+      if (info /= psb_success_) then 
+        ch_err="psb_spsm"
+        goto 9999
+      end if
+      
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='Invalid factorization')
       goto 9999
-    end if
+    end select
-
-
-  case default
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='Invalid factorization')
-    goto 9999
-  end select
 
+  end associate
+  
   call psb_halo(y,desc_data,info,data=psb_comm_mov_)
 
-  call wv%free(info)
+  if (do_alloc_wrk) call prec%free_wrk(info)
-  call wv1%free(info)
   if (n_col <= size(work)) then 
     if ((4*n_col+n_col) <= size(work)) then 
     else

prec/impl/psb_z_bjacprec_impl.f90

@@ -90,6 +90,7 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
   integer(psb_ipk_) :: ictxt,np,me
   integer(psb_ipk_) :: err_act, ierr(5)
   integer(psb_ipk_) :: debug_level, debug_unit
+  logical            :: do_alloc_wrk
   character          :: trans_
   character(len=20)  :: name='z_bjac_prec_apply'
   character(len=20)  :: ch_err
@@ -154,55 +155,57 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
     goto 9999      
   end if
 
-  call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) 
+  do_alloc_wrk = .not.prec%is_allocated_wrk()
-  call psb_geasb(wv1,desc_data,info,mold=x%v,scratch=.true.) 
+  if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
-
+
-  select case(prec%iprcparm(psb_f_type_))
+  associate (wv => prec%wrk(1), wv1 => prec%wrk(2))
-  case(psb_f_ilu_n_) 
+    
-
+    select case(prec%iprcparm(psb_f_type_))
-    select case(trans_)
+    case(psb_f_ilu_n_) 
-    case('N')
+      
-      call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,&
+      select case(trans_)
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
+      case('N')
-      if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
+        call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_,work=aux)
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('T')
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-      call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,&
+        
-           & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
+      case('T')
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
+        call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,&
-           & beta,y,desc_data,info,&
+             & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux)
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv,&
-
+             & beta,y,desc_data,info,&
-    case('C')
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-
+        
-      call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,&
+      case('C')
-           & trans=trans_,scale='U',choice=psb_none_, work=aux)
+        
-
+        call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,&
-      call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_)
+             & trans=trans_,scale='U',choice=psb_none_, work=aux)
-  
+        
-      if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
+        call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_)
-           & beta,y,desc_data,info,&
+        
-           & trans=trans_,scale='U',choice=psb_none_,work=aux)
+        if(info == psb_success_)  call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
-
+             & beta,y,desc_data,info,&
-    end select
+             & trans=trans_,scale='U',choice=psb_none_,work=aux)
-    if (info /= psb_success_) then 
+        
-      ch_err="psb_spsm"
+      end select
+      if (info /= psb_success_) then 
+        ch_err="psb_spsm"
+        goto 9999
+      end if
+      
+    case default
+      info = psb_err_internal_error_
+      call psb_errpush(info,name,a_err='Invalid factorization')
       goto 9999
-    end if
+    end select
-
-
-  case default
-    info = psb_err_internal_error_
-    call psb_errpush(info,name,a_err='Invalid factorization')
-    goto 9999
-  end select
 
+  end associate
+  
   call psb_halo(y,desc_data,info,data=psb_comm_mov_)
 
-  call wv%free(info)
+  if (do_alloc_wrk) call prec%free_wrk(info)
-  call wv1%free(info)
   if (n_col <= size(work)) then 
     if ((4*n_col+n_col) <= size(work)) then 
     else

prec/psb_c_base_prec_mod.f90

@@ -61,6 +61,9 @@ module psb_c_base_prec_mod
     generic, public       :: build     => precbld
     generic, public       :: descr     => precdescr
     procedure, pass(prec) :: desc_prefix => psb_c_base_desc_prefix
+    procedure, pass(prec) :: allocate_wrk => psb_c_base_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_c_base_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_c_base_is_allocated_wrk
     procedure(psb_c_base_precbld), pass(prec), deferred :: precbld    
     procedure(psb_c_base_sizeof), pass(prec), deferred :: sizeof     
     procedure(psb_c_base_precinit), pass(prec), deferred :: precinit   
@@ -258,6 +261,87 @@ contains
 
   end subroutine psb_c_base_precsetc
 
+  subroutine psb_c_base_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_c_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_c_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_c_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_c_base_allocate_wrk
+
+  subroutine psb_c_base_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_c_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_c_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_c_base_free_wrk
+
+  function psb_c_base_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_c_base_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = .true.
+    
+  end function psb_c_base_is_allocated_wrk
+  
   subroutine psb_c_base_set_ctxt(prec,ictxt)
     implicit none 
     class(psb_c_base_prec_type), intent(inout) :: prec

prec/psb_c_bjacprec.f90

@@ -36,7 +36,7 @@ module psb_c_bjacprec
   type, extends(psb_c_base_prec_type)   :: psb_c_bjac_prec_type
     integer(psb_ipk_), allocatable      :: iprcparm(:)
     type(psb_cspmat_type), allocatable  :: av(:)
-    type(psb_c_vect_type), allocatable  :: dv
+    type(psb_c_vect_type), allocatable  :: dv, wrk(:)
   contains
     procedure, pass(prec) :: c_apply_v => psb_c_bjac_apply_vect
     procedure, pass(prec) :: c_apply   => psb_c_bjac_apply
@@ -49,6 +49,9 @@ module psb_c_bjacprec
     procedure, pass(prec) :: free      => psb_c_bjac_precfree
     procedure, pass(prec) :: sizeof    => psb_c_bjac_sizeof
     procedure, pass(prec) :: get_nzeros => psb_c_bjac_get_nzeros
+    procedure, pass(prec) :: allocate_wrk => psb_c_bjac_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_c_bjac_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_c_bjac_is_allocated_wrk
   end type psb_c_bjac_prec_type
 
   private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros
@@ -308,4 +311,111 @@ contains
 
   end subroutine psb_c_bjac_clone
 
+  subroutine psb_c_bjac_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_c_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_c_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act, i
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_c_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+    if (allocated(prec%wrk)) then
+      if (size(prec%wrk)<2) then
+        do i=1,size(prec%wrk)
+          if (info == 0) call prec%wrk(i)%free(info)
+        end do
+        if (info == 0) deallocate(prec%wrk,stat=info)
+      end if
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    if (.not.allocated(prec%wrk)) then
+      if (.not.present(desc)) then 
+        info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?");   goto 9999
+      end if      
+      allocate(prec%wrk(2),stat=info)
+      do i=1, 2
+        if (info == 0)   call psb_geall(prec%wrk(i),desc,info)
+        if (info == 0)   call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.)
+      end do
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_c_bjac_allocate_wrk
+
+  subroutine psb_c_bjac_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_c_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    integer(psb_ipk_) :: i
+    character(len=20) :: name
+    
+    info=psb_success_
+    name = 'psb_c_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    info = psb_success_ 
+    if (allocated(prec%wrk)) then
+      do i=1,size(prec%wrk)
+        if (info == 0) call prec%wrk(i)%free(info)
+      end do
+      if (info == 0) deallocate(prec%wrk,stat=info)
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_c_bjac_free_wrk
+
+  function psb_c_bjac_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_c_bjac_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = allocated(prec%wrk)
+    
+  end function psb_c_bjac_is_allocated_wrk
+
+  
 end module psb_c_bjacprec

prec/psb_c_prec_type.f90

@@ -54,6 +54,9 @@ module psb_c_prec_type
     procedure, pass(prec)               :: build  => psb_cprecbld
     procedure, pass(prec)               :: init   => psb_cprecinit
     procedure, pass(prec)               :: descr  => psb_cfile_prec_descr
+    procedure, pass(prec)               :: allocate_wrk => psb_c_allocate_wrk
+    procedure, pass(prec)               :: free_wrk => psb_c_free_wrk
+    procedure, pass(prec)               :: is_allocated_wrk => psb_c_is_allocated_wrk
   end type psb_cprec_type
 
   interface psb_precfree
@@ -193,6 +196,90 @@ contains
     
   end subroutine psb_c_prec_dump
 
+  subroutine psb_c_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_cprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_c_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_c_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_c_allocate_wrk
+
+  subroutine psb_c_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_cprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_c_free_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to free a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%free_wrk(info)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_c_free_wrk
+
+  function psb_c_is_allocated_wrk(prec) result(res)
+    implicit none
+    
+    ! Arguments
+    class(psb_cprec_type), intent(in) :: prec
+    logical :: res
+
+    if (.not.allocated(prec%prec)) then
+      res = .false.
+    else 
+      res = prec%prec%is_allocated_wrk()
+    end if
+    
+  end function psb_c_is_allocated_wrk
 
   subroutine psb_c_precfree(p,info)
     use psb_base_mod

prec/psb_d_base_prec_mod.f90

@@ -61,6 +61,9 @@ module psb_d_base_prec_mod
     generic, public       :: build     => precbld
     generic, public       :: descr     => precdescr
     procedure, pass(prec) :: desc_prefix => psb_d_base_desc_prefix
+    procedure, pass(prec) :: allocate_wrk => psb_d_base_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_d_base_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_d_base_is_allocated_wrk
     procedure(psb_d_base_precbld), pass(prec), deferred :: precbld    
     procedure(psb_d_base_sizeof), pass(prec), deferred :: sizeof     
     procedure(psb_d_base_precinit), pass(prec), deferred :: precinit   
@@ -258,6 +261,87 @@ contains
 
   end subroutine psb_d_base_precsetc
 
+  subroutine psb_d_base_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_d_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_d_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_d_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_d_base_allocate_wrk
+
+  subroutine psb_d_base_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_d_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_d_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_d_base_free_wrk
+
+  function psb_d_base_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_d_base_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = .true.
+    
+  end function psb_d_base_is_allocated_wrk
+  
   subroutine psb_d_base_set_ctxt(prec,ictxt)
     implicit none 
     class(psb_d_base_prec_type), intent(inout) :: prec

prec/psb_d_bjacprec.f90

@@ -36,7 +36,7 @@ module psb_d_bjacprec
   type, extends(psb_d_base_prec_type)   :: psb_d_bjac_prec_type
     integer(psb_ipk_), allocatable      :: iprcparm(:)
     type(psb_dspmat_type), allocatable  :: av(:)
-    type(psb_d_vect_type), allocatable  :: dv
+    type(psb_d_vect_type), allocatable  :: dv, wrk(:)
   contains
     procedure, pass(prec) :: d_apply_v => psb_d_bjac_apply_vect
     procedure, pass(prec) :: d_apply   => psb_d_bjac_apply
@@ -49,6 +49,9 @@ module psb_d_bjacprec
     procedure, pass(prec) :: free      => psb_d_bjac_precfree
     procedure, pass(prec) :: sizeof    => psb_d_bjac_sizeof
     procedure, pass(prec) :: get_nzeros => psb_d_bjac_get_nzeros
+    procedure, pass(prec) :: allocate_wrk => psb_d_bjac_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_d_bjac_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_d_bjac_is_allocated_wrk
   end type psb_d_bjac_prec_type
 
   private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros
@@ -308,4 +311,111 @@ contains
 
   end subroutine psb_d_bjac_clone
 
+  subroutine psb_d_bjac_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_d_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_d_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act, i
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_d_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+    if (allocated(prec%wrk)) then
+      if (size(prec%wrk)<2) then
+        do i=1,size(prec%wrk)
+          if (info == 0) call prec%wrk(i)%free(info)
+        end do
+        if (info == 0) deallocate(prec%wrk,stat=info)
+      end if
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    if (.not.allocated(prec%wrk)) then
+      if (.not.present(desc)) then 
+        info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?");   goto 9999
+      end if      
+      allocate(prec%wrk(2),stat=info)
+      do i=1, 2
+        if (info == 0)   call psb_geall(prec%wrk(i),desc,info)
+        if (info == 0)   call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.)
+      end do
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_d_bjac_allocate_wrk
+
+  subroutine psb_d_bjac_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_d_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    integer(psb_ipk_) :: i
+    character(len=20) :: name
+    
+    info=psb_success_
+    name = 'psb_d_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    info = psb_success_ 
+    if (allocated(prec%wrk)) then
+      do i=1,size(prec%wrk)
+        if (info == 0) call prec%wrk(i)%free(info)
+      end do
+      if (info == 0) deallocate(prec%wrk,stat=info)
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_d_bjac_free_wrk
+
+  function psb_d_bjac_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_d_bjac_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = allocated(prec%wrk)
+    
+  end function psb_d_bjac_is_allocated_wrk
+
+  
 end module psb_d_bjacprec

prec/psb_d_prec_type.f90

@@ -54,6 +54,9 @@ module psb_d_prec_type
     procedure, pass(prec)               :: build  => psb_dprecbld
     procedure, pass(prec)               :: init   => psb_dprecinit
     procedure, pass(prec)               :: descr  => psb_dfile_prec_descr
+    procedure, pass(prec)               :: allocate_wrk => psb_d_allocate_wrk
+    procedure, pass(prec)               :: free_wrk => psb_d_free_wrk
+    procedure, pass(prec)               :: is_allocated_wrk => psb_d_is_allocated_wrk
   end type psb_dprec_type
 
   interface psb_precfree
@@ -193,6 +196,90 @@ contains
     
   end subroutine psb_d_prec_dump
 
+  subroutine psb_d_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_dprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_d_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_d_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_d_allocate_wrk
+
+  subroutine psb_d_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_dprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_d_free_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to free a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%free_wrk(info)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_d_free_wrk
+
+  function psb_d_is_allocated_wrk(prec) result(res)
+    implicit none
+    
+    ! Arguments
+    class(psb_dprec_type), intent(in) :: prec
+    logical :: res
+
+    if (.not.allocated(prec%prec)) then
+      res = .false.
+    else 
+      res = prec%prec%is_allocated_wrk()
+    end if
+    
+  end function psb_d_is_allocated_wrk
 
   subroutine psb_d_precfree(p,info)
     use psb_base_mod

prec/psb_s_base_prec_mod.f90

@@ -61,6 +61,9 @@ module psb_s_base_prec_mod
     generic, public       :: build     => precbld
     generic, public       :: descr     => precdescr
     procedure, pass(prec) :: desc_prefix => psb_s_base_desc_prefix
+    procedure, pass(prec) :: allocate_wrk => psb_s_base_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_s_base_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_s_base_is_allocated_wrk
     procedure(psb_s_base_precbld), pass(prec), deferred :: precbld    
     procedure(psb_s_base_sizeof), pass(prec), deferred :: sizeof     
     procedure(psb_s_base_precinit), pass(prec), deferred :: precinit   
@@ -258,6 +261,87 @@ contains
 
   end subroutine psb_s_base_precsetc
 
+  subroutine psb_s_base_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_s_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_s_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_s_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_s_base_allocate_wrk
+
+  subroutine psb_s_base_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_s_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_s_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_s_base_free_wrk
+
+  function psb_s_base_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_s_base_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = .true.
+    
+  end function psb_s_base_is_allocated_wrk
+  
   subroutine psb_s_base_set_ctxt(prec,ictxt)
     implicit none 
     class(psb_s_base_prec_type), intent(inout) :: prec

prec/psb_s_bjacprec.f90

@@ -36,7 +36,7 @@ module psb_s_bjacprec
   type, extends(psb_s_base_prec_type)   :: psb_s_bjac_prec_type
     integer(psb_ipk_), allocatable      :: iprcparm(:)
     type(psb_sspmat_type), allocatable  :: av(:)
-    type(psb_s_vect_type), allocatable  :: dv
+    type(psb_s_vect_type), allocatable  :: dv, wrk(:)
   contains
     procedure, pass(prec) :: s_apply_v => psb_s_bjac_apply_vect
     procedure, pass(prec) :: s_apply   => psb_s_bjac_apply
@@ -49,6 +49,9 @@ module psb_s_bjacprec
     procedure, pass(prec) :: free      => psb_s_bjac_precfree
     procedure, pass(prec) :: sizeof    => psb_s_bjac_sizeof
     procedure, pass(prec) :: get_nzeros => psb_s_bjac_get_nzeros
+    procedure, pass(prec) :: allocate_wrk => psb_s_bjac_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_s_bjac_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_s_bjac_is_allocated_wrk
   end type psb_s_bjac_prec_type
 
   private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros
@@ -308,4 +311,111 @@ contains
 
   end subroutine psb_s_bjac_clone
 
+  subroutine psb_s_bjac_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_s_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_s_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act, i
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_s_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+    if (allocated(prec%wrk)) then
+      if (size(prec%wrk)<2) then
+        do i=1,size(prec%wrk)
+          if (info == 0) call prec%wrk(i)%free(info)
+        end do
+        if (info == 0) deallocate(prec%wrk,stat=info)
+      end if
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    if (.not.allocated(prec%wrk)) then
+      if (.not.present(desc)) then 
+        info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?");   goto 9999
+      end if      
+      allocate(prec%wrk(2),stat=info)
+      do i=1, 2
+        if (info == 0)   call psb_geall(prec%wrk(i),desc,info)
+        if (info == 0)   call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.)
+      end do
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_s_bjac_allocate_wrk
+
+  subroutine psb_s_bjac_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_s_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    integer(psb_ipk_) :: i
+    character(len=20) :: name
+    
+    info=psb_success_
+    name = 'psb_s_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    info = psb_success_ 
+    if (allocated(prec%wrk)) then
+      do i=1,size(prec%wrk)
+        if (info == 0) call prec%wrk(i)%free(info)
+      end do
+      if (info == 0) deallocate(prec%wrk,stat=info)
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_s_bjac_free_wrk
+
+  function psb_s_bjac_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_s_bjac_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = allocated(prec%wrk)
+    
+  end function psb_s_bjac_is_allocated_wrk
+
+  
 end module psb_s_bjacprec

prec/psb_s_prec_type.f90

@@ -54,6 +54,9 @@ module psb_s_prec_type
     procedure, pass(prec)               :: build  => psb_sprecbld
     procedure, pass(prec)               :: init   => psb_sprecinit
     procedure, pass(prec)               :: descr  => psb_sfile_prec_descr
+    procedure, pass(prec)               :: allocate_wrk => psb_s_allocate_wrk
+    procedure, pass(prec)               :: free_wrk => psb_s_free_wrk
+    procedure, pass(prec)               :: is_allocated_wrk => psb_s_is_allocated_wrk
   end type psb_sprec_type
 
   interface psb_precfree
@@ -193,6 +196,90 @@ contains
     
   end subroutine psb_s_prec_dump
 
+  subroutine psb_s_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_sprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_s_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_s_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_s_allocate_wrk
+
+  subroutine psb_s_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_sprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_s_free_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to free a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%free_wrk(info)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_s_free_wrk
+
+  function psb_s_is_allocated_wrk(prec) result(res)
+    implicit none
+    
+    ! Arguments
+    class(psb_sprec_type), intent(in) :: prec
+    logical :: res
+
+    if (.not.allocated(prec%prec)) then
+      res = .false.
+    else 
+      res = prec%prec%is_allocated_wrk()
+    end if
+    
+  end function psb_s_is_allocated_wrk
 
   subroutine psb_s_precfree(p,info)
     use psb_base_mod

prec/psb_z_base_prec_mod.f90

@@ -61,6 +61,9 @@ module psb_z_base_prec_mod
     generic, public       :: build     => precbld
     generic, public       :: descr     => precdescr
     procedure, pass(prec) :: desc_prefix => psb_z_base_desc_prefix
+    procedure, pass(prec) :: allocate_wrk => psb_z_base_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_z_base_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_z_base_is_allocated_wrk
     procedure(psb_z_base_precbld), pass(prec), deferred :: precbld    
     procedure(psb_z_base_sizeof), pass(prec), deferred :: sizeof     
     procedure(psb_z_base_precinit), pass(prec), deferred :: precinit   
@@ -258,6 +261,87 @@ contains
 
   end subroutine psb_z_base_precsetc
 
+  subroutine psb_z_base_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_z_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_z_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_z_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_z_base_allocate_wrk
+
+  subroutine psb_z_base_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_z_base_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_z_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    !
+    ! Base version does nothing.
+    !
+
+    info = psb_success_ 
+
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_z_base_free_wrk
+
+  function psb_z_base_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_z_base_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = .true.
+    
+  end function psb_z_base_is_allocated_wrk
+  
   subroutine psb_z_base_set_ctxt(prec,ictxt)
     implicit none 
     class(psb_z_base_prec_type), intent(inout) :: prec

prec/psb_z_bjacprec.f90

@@ -36,7 +36,7 @@ module psb_z_bjacprec
   type, extends(psb_z_base_prec_type)   :: psb_z_bjac_prec_type
     integer(psb_ipk_), allocatable      :: iprcparm(:)
     type(psb_zspmat_type), allocatable  :: av(:)
-    type(psb_z_vect_type), allocatable  :: dv
+    type(psb_z_vect_type), allocatable  :: dv, wrk(:)
   contains
     procedure, pass(prec) :: z_apply_v => psb_z_bjac_apply_vect
     procedure, pass(prec) :: z_apply   => psb_z_bjac_apply
@@ -49,6 +49,9 @@ module psb_z_bjacprec
     procedure, pass(prec) :: free      => psb_z_bjac_precfree
     procedure, pass(prec) :: sizeof    => psb_z_bjac_sizeof
     procedure, pass(prec) :: get_nzeros => psb_z_bjac_get_nzeros
+    procedure, pass(prec) :: allocate_wrk => psb_z_bjac_allocate_wrk
+    procedure, pass(prec) :: free_wrk     => psb_z_bjac_free_wrk
+    procedure, pass(prec) :: is_allocated_wrk => psb_z_bjac_is_allocated_wrk
   end type psb_z_bjac_prec_type
 
   private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros
@@ -308,4 +311,111 @@ contains
 
   end subroutine psb_z_bjac_clone
 
+  subroutine psb_z_bjac_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_z_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_z_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act, i
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_z_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+    if (allocated(prec%wrk)) then
+      if (size(prec%wrk)<2) then
+        do i=1,size(prec%wrk)
+          if (info == 0) call prec%wrk(i)%free(info)
+        end do
+        if (info == 0) deallocate(prec%wrk,stat=info)
+      end if
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    if (.not.allocated(prec%wrk)) then
+      if (.not.present(desc)) then 
+        info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?");   goto 9999
+      end if      
+      allocate(prec%wrk(2),stat=info)
+      do i=1, 2
+        if (info == 0)   call psb_geall(prec%wrk(i),desc,info)
+        if (info == 0)   call psb_geasb(prec%wrk(i),desc,info,mold=vmold,scratch=.true.)
+      end do
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_z_bjac_allocate_wrk
+
+  subroutine psb_z_bjac_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_z_bjac_prec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    integer(psb_ipk_) :: i
+    character(len=20) :: name
+    
+    info=psb_success_
+    name = 'psb_z_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    info = psb_success_ 
+    if (allocated(prec%wrk)) then
+      do i=1,size(prec%wrk)
+        if (info == 0) call prec%wrk(i)%free(info)
+      end do
+      if (info == 0) deallocate(prec%wrk,stat=info)
+    end if
+    if (info /= 0) then
+      info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate");   goto 9999
+    end if
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_z_bjac_free_wrk
+
+  function psb_z_bjac_is_allocated_wrk(prec) result(res)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_z_bjac_prec_type), intent(in) :: prec
+    logical :: res
+
+    ! In the base version we can say yes, because 
+    ! there is nothing to allocate
+
+    res = allocated(prec%wrk)
+    
+  end function psb_z_bjac_is_allocated_wrk
+
+  
 end module psb_z_bjacprec

prec/psb_z_prec_type.f90

@@ -54,6 +54,9 @@ module psb_z_prec_type
     procedure, pass(prec)               :: build  => psb_zprecbld
     procedure, pass(prec)               :: init   => psb_zprecinit
     procedure, pass(prec)               :: descr  => psb_zfile_prec_descr
+    procedure, pass(prec)               :: allocate_wrk => psb_z_allocate_wrk
+    procedure, pass(prec)               :: free_wrk => psb_z_free_wrk
+    procedure, pass(prec)               :: is_allocated_wrk => psb_z_is_allocated_wrk
   end type psb_zprec_type
 
   interface psb_precfree
@@ -193,6 +196,90 @@ contains
     
   end subroutine psb_z_prec_dump
 
+  subroutine psb_z_allocate_wrk(prec,info,vmold,desc)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_zprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_z_base_vect_type), intent(in), optional  :: vmold
+    type(psb_desc_type), intent(in), optional :: desc
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_z_allocate_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to allocate wrk to a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%allocate_wrk(info,vmold=vmold,desc=desc)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_z_allocate_wrk
+
+  subroutine psb_z_free_wrk(prec,info)
+    use psb_base_mod
+    implicit none
+    
+    ! Arguments
+    class(psb_zprec_type), intent(inout) :: prec
+    integer(psb_ipk_), intent(out)        :: info
+
+    ! Local variables
+    integer(psb_ipk_) :: err_act
+    character(len=20)   :: name
+    
+    info=psb_success_
+    name = 'psb_z_free_wrk'
+    call psb_erractionsave(err_act)
+    
+    if (psb_get_errstatus().ne.0) goto 9999
+
+    if (.not.allocated(prec%prec)) then 
+      info = -1
+      write(psb_err_unit,*) 'Trying to free a non-built preconditioner'
+      return
+    end if
+    
+    call prec%prec%free_wrk(info)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 call psb_error_handler(err_act)
+    return
+    
+  end subroutine psb_z_free_wrk
+
+  function psb_z_is_allocated_wrk(prec) result(res)
+    implicit none
+    
+    ! Arguments
+    class(psb_zprec_type), intent(in) :: prec
+    logical :: res
+
+    if (.not.allocated(prec%prec)) then
+      res = .false.
+    else 
+      res = prec%prec%is_allocated_wrk()
+    end if
+    
+  end function psb_z_is_allocated_wrk
 
   subroutine psb_z_precfree(p,info)
     use psb_base_mod