From 0c9098065a83eb02005facc1d1c0619973945189 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Fri, 20 Nov 2020 14:08:50 +0100 Subject: [PATCH] Integrated INVT and INVK in BJAC prec --- prec/impl/psb_c_bjacprec_impl.f90 | 155 ++++++++++++++++++++++++++++-- prec/impl/psb_d_bjacprec_impl.f90 | 155 ++++++++++++++++++++++++++++-- prec/impl/psb_s_bjacprec_impl.f90 | 155 ++++++++++++++++++++++++++++-- prec/impl/psb_z_bjacprec_impl.f90 | 155 ++++++++++++++++++++++++++++-- prec/psb_c_bjacprec.f90 | 6 +- prec/psb_d_bjacprec.f90 | 6 +- prec/psb_prec_const_mod.f90 | 4 +- prec/psb_s_bjacprec.f90 | 6 +- prec/psb_z_bjacprec.f90 | 6 +- 9 files changed, 614 insertions(+), 34 deletions(-) diff --git a/prec/impl/psb_c_bjacprec_impl.f90 b/prec/impl/psb_c_bjacprec_impl.f90 index 315504da..45b96432 100644 --- a/prec/impl/psb_c_bjacprec_impl.f90 +++ b/prec/impl/psb_c_bjacprec_impl.f90 @@ -195,7 +195,8 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) + ! Application of approximate inverse preconditioner, just some spmm select case(trans_) case('N') @@ -373,11 +374,8 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) ! Application of approximate inverse preconditioner, just some spmm - ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & - ! & trans, work) - select case(trans_) @@ -866,7 +864,7 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 endif - ! Computin the factorization + ! Computing the factorization call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) if(info == psb_success_) then @@ -885,6 +883,151 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_invk_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(n) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invk_fact(a,fill_in, inv_fill,lf,dd,uf,desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_invt_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(eps) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invt_fact(a,fill_in,inv_fill,fact_eps,inv_thresh,lf,dd,uf,& + & desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' diff --git a/prec/impl/psb_d_bjacprec_impl.f90 b/prec/impl/psb_d_bjacprec_impl.f90 index d2b5daa4..f07e713e 100644 --- a/prec/impl/psb_d_bjacprec_impl.f90 +++ b/prec/impl/psb_d_bjacprec_impl.f90 @@ -195,7 +195,8 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) + ! Application of approximate inverse preconditioner, just some spmm select case(trans_) case('N') @@ -373,11 +374,8 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) ! Application of approximate inverse preconditioner, just some spmm - ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & - ! & trans, work) - select case(trans_) @@ -866,7 +864,7 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 endif - ! Computin the factorization + ! Computing the factorization call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) if(info == psb_success_) then @@ -885,6 +883,151 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_invk_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(n) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invk_fact(a,fill_in, inv_fill,lf,dd,uf,desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_invt_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(eps) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invt_fact(a,fill_in,inv_fill,fact_eps,inv_thresh,lf,dd,uf,& + & desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' diff --git a/prec/impl/psb_s_bjacprec_impl.f90 b/prec/impl/psb_s_bjacprec_impl.f90 index ec2fffd2..3abe5c5d 100644 --- a/prec/impl/psb_s_bjacprec_impl.f90 +++ b/prec/impl/psb_s_bjacprec_impl.f90 @@ -195,7 +195,8 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) + ! Application of approximate inverse preconditioner, just some spmm select case(trans_) case('N') @@ -373,11 +374,8 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) ! Application of approximate inverse preconditioner, just some spmm - ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & - ! & trans, work) - select case(trans_) @@ -866,7 +864,7 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 endif - ! Computin the factorization + ! Computing the factorization call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) if(info == psb_success_) then @@ -885,6 +883,151 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_invk_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(n) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invk_fact(a,fill_in, inv_fill,lf,dd,uf,desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_invt_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(eps) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invt_fact(a,fill_in,inv_fill,fact_eps,inv_thresh,lf,dd,uf,& + & desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' diff --git a/prec/impl/psb_z_bjacprec_impl.f90 b/prec/impl/psb_z_bjacprec_impl.f90 index 4ceb489f..321a9768 100644 --- a/prec/impl/psb_z_bjacprec_impl.f90 +++ b/prec/impl/psb_z_bjacprec_impl.f90 @@ -195,7 +195,8 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) + ! Application of approximate inverse preconditioner, just some spmm select case(trans_) case('N') @@ -373,11 +374,8 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) goto 9999 end if - case(psb_f_ainv_) + case(psb_f_ainv_,psb_f_invt_,psb_f_invk_) ! Application of approximate inverse preconditioner, just some spmm - ! call psb_spmm(alpha, a, x, beta, y,desc_a, info, & - ! & trans, work) - select case(trans_) @@ -866,7 +864,7 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') goto 9999 endif - ! Computin the factorization + ! Computing the factorization call psb_ainv_fact(a,iinvalg,inv_fill,inv_thresh,lf,dd,uf,desc_a,info,iscale=iscale) if(info == psb_success_) then @@ -885,6 +883,151 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end if + case(psb_f_invk_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(n) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invk_fact(a,fill_in, inv_fill,lf,dd,uf,desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_invt_) + ! Approximate Inverse Factorizations based on the sparse inversion of + ! triangular factors of an ILU(eps) factorization + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! Computing the factorization + call psb_invt_fact(a,fill_in,inv_fill,fact_eps,inv_thresh,lf,dd,uf,& + & desc_a,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilut_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' diff --git a/prec/psb_c_bjacprec.f90 b/prec/psb_c_bjacprec.f90 index 847baf00..c70d57e7 100644 --- a/prec/psb_c_bjacprec.f90 +++ b/prec/psb_c_bjacprec.f90 @@ -34,6 +34,8 @@ module psb_c_bjacprec use psb_c_base_prec_mod use psb_c_ilu_fact_mod use psb_c_ainv_fact_mod + use psb_c_invk_fact_mod + use psb_c_invt_fact_mod type, extends(psb_c_base_prec_type) :: psb_c_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -59,9 +61,9 @@ module psb_c_bjacprec end type psb_c_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:5)=(/'None ','ILU(0) ',& + & fact_names(0:6)=(/'None ','ILU(0) ',& & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& - & 'INVT '/) + & 'INVT ','INVK '/) private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros diff --git a/prec/psb_d_bjacprec.f90 b/prec/psb_d_bjacprec.f90 index 1cb90e05..2a80ec2f 100644 --- a/prec/psb_d_bjacprec.f90 +++ b/prec/psb_d_bjacprec.f90 @@ -34,6 +34,8 @@ module psb_d_bjacprec use psb_d_base_prec_mod use psb_d_ilu_fact_mod use psb_d_ainv_fact_mod + use psb_d_invk_fact_mod + use psb_d_invt_fact_mod type, extends(psb_d_base_prec_type) :: psb_d_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -59,9 +61,9 @@ module psb_d_bjacprec end type psb_d_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:5)=(/'None ','ILU(0) ',& + & fact_names(0:6)=(/'None ','ILU(0) ',& & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& - & 'INVT '/) + & 'INVT ','INVK '/) private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros diff --git a/prec/psb_prec_const_mod.f90 b/prec/psb_prec_const_mod.f90 index d8d9a43f..46c0f11b 100644 --- a/prec/psb_prec_const_mod.f90 +++ b/prec/psb_prec_const_mod.f90 @@ -57,8 +57,8 @@ module psb_prec_const_mod integer(psb_ipk_), parameter :: psb_rfpsz=8 ! Factorization types: none, ILU(0), ILU(N), ILU(N,E) integer(psb_ipk_), parameter :: psb_f_none_=0,psb_f_ilu_n_=1,psb_f_ilu_k_=2,psb_f_ilu_t_=3 - ! Approximate Inverse factorization type: AINV - integer(psb_ipk_), parameter :: psb_f_ainv_=4, psb_f_invt_=5 + ! Approximate Inverse factorization type: AINV, INVT, INVK + integer(psb_ipk_), parameter :: psb_f_ainv_=4, psb_f_invt_=5, psb_f_invk_=6 ! Fields for sparse matrices ensembles: integer(psb_ipk_), parameter :: psb_l_pr_=1, psb_u_pr_=2, psb_bp_ilu_avsz=2 integer(psb_ipk_), parameter :: psb_max_avsz=psb_bp_ilu_avsz diff --git a/prec/psb_s_bjacprec.f90 b/prec/psb_s_bjacprec.f90 index 028aabe8..2d42d615 100644 --- a/prec/psb_s_bjacprec.f90 +++ b/prec/psb_s_bjacprec.f90 @@ -34,6 +34,8 @@ module psb_s_bjacprec use psb_s_base_prec_mod use psb_s_ilu_fact_mod use psb_s_ainv_fact_mod + use psb_s_invk_fact_mod + use psb_s_invt_fact_mod type, extends(psb_s_base_prec_type) :: psb_s_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -59,9 +61,9 @@ module psb_s_bjacprec end type psb_s_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:5)=(/'None ','ILU(0) ',& + & fact_names(0:6)=(/'None ','ILU(0) ',& & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& - & 'INVT '/) + & 'INVT ','INVK '/) private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros diff --git a/prec/psb_z_bjacprec.f90 b/prec/psb_z_bjacprec.f90 index 4ca879db..7419359b 100644 --- a/prec/psb_z_bjacprec.f90 +++ b/prec/psb_z_bjacprec.f90 @@ -34,6 +34,8 @@ module psb_z_bjacprec use psb_z_base_prec_mod use psb_z_ilu_fact_mod use psb_z_ainv_fact_mod + use psb_z_invk_fact_mod + use psb_z_invt_fact_mod type, extends(psb_z_base_prec_type) :: psb_z_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) @@ -59,9 +61,9 @@ module psb_z_bjacprec end type psb_z_bjac_prec_type character(len=15), parameter, private :: & - & fact_names(0:5)=(/'None ','ILU(0) ',& + & fact_names(0:6)=(/'None ','ILU(0) ',& & 'ILU(n) ','ILU(eps) ','AINV(eps) ',& - & 'INVT '/) + & 'INVT ','INVK '/) private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros