Corrected bug in application phase

implement-ainv
Cirdans-Home 4 years ago
parent 3e22e9e963
commit 9d1390ad94

@ -13,7 +13,7 @@ MODOBJS=psb_prec_const_mod.o\
psb_s_diagprec.o psb_s_nullprec.o psb_s_bjacprec.o psb_d_ilu_fact_mod.o \
psb_c_diagprec.o psb_c_nullprec.o psb_c_bjacprec.o psb_c_ilu_fact_mod.o \
psb_z_diagprec.o psb_z_nullprec.o psb_z_bjacprec.o psb_z_ilu_fact_mod.o \
psb_base_ainv_mod.o \
psb_c_ainv_fact_mod.o psb_d_ainv_fact_mod.o psb_s_ainv_fact_mod.o psb_z_ainv_fact_mod.o \
psb_c_ainv_tools_mod.o psb_d_ainv_tools_mod.o psb_s_ainv_tools_mod.o psb_z_ainv_tools_mod.o \
psb_ainv_tools_mod.o \
psb_biconjg_mod.o psb_c_biconjg_mod.o psb_d_biconjg_mod.o psb_s_biconjg_mod.o \
@ -56,7 +56,10 @@ psb_s_bjacprec.o: psb_s_ilu_fact_mod.o psb_s_ainv_fact_mod.o
psb_d_bjacprec.o: psb_d_ilu_fact_mod.o psb_d_ainv_fact_mod.o
psb_c_bjacprec.o: psb_c_ilu_fact_mod.o psb_c_ainv_fact_mod.o
psb_z_bjacprec.o: psb_z_ilu_fact_mod.o psb_z_ainv_fact_mod.o
psb_c_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o psb_s_ainv_fact_mod.o psb_d_ainv_fact_mod.o psb_c_ainv_fact_mod.o psb_z_ainv_fact_mod.o
psb_d_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o
psb_s_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o
psb_c_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o
psb_z_ainv_fact_mod.o: psb_prec_const_mod.o psb_ainv_tools_mod.o
psb_ainv_tools_mod.o: psb_c_ainv_tools_mod.o psb_d_ainv_tools_mod.o psb_s_ainv_tools_mod.o psb_z_ainv_tools_mod.o
psb_biconjg_mod.o: psb_prec_const_mod.o psb_c_biconjg_mod.o \
psb_d_biconjg_mod.o psb_s_biconjg_mod.o psb_z_biconjg_mod.o

@ -200,34 +200,20 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N')
call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,&
& trans=trans_, work=aux)
& trans=trans_,work=aux,doswap=.false.)
call wv1%mlt(cone,prec%dv,wv,czero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
& beta,y,desc_data,info,&
& trans=trans_, work=aux)
case('T')
call psb_spmm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,&
& trans=trans_, work=aux)
call wv1%mlt(cone,prec%dv,wv,czero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('C')
call psb_spmm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,&
& trans=trans_,work=aux)
call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_)
if (info == psb_success_) call wv1%mlt(cone,prec%dv,wv,czero,info)
if(info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1,&
& beta,y,desc_data,info, trans=trans_, work=aux,doswap=.false.)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('T','C')
call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
if (info == psb_success_) call wv1%mlt(cone,prec%dv,wv,czero,info)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1, &
& beta,y,desc_data,info,trans=trans_,work=aux,doswap=.false.)
end select
if (info /= psb_success_) then
@ -394,19 +380,30 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N','T')
call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('C')
call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*conjg(prec%dv%v%v(1:n_row))
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('N')
call psb_spmm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('T')
call psb_spmm(cone,prec%av(psb_u_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('C')
call psb_spmm(cone,prec%av(psb_u_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * conjg(prec%dv%v%v(1:n_row))
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
end select
@ -592,6 +589,8 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold)
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
end if
end if
inv_fill = prec%iprcparm(psb_inv_fillin_)
if (inv_fill <= 0) inv_fill = m ! If no limit on the fill_in is required we allow everything
! Select on the type of factorization to be used
select case(prec%iprcparm(psb_f_type_))

@ -453,6 +453,14 @@ subroutine psb_ccprecsetc(prec,what,string,info,ilev,ilmax,pos,idx)
select case (psb_toupper(string))
case ("MAXVAL")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info)
case ("DIAG")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_diag_,info)
case ("ARWSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arwsum_,info)
case ("ARCSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arcsum_,info)
case ("ACLSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_aclsum_,info)
case default
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info)
end select

@ -100,7 +100,6 @@ subroutine psb_c_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe
end do
else
allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info)
if (info /= psb_success_) then
return

@ -224,7 +224,6 @@ subroutine psb_csparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info)
!
! Sparsify current ZVAL and put into ZMAT
!
write(psb_out_unit,'("z(1) = ",f16.14)') zval(1)
call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr)
if (info /= psb_success_) then
info = psb_err_internal_error_

@ -200,34 +200,20 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N')
call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,&
& trans=trans_, work=aux)
& trans=trans_,work=aux,doswap=.false.)
call wv1%mlt(done,prec%dv,wv,dzero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
& beta,y,desc_data,info,&
& trans=trans_, work=aux)
case('T')
call psb_spmm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,&
& trans=trans_, work=aux)
call wv1%mlt(done,prec%dv,wv,dzero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('C')
call psb_spmm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,&
& trans=trans_,work=aux)
call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_)
if (info == psb_success_) call wv1%mlt(done,prec%dv,wv,dzero,info)
if(info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1,&
& beta,y,desc_data,info, trans=trans_, work=aux,doswap=.false.)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('T','C')
call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
if (info == psb_success_) call wv1%mlt(done,prec%dv,wv,dzero,info)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1, &
& beta,y,desc_data,info,trans=trans_,work=aux,doswap=.false.)
end select
if (info /= psb_success_) then
@ -394,19 +380,30 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N','T')
call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('C')
call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*(prec%dv%v%v(1:n_row))
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('N')
call psb_spmm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('T')
call psb_spmm(done,prec%av(psb_u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('C')
call psb_spmm(done,prec%av(psb_u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * (prec%dv%v%v(1:n_row))
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
end select
@ -592,6 +589,8 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold)
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
end if
end if
inv_fill = prec%iprcparm(psb_inv_fillin_)
if (inv_fill <= 0) inv_fill = m ! If no limit on the fill_in is required we allow everything
! Select on the type of factorization to be used
select case(prec%iprcparm(psb_f_type_))

@ -453,6 +453,14 @@ subroutine psb_dcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx)
select case (psb_toupper(string))
case ("MAXVAL")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info)
case ("DIAG")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_diag_,info)
case ("ARWSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arwsum_,info)
case ("ARCSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arcsum_,info)
case ("ACLSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_aclsum_,info)
case default
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info)
end select

@ -100,7 +100,6 @@ subroutine psb_d_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe
end do
else
allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info)
if (info /= psb_success_) then
return

@ -224,7 +224,6 @@ subroutine psb_dsparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info)
!
! Sparsify current ZVAL and put into ZMAT
!
write(psb_out_unit,'("z(1) = ",f16.14)') zval(1)
call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr)
if (info /= psb_success_) then
info = psb_err_internal_error_

@ -200,34 +200,20 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N')
call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,&
& trans=trans_, work=aux)
& trans=trans_,work=aux,doswap=.false.)
call wv1%mlt(sone,prec%dv,wv,szero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
& beta,y,desc_data,info,&
& trans=trans_, work=aux)
case('T')
call psb_spmm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,&
& trans=trans_, work=aux)
call wv1%mlt(sone,prec%dv,wv,szero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('C')
call psb_spmm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,&
& trans=trans_,work=aux)
call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_)
if (info == psb_success_) call wv1%mlt(sone,prec%dv,wv,szero,info)
if(info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1,&
& beta,y,desc_data,info, trans=trans_, work=aux,doswap=.false.)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('T','C')
call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
if (info == psb_success_) call wv1%mlt(sone,prec%dv,wv,szero,info)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1, &
& beta,y,desc_data,info,trans=trans_,work=aux,doswap=.false.)
end select
if (info /= psb_success_) then
@ -394,19 +380,30 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N','T')
call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('C')
call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*(prec%dv%v%v(1:n_row))
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('N')
call psb_spmm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('T')
call psb_spmm(sone,prec%av(psb_u_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('C')
call psb_spmm(sone,prec%av(psb_u_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * (prec%dv%v%v(1:n_row))
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
end select
@ -592,6 +589,8 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold)
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
end if
end if
inv_fill = prec%iprcparm(psb_inv_fillin_)
if (inv_fill <= 0) inv_fill = m ! If no limit on the fill_in is required we allow everything
! Select on the type of factorization to be used
select case(prec%iprcparm(psb_f_type_))

@ -453,6 +453,14 @@ subroutine psb_scprecsetc(prec,what,string,info,ilev,ilmax,pos,idx)
select case (psb_toupper(string))
case ("MAXVAL")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info)
case ("DIAG")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_diag_,info)
case ("ARWSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arwsum_,info)
case ("ARCSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arcsum_,info)
case ("ACLSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_aclsum_,info)
case default
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info)
end select

@ -100,7 +100,6 @@ subroutine psb_s_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe
end do
else
allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info)
if (info /= psb_success_) then
return

@ -224,7 +224,6 @@ subroutine psb_ssparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info)
!
! Sparsify current ZVAL and put into ZMAT
!
write(psb_out_unit,'("z(1) = ",f16.14)') zval(1)
call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr)
if (info /= psb_success_) then
info = psb_err_internal_error_

@ -200,34 +200,20 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N')
call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,&
& trans=trans_, work=aux)
& trans=trans_,work=aux,doswap=.false.)
call wv1%mlt(zone,prec%dv,wv,zzero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,&
& beta,y,desc_data,info,&
& trans=trans_, work=aux)
case('T')
call psb_spmm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,&
& trans=trans_, work=aux)
call wv1%mlt(zone,prec%dv,wv,zzero,info)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('C')
call psb_spmm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,&
& trans=trans_,work=aux)
call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_)
if (info == psb_success_) call wv1%mlt(zone,prec%dv,wv,zzero,info)
if(info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1,&
& beta,y,desc_data,info, trans=trans_, work=aux,doswap=.false.)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,&
& beta,y,desc_data,info,&
& trans=trans_,work=aux)
case('T','C')
call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
if (info == psb_success_) call wv1%mlt(zone,prec%dv,wv,zzero,info)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),wv1, &
& beta,y,desc_data,info,trans=trans_,work=aux,doswap=.false.)
end select
if (info /= psb_success_) then
@ -394,19 +380,30 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
select case(trans_)
case('N','T')
call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*prec%dv%v%v(1:n_row)
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('C')
call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_, work=aux)
ww(1:n_row) = ww(1:n_row)*conjg(prec%dv%v%v(1:n_row))
if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,&
& beta,y,desc_data,info, trans=trans_, work=aux)
case('N')
call psb_spmm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('T')
call psb_spmm(zone,prec%av(psb_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * prec%dv%v%v(1:n_row)
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
case('C')
call psb_spmm(zone,prec%av(psb_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
ww(1:n_row) = ww(1:n_row) * conjg(prec%dv%v%v(1:n_row))
if (info == psb_success_) &
& call psb_spmm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,work=aux,doswap=.false.)
end select
@ -592,6 +589,8 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold)
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
end if
end if
inv_fill = prec%iprcparm(psb_inv_fillin_)
if (inv_fill <= 0) inv_fill = m ! If no limit on the fill_in is required we allow everything
! Select on the type of factorization to be used
select case(prec%iprcparm(psb_f_type_))

@ -453,6 +453,14 @@ subroutine psb_zcprecsetc(prec,what,string,info,ilev,ilmax,pos,idx)
select case (psb_toupper(string))
case ("MAXVAL")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_maxval_,info)
case ("DIAG")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_diag_,info)
case ("ARWSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arwsum_,info)
case ("ARCSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_arcsum_,info)
case ("ACLSUM")
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_aclsum_,info)
case default
call prec%prec%precset(psb_ilu_scale_,psb_ilu_scale_none_,info)
end select

@ -100,7 +100,6 @@ subroutine psb_z_sparsify(idiag,nzrmax,sp_thresh,n,zw,nz,iz,valz,info,istart,ihe
end do
else
allocate(xw(nzrmax),xwid(nzrmax),indx(nzrmax),stat=info)
if (info /= psb_success_) then
return

@ -224,7 +224,6 @@ subroutine psb_zsparse_biconjg_llk(n,a,p,z,w,nzrmax,sp_thresh,info)
!
! Sparsify current ZVAL and put into ZMAT
!
write(psb_out_unit,'("z(1) = ",f16.14)') zval(1)
call sparsify(i,nzrmax,sp_thresh,n,zval,nzrz,ia,val,info,iheap=heap,ikr=izkr)
if (info /= psb_success_) then
info = psb_err_internal_error_

@ -166,6 +166,8 @@ contains
goto 9999
end if
if (nzrmax <= 0) write(psb_out_unit,'("Out nzrmax = ",i0)') nzrmax
select case(alg)
case (psb_ainv_llk_)
call psb_csparse_biconjg_llk(n,acsr,p,zcsc,wcsc,nzrmax,sp_thresh,info)

@ -166,6 +166,8 @@ contains
goto 9999
end if
if (nzrmax <= 0) write(psb_out_unit,'("Out nzrmax = ",i0)') nzrmax
select case(alg)
case (psb_ainv_llk_)
call psb_dsparse_biconjg_llk(n,acsr,p,zcsc,wcsc,nzrmax,sp_thresh,info)

@ -166,6 +166,8 @@ contains
goto 9999
end if
if (nzrmax <= 0) write(psb_out_unit,'("Out nzrmax = ",i0)') nzrmax
select case(alg)
case (psb_ainv_llk_)
call psb_ssparse_biconjg_llk(n,acsr,p,zcsc,wcsc,nzrmax,sp_thresh,info)

@ -166,6 +166,8 @@ contains
goto 9999
end if
if (nzrmax <= 0) write(psb_out_unit,'("Out nzrmax = ",i0)') nzrmax
select case(alg)
case (psb_ainv_llk_)
call psb_zsparse_biconjg_llk(n,acsr,p,zcsc,wcsc,nzrmax,sp_thresh,info)

@ -678,7 +678,9 @@ program psb_d_pde3d
call prec%set('sub_iluthrs', parms%thresh, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case ("AINV")
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case default
! Do nothing, use default setting in the init routine
end select
@ -888,7 +890,9 @@ contains
write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh
case ('AINV','AORTH')
write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh
write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill
write(psb_out_unit,'("Orthogonalization : ",a)') parms%orth_alg
write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale
case default
write(psb_out_unit,'("Unknown diagonal solver")')
end select

@ -678,7 +678,9 @@ program psb_s_pde3d
call prec%set('sub_iluthrs', parms%thresh, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case ("AINV")
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case default
! Do nothing, use default setting in the init routine
end select
@ -888,7 +890,9 @@ contains
write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh
case ('AINV','AORTH')
write(psb_out_unit,'("Inverse Threshold : ",es12.5)') parms%inv_thresh
write(psb_out_unit,'("Invese Fill in : ",i0)') parms%inv_fill
write(psb_out_unit,'("Orthogonalization : ",a)') parms%orth_alg
write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale
case default
write(psb_out_unit,'("Unknown diagonal solver")')
end select

Loading…
Cancel
Save