diff --git a/prec/Makefile b/prec/Makefile index 5fe1b595..eccd6814 100644 --- a/prec/Makefile +++ b/prec/Makefile @@ -2,12 +2,10 @@ include ../Make.inc LIBDIR=../lib HERE=. -MODOBJS= psb_prec_type.o psb_prec_mod.o -F90OBJS= psb_dbjac_bld.o psb_dilu_fct.o\ +MODOBJS= psb_prec_type.o psb_prec_mod.o psb_d_diagprec.o psb_d_nullprec.o \ + psb_d_bjacprec.o +F90OBJS= psb_dilu_fct.o\ psb_dprecbld.o psb_dprecset.o psb_dprecinit.o \ - psb_ddiagsc_bld.o \ - psb_dprc_aply.o \ - psb_dgprec_aply.o psb_dbjac_aply.o\ psb_sbjac_bld.o psb_silu_fct.o\ psb_sprecbld.o psb_sprecset.o psb_sprecinit.o \ psb_sdiagsc_bld.o \ @@ -39,6 +37,7 @@ $(OBJS): $(LIBDIR)/psb_base_mod$(.mod) $(F90OBJS): $(MODOBJS) psb_prec_mod.o: psb_prec_type.o +psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o: psb_prec_type.o psb_prec_mod.o veryclean: clean /bin/rm -f $(LIBNAME) diff --git a/prec/psb_d_bjacprec.f03 b/prec/psb_d_bjacprec.f03 new file mode 100644 index 00000000..724e03e4 --- /dev/null +++ b/prec/psb_d_bjacprec.f03 @@ -0,0 +1,560 @@ +module psb_d_bjacprec + use psb_prec_type + + + type, extends(psb_d_base_prec_type) :: psb_d_bjac_prec_type + integer, allocatable :: iprcparm(:) + type(psb_d_sparse_mat), allocatable :: av(:) + real(psb_dpk_), allocatable :: d(:) + contains + procedure, pass(prec) :: apply => d_bjac_apply + procedure, pass(prec) :: precbld => d_bjac_precbld + procedure, pass(prec) :: precinit => d_bjac_precinit + procedure, pass(prec) :: d_base_precseti => d_bjac_precseti + procedure, pass(prec) :: d_base_precsetr => d_bjac_precsetr + procedure, pass(prec) :: d_base_precsetc => d_bjac_precsetc + procedure, pass(prec) :: precfree => d_bjac_precfree + procedure, pass(prec) :: precdescr => d_bjac_precdescr + procedure, pass(prec) :: sizeof => d_bjac_sizeof + end type psb_d_bjac_prec_type + + + character(len=15), parameter, private :: & + & fact_names(0:2)=(/'None ','ILU(n) ',& + & 'ILU(eps) '/) + +contains + + + subroutine d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(psb_d_bjac_prec_type), intent(in) :: prec + real(psb_dpk_),intent(in) :: alpha,beta + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + + ! Local variables + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:), aux(:) + integer :: ictxt,np,me, err_act, int_err(5) + integer :: debug_level, debug_unit + character :: trans_ + character(len=20) :: name='d_bjac_prec_apply' + character(len=20) :: ch_err + + info = 0 + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = psb_cd_get_context(desc_data) + call psb_info(ictxt, me, np) + + + trans_ = psb_toupper(trans) + select case(trans_) + case('N','T','C') + ! Ok + case default + call psb_errpush(40,name) + goto 9999 + end select + + + n_row = psb_cd_get_local_rows(desc_data) + n_col = psb_cd_get_local_cols(desc_data) + + if (size(x) < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/)) + goto 9999 + end if + if (size(y) < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/)) + goto 9999 + end if + if (.not.allocated(prec%d)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (size(prec%d) < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + endif + + + select case(prec%iprcparm(psb_f_type_)) + case(psb_f_ilu_n_) + + select case(trans_) + case('N') + + call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,& + & trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux) + if(info /=0) goto 9999 + call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + if(info /=0) goto 9999 + + case('T','C') + call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,ww,desc_data,info,& + & trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux) + if(info /=0) goto 9999 + call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + if(info /=0) goto 9999 + + end select + + + case default + info = 4001 + call psb_errpush(info,name,a_err='Invalid factorization') + goto 9999 + end select + + call psb_halo(y,desc_data,info,data=psb_comm_mov_) + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + + end subroutine d_bjac_apply + + subroutine d_bjac_precinit(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_null_precinit' + + call psb_erractionsave(err_act) + + info = 0 + call psb_realloc(psb_ifpsz,prec%iprcparm,info) + if (info /= 0) then + info = 4000 + call psb_Errpush(info,name) + goto 9999 + end if + + prec%iprcparm(:) = 0 + prec%iprcparm(psb_p_type_) = psb_bjac_ + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + prec%iprcparm(psb_ilu_fill_in_) = 0 + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_bjac_precinit + + + subroutine d_bjac_precbld(a,desc_a,prec,info,upd) + + use psb_base_mod + use psb_prec_mod + Implicit None + + type(psb_d_sparse_mat), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer, intent(out) :: info + character, intent(in), optional :: upd + + ! .. Local Scalars .. + integer :: i, m + integer :: int_err(5) + character :: trans, unitd + type(psb_d_csr_sparse_mat), allocatable :: lf, uf + real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8 + integer nztota, err_act, n_row, nrow_a,n_col, nhalo + integer :: ictxt,np,me + character(len=20) :: name='d_bjac_precbld' + character(len=20) :: ch_err + + + if(psb_get_errstatus() /= 0) return + info = 0 + + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt, me, np) + + m = a%get_nrows() + if (m < 0) then + info = 10 + int_err(1) = 1 + int_err(2) = m + call psb_errpush(info,name,i_err=int_err) + goto 9999 + endif + trans = 'N' + unitd = 'U' + + select case(prec%iprcparm(psb_f_type_)) + + case(psb_f_ilu_n_) + + 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 /= 0) then + call psb_errpush(4000,name) + goto 9999 + end if + endif + + nrow_a = psb_cd_get_local_rows(desc_a) + nztota = a%get_nzeros() + + n_col = psb_cd_get_local_cols(desc_a) + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == 0) call lf%allocate(n_row,n_row,nztota) + if (info == 0) call uf%allocate(n_row,n_row,nztota) + + if(info/=0) then + info=4010 + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (allocated(prec%d)) then + if (size(prec%d) < n_row) then + deallocate(prec%d) + endif + endif + if (.not.allocated(prec%d)) then + allocate(prec%d(n_row),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + endif + t3 = psb_wtime() + ! This is where we have no renumbering, thus no need + call psb_ilu_fct(a,lf,uf,prec%d,info) + + if(info==0) then + call prec%av(psb_l_pr_)%mv_from(lf) + call prec%av(psb_u_pr_)%mv_from(uf) + 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() + else + info=4010 + ch_err='psb_ilu_fct' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_none_) + info=4010 + ch_err='Inconsistent prec psb_f_none_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + + case default + info=4010 + ch_err='Unknown psb_f_type_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_bjac_precbld + + subroutine d_bjac_precseti(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer, intent(in) :: what + integer, intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_bjac_precset' + + call psb_erractionsave(err_act) + + info = 0 + if (.not.allocated(prec%iprcparm)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + + select case(what) + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_f_type_) = val + + case (psb_ilu_fill_in_) + if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.(prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_fill_in_) = val + + case default + write(0,*) 'WHAT is invalid, ignoring user specification' + + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_bjac_precseti + + subroutine d_bjac_precsetr(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer, intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_bjac_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_bjac_precsetr + + subroutine d_bjac_precsetc(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer, intent(in) :: what + character(len=*), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_bjac_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_bjac_precsetc + + subroutine d_bjac_precfree(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_bjac_prec_type), intent(inout) :: prec + integer, intent(out) :: info + + Integer :: err_act, i + character(len=20) :: name='d_bjac_precfree' + + call psb_erractionsave(err_act) + + info = 0 + if (allocated(prec%av)) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + end if + if (allocated(prec%d)) then + deallocate(prec%d,stat=info) + end if + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_bjac_precfree + + + subroutine d_bjac_precdescr(prec,iout) + + use psb_base_mod + Implicit None + + class(psb_d_bjac_prec_type), intent(in) :: prec + integer, intent(in), optional :: iout + + Integer :: err_act, nrow, info + character(len=20) :: name='d_bjac_precdescr' + integer :: iout_ + + call psb_erractionsave(err_act) + + info = 0 + + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + + if (.not.allocated(prec%iprcparm)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + + write(iout_,*) 'Block Jacobi with: ',& + & fact_names(prec%iprcparm(psb_f_type_)) + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_bjac_precdescr + + function d_bjac_sizeof(prec) result(val) + use psb_base_mod + class(psb_d_bjac_prec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + + val = 0 + if (allocated(prec%d)) then + val = val + psb_sizeof_dp * size(prec%d) + endif + if (allocated(prec%av)) then + val = val + psb_sizeof(prec%av(psb_l_pr_)) + val = val + psb_sizeof(prec%av(psb_u_pr_)) + endif + return + end function d_bjac_sizeof + +end module psb_d_bjacprec diff --git a/prec/psb_d_diagprec.f03 b/prec/psb_d_diagprec.f03 new file mode 100644 index 00000000..f3d2f651 --- /dev/null +++ b/prec/psb_d_diagprec.f03 @@ -0,0 +1,352 @@ +module psb_d_diagprec + use psb_prec_type + + + type, extends(psb_d_base_prec_type) :: psb_d_diag_prec_type + real(psb_dpk_), allocatable :: d(:) + contains + procedure, pass(prec) :: apply => d_diag_apply + procedure, pass(prec) :: precbld => d_diag_precbld + procedure, pass(prec) :: precinit => d_diag_precinit + procedure, pass(prec) :: d_base_precseti => d_diag_precseti + procedure, pass(prec) :: d_base_precsetr => d_diag_precsetr + procedure, pass(prec) :: d_base_precsetc => d_diag_precsetc + procedure, pass(prec) :: precfree => d_diag_precfree + procedure, pass(prec) :: precdescr => d_diag_precdescr + procedure, pass(prec) :: sizeof => d_diag_sizeof + end type psb_d_diag_prec_type + + +contains + + + subroutine d_diag_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(psb_d_diag_prec_type), intent(in) :: prec + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(in) :: alpha, beta + real(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + Integer :: err_act, nrow + character(len=20) :: name='d_diag_prec_apply' + real(psb_dpk_), pointer :: ww(:) + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the DIAG preonditioner??? + ! + info = 0 + + nrow = psb_cd_get_local_rows(desc_data) + if (size(x) < nrow) then + info = 36 + call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) + goto 9999 + end if + if (size(y) < nrow) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) + goto 9999 + end if + if (.not.allocated(prec%d)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (size(prec%d) < nrow) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + if (size(work) >= size(x)) then + ww => work + else + allocate(ww(size(x)),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/size(x),0,0,0,0/),a_err='real(psb_dpk_)') + goto 9999 + end if + end if + + ww(1:nrow) = x(1:nrow)*prec%d(1:nrow) + call psb_geaxpby(alpha,ww,beta,y,desc_data,info) + + if (size(work) < size(x)) then + deallocate(ww,stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Deallocate') + goto 9999 + end if + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_diag_apply + + subroutine d_diag_precinit(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_diag_prec_type),intent(inout) :: prec + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_diag_precinit' + + call psb_erractionsave(err_act) + + info = 0 + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_diag_precinit + + + subroutine d_diag_precbld(a,desc_a,prec,info,upd) + + use psb_base_mod + Implicit None + + type(psb_d_sparse_mat), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + class(psb_d_diag_prec_type),intent(inout) :: prec + integer, intent(out) :: info + character, intent(in), optional :: upd + Integer :: err_act, nrow,i + character(len=20) :: name='d_diag_precbld' + + call psb_erractionsave(err_act) + + info = 0 + nrow = psb_cd_get_local_cols(desc_a) + if (allocated(prec%d)) then + if (size(prec%d) < nrow) then + deallocate(prec%d,stat=info) + end if + end if + if ((info == 0).and.(.not.allocated(prec%d))) then + allocate(prec%d(nrow), stat=info) + end if + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + + call a%get_diag(prec%d,info) + if (info /= 0) then + info = 4010 + call psb_errpush(info,name, a_err='get_diag') + goto 9999 + end if + + do i=1,nrow + if (prec%d(i) == dzero) then + prec%d(i) = done + else + prec%d(i) = done/prec%d(i) + endif + end do + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_diag_precbld + + subroutine d_diag_precseti(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_diag_prec_type),intent(inout) :: prec + integer, intent(in) :: what + integer, intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_diag_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_diag_precseti + + subroutine d_diag_precsetr(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_diag_prec_type),intent(inout) :: prec + integer, intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_diag_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_diag_precsetr + + subroutine d_diag_precsetc(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_diag_prec_type),intent(inout) :: prec + integer, intent(in) :: what + character(len=*), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_diag_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_diag_precsetc + + subroutine d_diag_precfree(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_diag_prec_type), intent(inout) :: prec + integer, intent(out) :: info + + Integer :: err_act, nrow + character(len=20) :: name='d_diag_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_diag_precfree + + + subroutine d_diag_precdescr(prec,iout) + + use psb_base_mod + Implicit None + + class(psb_d_diag_prec_type), intent(in) :: prec + integer, intent(in), optional :: iout + + Integer :: err_act, nrow, info + character(len=20) :: name='d_diag_precdescr' + + integer :: iout_ + + call psb_erractionsave(err_act) + + info = 0 + + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + + write(iout_,*) 'Diagonal scaling' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_diag_precdescr + + function d_diag_sizeof(prec) result(val) + use psb_base_mod + class(psb_d_diag_prec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + + val = 0 + val = val + psb_sizeof_dp * size(prec%d) + return + end function d_diag_sizeof + +end module psb_d_diagprec diff --git a/prec/psb_d_nullprec.f03 b/prec/psb_d_nullprec.f03 new file mode 100644 index 00000000..e648be62 --- /dev/null +++ b/prec/psb_d_nullprec.f03 @@ -0,0 +1,293 @@ +module psb_d_nullprec + use psb_prec_type + + + type, extends(psb_d_base_prec_type) :: psb_d_null_prec_type + contains + procedure, pass(prec) :: apply => d_null_apply + procedure, pass(prec) :: precbld => d_null_precbld + procedure, pass(prec) :: precinit => d_null_precinit + procedure, pass(prec) :: d_base_precseti => d_null_precseti + procedure, pass(prec) :: d_base_precsetr => d_null_precsetr + procedure, pass(prec) :: d_base_precsetc => d_null_precsetc + procedure, pass(prec) :: precfree => d_null_precfree + procedure, pass(prec) :: precdescr => d_null_precdescr + procedure, pass(prec) :: sizeof => d_null_sizeof + end type psb_d_null_prec_type + + +contains + + + subroutine d_null_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(psb_d_null_prec_type), intent(in) :: prec + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(in) :: alpha, beta + real(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + Integer :: err_act, nrow + character(len=20) :: name='d_null_prec_apply' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 0 + + nrow = psb_cd_get_local_rows(desc_data) + if (size(x) < nrow) then + info = 36 + call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) + goto 9999 + end if + if (size(y) < nrow) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) + goto 9999 + end if + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + if (info /= 0 ) then + info = 4010 + call psb_errpush(infoi,name,a_err="psb_geaxpby") + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_null_apply + + + subroutine d_null_precinit(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_null_prec_type),intent(inout) :: prec + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_null_precinit' + + call psb_erractionsave(err_act) + + info = 0 + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_null_precinit + + subroutine d_null_precbld(a,desc_a,prec,info,upd) + + use psb_base_mod + Implicit None + + type(psb_d_sparse_mat), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + class(psb_d_null_prec_type),intent(inout) :: prec + integer, intent(out) :: info + character, intent(in), optional :: upd + Integer :: err_act, nrow + character(len=20) :: name='d_null_precbld' + + call psb_erractionsave(err_act) + + info = 0 + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_null_precbld + + subroutine d_null_precseti(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_null_prec_type),intent(inout) :: prec + integer, intent(in) :: what + integer, intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_null_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_null_precseti + + subroutine d_null_precsetr(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_null_prec_type),intent(inout) :: prec + integer, intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_null_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_null_precsetr + + subroutine d_null_precsetc(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_null_prec_type),intent(inout) :: prec + integer, intent(in) :: what + character(len=*), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_null_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_null_precsetc + + subroutine d_null_precfree(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_null_prec_type), intent(inout) :: prec + integer, intent(out) :: info + + Integer :: err_act, nrow + character(len=20) :: name='d_null_precset' + + call psb_erractionsave(err_act) + + info = 0 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_null_precfree + + + subroutine d_null_precdescr(prec,iout) + + use psb_base_mod + Implicit None + + class(psb_d_null_prec_type), intent(in) :: prec + integer, intent(in), optional :: iout + + Integer :: err_act, nrow, info + character(len=20) :: name='d_null_precset' + integer :: iout_ + + call psb_erractionsave(err_act) + + info = 0 + + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + end if + + write(iout_,*) 'No preconditioning' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_null_precdescr + + function d_null_sizeof(prec) result(val) + use psb_base_mod + class(psb_d_null_prec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + + val = 0 + + return + end function d_null_sizeof + +end module psb_d_nullprec diff --git a/prec/psb_dprc_aply.f90 b/prec/psb_dprc_aply.f90 index d2e4692b..cb3b1124 100644 --- a/prec/psb_dprc_aply.f90 +++ b/prec/psb_dprc_aply.f90 @@ -73,8 +73,13 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) end if - call psb_gprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info) - + if (.not.allocated(prec%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if +!!$ call psb_gprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info) + call prec%dprec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_) if (present(work)) then else deallocate(work_) @@ -139,7 +144,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) type(psb_desc_type),intent(in) :: desc_data type(psb_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: x(:) integer, intent(out) :: info character(len=1), optional :: trans @@ -161,12 +166,17 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) trans_='N' end if + if (.not.allocated(prec%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if allocate(ww(size(x)),w1(size(x)),stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - call psb_dprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) + call prec%dprec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1) if(info /=0) goto 9999 x(:) = ww(:) deallocate(ww,W1) diff --git a/prec/psb_dprecbld.f90 b/prec/psb_dprecbld.f90 index af7690fc..7e5bfe23 100644 --- a/prec/psb_dprecbld.f90 +++ b/prec/psb_dprecbld.f90 @@ -81,63 +81,74 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) ! ALso should define symbolic names for the preconditioners. ! - call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',& - & psb_diag_,is_legal_prec) - - call psb_nullify_desc(p%desc_data) - - select case(p%iprcparm(psb_p_type_)) - case (psb_noprec_) - ! Do nothing. - call psb_cdcpy(desc_a,p%desc_data,info) - if(info /= 0) then - info=4010 - ch_err='psb_cdcpy' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case (psb_diag_) - - call psb_diagsc_bld(a,desc_a,p,upd_,info) - if(info /= 0) then - info=4010 - ch_err='psb_diagsc_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case (psb_bjac_) - - call psb_check_def(p%iprcparm(psb_f_type_),'fact',& - & psb_f_ilu_n_,is_legal_ml_fact) - - call psb_bjac_bld(a,desc_a,p,upd_,info) - - if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_bjac_bld') + if (.false.) then +!!$ call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',& +!!$ & psb_diag_,is_legal_prec) +!!$ +!!$ call psb_nullify_desc(p%desc_data) +!!$ +!!$ select case(p%iprcparm(psb_p_type_)) +!!$ case (psb_noprec_) +!!$ ! Do nothing. +!!$ call psb_cdcpy(desc_a,p%desc_data,info) +!!$ if(info /= 0) then +!!$ info=4010 +!!$ ch_err='psb_cdcpy' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ +!!$ case (psb_diag_) +!!$ +!!$ call psb_diagsc_bld(a,desc_a,p,upd_,info) +!!$ if(info /= 0) then +!!$ info=4010 +!!$ ch_err='psb_diagsc_bld' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ end if +!!$ +!!$ case (psb_bjac_) +!!$ +!!$ call psb_check_def(p%iprcparm(psb_f_type_),'fact',& +!!$ & psb_f_ilu_n_,is_legal_ml_fact) +!!$ +!!$ call psb_bjac_bld(a,desc_a,p,upd_,info) +!!$ +!!$ if(info /= 0) then +!!$ call psb_errpush(4010,name,a_err='psb_bjac_bld') +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ info=4010 +!!$ ch_err='Unknown psb_p_type_' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ +!!$ end select + else + if (.not.allocated(p%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if + + call p%dprec%precbld(a,desc_a,info,upd) + if (info /= 0) goto 9999 - case default - info=4010 - ch_err='Unknown psb_p_type_' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - - end select - + end if call psb_erractionrestore(err_act) return 9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if return - end if - return -end subroutine psb_dprecbld + end subroutine psb_dprecbld diff --git a/prec/psb_dprecinit.f90 b/prec/psb_dprecinit.f90 index 2ccfc58e..381fec46 100644 --- a/prec/psb_dprecinit.f90 +++ b/prec/psb_dprecinit.f90 @@ -33,6 +33,9 @@ subroutine psb_dprecinit(p,ptype,info) use psb_base_mod use psb_prec_mod, psb_protect_name => psb_dprecinit + use psb_d_nullprec + use psb_d_diagprec + use psb_d_bjacprec implicit none type(psb_dprec_type), intent(inout) :: p character(len=*), intent(in) :: ptype @@ -40,31 +43,57 @@ subroutine psb_dprecinit(p,ptype,info) info = 0 - call psb_realloc(psb_ifpsz,p%iprcparm,info) - if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info) - if (info /= 0) return + if (.false.) then +!!$ call psb_realloc(psb_ifpsz,p%iprcparm,info) +!!$ if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info) +!!$ if (info /= 0) return +!!$ +!!$ select case(psb_toupper(ptype(1:len_trim(ptype)))) +!!$ case ('NONE','NOPREC') +!!$ p%iprcparm(:) = 0 +!!$ p%iprcparm(psb_p_type_) = psb_noprec_ +!!$ p%iprcparm(psb_f_type_) = psb_f_none_ +!!$ +!!$ case ('DIAG') +!!$ p%iprcparm(:) = 0 +!!$ p%iprcparm(psb_p_type_) = psb_diag_ +!!$ p%iprcparm(psb_f_type_) = psb_f_none_ +!!$ +!!$ case ('BJAC') +!!$ p%iprcparm(:) = 0 +!!$ p%iprcparm(psb_p_type_) = psb_bjac_ +!!$ p%iprcparm(psb_f_type_) = psb_f_ilu_n_ +!!$ p%iprcparm(psb_ilu_fill_in_) = 0 +!!$ +!!$ case default +!!$ write(0,*) 'Unknown preconditioner type request "',ptype,'"' +!!$ info = 2 +!!$ +!!$ end select + else + if (allocated(p%dprec) ) then + call p%dprec%precfree(info) + if (info == 0) deallocate(p%dprec,stat=info) + if (info /= 0) return + end if - select case(psb_toupper(ptype(1:len_trim(ptype)))) - case ('NONE','NOPREC') - p%iprcparm(:) = 0 - p%iprcparm(psb_p_type_) = psb_noprec_ - p%iprcparm(psb_f_type_) = psb_f_none_ + select case(psb_toupper(ptype(1:len_trim(ptype)))) + case ('NONE','NOPREC') + + allocate(psb_d_null_prec_type :: p%dprec, stat=info) + + case ('DIAG') + allocate(psb_d_diag_prec_type :: p%dprec, stat=info) - case ('DIAG') - p%iprcparm(:) = 0 - p%iprcparm(psb_p_type_) = psb_diag_ - p%iprcparm(psb_f_type_) = psb_f_none_ + case ('BJAC') + allocate(psb_d_bjac_prec_type :: p%dprec, stat=info) - case ('BJAC') - p%iprcparm(:) = 0 - p%iprcparm(psb_p_type_) = psb_bjac_ - p%iprcparm(psb_f_type_) = psb_f_ilu_n_ - p%iprcparm(psb_ilu_fill_in_) = 0 + case default + write(0,*) 'Unknown preconditioner type request "',ptype,'"' + info = 2 - case default - write(0,*) 'Unknown preconditioner type request "',ptype,'"' - info = 2 - - end select + end select + if (info == 0) call p%dprec%precinit(info) + end if end subroutine psb_dprecinit diff --git a/prec/psb_dprecset.f90 b/prec/psb_dprecset.f90 index 06070d64..54989e56 100644 --- a/prec/psb_dprecset.f90 +++ b/prec/psb_dprecset.f90 @@ -37,30 +37,38 @@ subroutine psb_dprecseti(p,what,val,info) type(psb_dprec_type), intent(inout) :: p integer :: what, val integer, intent(out) :: info + character(len=20) :: name='precset' info = 0 - return - select case(what) - case (psb_f_type_) - if (p%iprcparm(psb_p_type_) /= psb_bjac_) then - write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif - p%iprcparm(psb_f_type_) = val - - case (psb_ilu_fill_in_) - if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then - write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),& - & 'ignoring user specification' - return - endif - p%iprcparm(psb_ilu_fill_in_) = val - - case default - write(0,*) 'WHAT is invalid, ignoring user specification' + if (.not.allocated(p%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + return +!!$ goto 9999 + end if - end select + call p%dprec%precset(what,val,info) +!!$ select case(what) +!!$ case (psb_f_type_) +!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then +!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),& +!!$ & 'ignoring user specification' +!!$ return +!!$ endif +!!$ p%iprcparm(psb_f_type_) = val +!!$ +!!$ case (psb_ilu_fill_in_) +!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then +!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),& +!!$ & 'ignoring user specification' +!!$ return +!!$ endif +!!$ p%iprcparm(psb_ilu_fill_in_) = val +!!$ +!!$ case default +!!$ write(0,*) 'WHAT is invalid, ignoring user specification' +!!$ +!!$ end select return end subroutine psb_dprecseti @@ -75,12 +83,22 @@ subroutine psb_dprecsetd(p,what,val,info) integer :: what real(psb_dpk_) :: val integer, intent(out) :: info + character(len=20) :: name='precset' + + info = 0 + if (.not.allocated(p%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + return +!!$ goto 9999 + end if -! -! This will have to be changed if/when we put together an ILU(eps) -! factorization. -! - select case(what) + call p%dprec%precset(what,val,info) +!!$! +!!$! This will have to be changed if/when we put together an ILU(eps) +!!$! factorization. +!!$! +!!$ select case(what) !!$ case (psb_f_type_) !!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then !!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),& @@ -96,11 +114,11 @@ subroutine psb_dprecsetd(p,what,val,info) !!$ return !!$ endif !!$ p%iprcparm(psb_ilu_fill_in_) = val - - case default - write(0,*) 'WHAT is invalid, ignoring user specification' - - end select +!!$ +!!$ case default +!!$ write(0,*) 'WHAT is invalid, ignoring user specification' +!!$ +!!$ end select return end subroutine psb_dprecsetd diff --git a/prec/psb_prec_type.f03 b/prec/psb_prec_type.f03 index 0e8da866..a12c1ac0 100644 --- a/prec/psb_prec_type.f03 +++ b/prec/psb_prec_type.f03 @@ -78,14 +78,23 @@ module psb_prec_type generic, public :: apply => s_apply2v, s_apply1v end type psb_sprec_type + + type psb_d_base_prec_type + contains + procedure, pass(prec) :: apply => d_base_apply + procedure, pass(prec) :: precbld => d_base_precbld + procedure, pass(prec) :: d_base_precseti + procedure, pass(prec) :: d_base_precsetr + procedure, pass(prec) :: d_base_precsetc + procedure, pass(prec) :: sizeof => d_base_sizeof + generic, public :: precset => d_base_precseti, d_base_precsetr, d_base_precsetc + procedure, pass(prec) :: precinit => d_base_precinit + procedure, pass(prec) :: precfree => d_base_precfree + procedure, pass(prec) :: precdescr => d_base_precdescr + end type psb_d_base_prec_type + type psb_dprec_type - type(psb_d_sparse_mat), allocatable :: av(:) - real(psb_dpk_), allocatable :: d(:) - type(psb_desc_type) :: desc_data - integer, allocatable :: iprcparm(:) - real(psb_dpk_), allocatable :: rprcparm(:) - integer, allocatable :: perm(:), invperm(:) - integer :: prec + class(psb_d_base_prec_type), allocatable :: dprec contains procedure, pass(prec) :: d_apply2v procedure, pass(prec) :: d_apply1v @@ -148,8 +157,10 @@ module psb_prec_type end interface interface psb_sizeof - module procedure psb_sprec_sizeof, psb_dprec_sizeof,& - & psb_cprec_sizeof, psb_zprec_sizeof + module procedure psb_sprec_sizeof, & + & psb_dprec_sizeof,& + & psb_cprec_sizeof, & + & psb_zprec_sizeof end interface @@ -161,8 +172,8 @@ module psb_prec_type import psb_sprec_type type(psb_desc_type),intent(in) :: desc_data type(psb_sprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: x(:) + real(psb_spk_),intent(inout) :: y(:) integer, intent(out) :: info character(len=1), optional :: trans real(psb_spk_),intent(inout), optional, target :: work(:) @@ -192,7 +203,7 @@ module psb_prec_type import psb_dprec_type type(psb_desc_type),intent(in) :: desc_data type(psb_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(inout) :: x(:) + real(psb_dpk_),intent(inout) :: x(:) integer, intent(out) :: info character(len=1), optional :: trans end subroutine psb_dprc_aply1 @@ -201,8 +212,8 @@ module psb_prec_type import psb_cprec_type type(psb_desc_type),intent(in) :: desc_data type(psb_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) integer, intent(out) :: info character(len=1), optional :: trans complex(psb_spk_),intent(inout), optional, target :: work(:) @@ -212,7 +223,7 @@ module psb_prec_type import psb_cprec_type type(psb_desc_type),intent(in) :: desc_data type(psb_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: x(:) integer, intent(out) :: info character(len=1), optional :: trans end subroutine psb_cprc_aply1 @@ -221,8 +232,8 @@ module psb_prec_type import psb_zprec_type type(psb_desc_type),intent(in) :: desc_data type(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) integer, intent(out) :: info character(len=1), optional :: trans complex(psb_dpk_),intent(inout), optional, target :: work(:) @@ -232,7 +243,7 @@ module psb_prec_type import psb_zprec_type type(psb_desc_type),intent(in) :: desc_data type(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: x(:) integer, intent(out) :: info character(len=1), optional :: trans end subroutine psb_zprc_aply1 @@ -245,9 +256,11 @@ contains subroutine psb_file_prec_descr(p,iout) + use psb_base_mod type(psb_dprec_type), intent(in) :: p integer, intent(in), optional :: iout integer :: iout_ + character(len=20) :: name='prec_descr' if (present(iout)) then iout_ = iout @@ -255,17 +268,12 @@ contains iout_ = 6 end if - write(iout_,*) 'Preconditioner description' - select case(p%iprcparm(psb_p_type_)) - case(psb_noprec_) - write(iout_,*) 'No preconditioning' - case(psb_diag_) - write(iout_,*) 'Diagonal scaling' - case(psb_bjac_) - write(iout_,*) 'Block Jacobi with: ',& - & fact_names(p%iprcparm(psb_f_type_)) - end select - + if (.not.allocated(p%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + end if + call p%dprec%precdescr(iout) + end subroutine psb_file_prec_descr subroutine psb_sfile_prec_descr(p,iout) @@ -496,40 +504,13 @@ contains me=-1 - ! Actually we migh just deallocate the top level array, except - ! for the inner UMFPACK or SLU stuff - - if (allocated(p%d)) then - deallocate(p%d,stat=info) - end if - - if (allocated(p%av)) then - do i=1,size(p%av) - call p%av(i)%free() - enddo - deallocate(p%av,stat=info) - end if - - if (allocated(p%desc_data%matrix_data)) & - & call psb_cdfree(p%desc_data,info) - - if (allocated(p%rprcparm)) then - deallocate(p%rprcparm,stat=info) - end if - - if (allocated(p%perm)) then - deallocate(p%perm,stat=info) - endif - - if (allocated(p%invperm)) then - deallocate(p%invperm,stat=info) - endif - - if (allocated(p%iprcparm)) then - deallocate(p%iprcparm,stat=info) + + if (allocated(p%dprec)) then + call p%dprec%precfree(info) + if (info /= 0) goto 9999 + deallocate(p%dprec,stat=info) + if (info /= 0) goto 9999 end if - call psb_nullify_prec(p) - call psb_erractionrestore(err_act) return @@ -546,8 +527,6 @@ contains subroutine psb_nullify_dprec(p) type(psb_dprec_type), intent(inout) :: p -!!$ nullify(p%av,p%d,p%iprcparm,p%rprcparm,p%perm,p%invperm,p%mlia,& -!!$ & p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac) end subroutine psb_nullify_dprec @@ -695,18 +674,10 @@ contains integer(psb_long_int_k_) :: val integer :: i val = 0 - if (allocated(prec%iprcparm)) val = val + psb_sizeof_int * size(prec%iprcparm) - if (allocated(prec%rprcparm)) val = val + psb_sizeof_dp * size(prec%rprcparm) - if (allocated(prec%d)) val = val + psb_sizeof_dp * size(prec%d) - if (allocated(prec%perm)) val = val + psb_sizeof_int * size(prec%perm) - if (allocated(prec%invperm)) val = val + psb_sizeof_int * size(prec%invperm) - val = val + psb_sizeof(prec%desc_data) - if (allocated(prec%av)) then - do i=1,size(prec%av) - val = val + psb_sizeof(prec%av(i)) - end do - end if + if (allocated(prec%dprec)) then + val = val + prec%dprec%sizeof() + end if end function psb_dprec_sizeof function psb_sprec_sizeof(prec) result(val) @@ -851,23 +822,56 @@ contains integer, intent(out) :: info character(len=1), optional :: trans real(psb_dpk_),intent(inout), optional, target :: work(:) - Integer :: err_act - character(len=20) :: name='d_prec_apply' - + + character :: trans_ + real(psb_dpk_), pointer :: work_(:) + integer :: ictxt,np,me,err_act + character(len=20) :: name + + name='psb_d_prc_apply' + info = 0 call psb_erractionsave(err_act) - select type(prec) - type is (psb_dprec_type) - call psb_precaply(prec,x,y,desc_data,info,trans,work) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select + ictxt = psb_cd_get_context(desc_data) + call psb_info(ictxt, me, np) + + if (present(trans)) then + trans_=trans + else + trans_='N' + end if + + if (present(work)) then + work_ => work + else + allocate(work_(4*psb_cd_get_local_cols(desc_data)),stat=info) + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + + end if + + if (.not.allocated(prec%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + call prec%dprec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_) + if (present(work)) then + else + deallocate(work_,stat=info) + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err='DeAllocate') + goto 9999 + end if + end if call psb_erractionrestore(err_act) return - + 9999 continue call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then @@ -885,24 +889,51 @@ contains real(psb_dpk_),intent(inout) :: x(:) integer, intent(out) :: info character(len=1), optional :: trans - Integer :: err_act - character(len=20) :: name='d_prec_apply' + character :: trans_ + integer :: ictxt,np,me, err_act + real(psb_dpk_), pointer :: WW(:), w1(:) + character(len=20) :: name + name='psb_d_apply1' + info = 0 call psb_erractionsave(err_act) - select type(prec) - type is (psb_dprec_type) - call psb_precaply(prec,x,desc_data,info,trans) - class default - info = 700 - call psb_errpush(info,name) - goto 9999 - end select + + ictxt=psb_cd_get_context(desc_data) + call psb_info(ictxt, me, np) + if (present(trans)) then + trans_=psb_toupper(trans) + else + trans_='N' + end if + + if (.not.allocated(prec%dprec)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + allocate(ww(size(x)),w1(size(x)),stat=info) + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err='Allocate') + goto 9999 + end if + call prec%dprec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1) + if(info /=0) goto 9999 + x(:) = ww(:) + deallocate(ww,W1,stat=info) + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err='DeAllocate') + goto 9999 + end if + call psb_erractionrestore(err_act) return 9999 continue + call psb_errpush(info,name) call psb_erractionrestore(err_act) if (err_act == psb_act_abort_) then call psb_error() @@ -1052,4 +1083,290 @@ contains end subroutine z_apply1v + + subroutine d_base_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(psb_d_base_prec_type), intent(in) :: prec + real(psb_dpk_),intent(in) :: alpha, beta + real(psb_dpk_),intent(in) :: x(:) + real(psb_dpk_),intent(inout) :: y(:) + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + Integer :: err_act, nrow + character(len=20) :: name='d_base_prec_apply' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_base_apply + + subroutine d_base_precinit(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_base_prec_type),intent(inout) :: prec + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_base_precinit' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_base_precinit + + subroutine d_base_precbld(a,desc_a,prec,info,upd) + + use psb_base_mod + Implicit None + + type(psb_d_sparse_mat), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + class(psb_d_base_prec_type),intent(inout) :: prec + integer, intent(out) :: info + character, intent(in), optional :: upd + Integer :: err_act, nrow + character(len=20) :: name='d_base_precbld' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_base_precbld + + subroutine d_base_precseti(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_base_prec_type),intent(inout) :: prec + integer, intent(in) :: what + integer, intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_base_precseti' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_base_precseti + + subroutine d_base_precsetr(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_base_prec_type),intent(inout) :: prec + integer, intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_base_precsetr' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_base_precsetr + + subroutine d_base_precsetc(prec,what,val,info) + + use psb_base_mod + Implicit None + + class(psb_d_base_prec_type),intent(inout) :: prec + integer, intent(in) :: what + character(len=*), intent(in) :: val + integer, intent(out) :: info + Integer :: err_act, nrow + character(len=20) :: name='d_base_precsetc' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine d_base_precsetc + + subroutine d_base_precfree(prec,info) + + use psb_base_mod + Implicit None + + class(psb_d_base_prec_type), intent(inout) :: prec + integer, intent(out) :: info + + Integer :: err_act, nrow + character(len=20) :: name='d_base_precfree' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_base_precfree + + + subroutine d_base_precdescr(prec,iout) + + use psb_base_mod + Implicit None + + class(psb_d_base_prec_type), intent(in) :: prec + integer, intent(in), optional :: iout + + Integer :: err_act, nrow, info + character(len=20) :: name='d_base_precdescr' + + call psb_erractionsave(err_act) + + ! + ! This is the base version and we should throw an error. + ! Or should it be the NULL preonditioner??? + ! + info = 700 + call psb_errpush(info,name) + goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_base_precdescr + + + function d_base_sizeof(prec) result(val) + use psb_base_mod + class(psb_d_base_prec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + + val = 0 + return + end function d_base_sizeof + + end module psb_prec_type