diff --git a/mlprec/mld_d_as_smoother.f90 b/mlprec/mld_d_as_smoother.f90 index 11878173..7d905ac0 100644 --- a/mlprec/mld_d_as_smoother.f90 +++ b/mlprec/mld_d_as_smoother.f90 @@ -55,16 +55,17 @@ module mld_d_as_smoother type(psb_desc_type) :: desc_data integer :: novr, restr, prol, nd_nnz_tot contains - procedure, pass(sm) :: check => d_as_smoother_check - procedure, pass(sm) :: dump => d_as_smoother_dmp - procedure, pass(sm) :: build => d_as_smoother_bld - procedure, pass(sm) :: apply => d_as_smoother_apply - procedure, pass(sm) :: free => d_as_smoother_free - procedure, pass(sm) :: seti => d_as_smoother_seti - procedure, pass(sm) :: setc => d_as_smoother_setc - procedure, pass(sm) :: setr => d_as_smoother_setr - procedure, pass(sm) :: descr => d_as_smoother_descr - procedure, pass(sm) :: sizeof => d_as_smoother_sizeof + procedure, pass(sm) :: check => d_as_smoother_check + procedure, pass(sm) :: dump => d_as_smoother_dmp + procedure, pass(sm) :: build => d_as_smoother_bld + procedure, pass(sm) :: apply_v => d_as_smoother_apply_vect + procedure, pass(sm) :: apply_a => d_as_smoother_apply + procedure, pass(sm) :: free => d_as_smoother_free + procedure, pass(sm) :: seti => d_as_smoother_seti + procedure, pass(sm) :: setc => d_as_smoother_setc + procedure, pass(sm) :: setr => d_as_smoother_setr + procedure, pass(sm) :: descr => d_as_smoother_descr + procedure, pass(sm) :: sizeof => d_as_smoother_sizeof procedure, pass(sm) :: default => d_as_smoother_default end type mld_d_as_smoother_type @@ -74,7 +75,7 @@ module mld_d_as_smoother & d_as_smoother_setc, d_as_smoother_setr,& & d_as_smoother_descr, d_as_smoother_sizeof, & & d_as_smoother_check, d_as_smoother_default,& - & d_as_smoother_dmp + & d_as_smoother_dmp, d_as_smoother_apply_vect character(len=6), parameter, private :: & & restrict_names(0:4)=(/'none ','halo ',' ',' ',' '/) @@ -151,6 +152,442 @@ contains return end subroutine d_as_smoother_check + subroutine d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_as_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col, nrow_d, i + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me, err_act,isz,int_err(5) + character :: trans_ + character(len=20) :: name='d_as_smoother_apply', ch_err + + call psb_erractionsave(err_act) + + info = psb_success_ + ictxt = desc_data%get_context() + call psb_info(ictxt,me,np) + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + n_row = sm%desc_data%get_local_rows() + n_col = sm%desc_data%get_local_cols() + nrow_d = desc_data%get_local_rows() + isz=max(n_row,N_COL) + if ((6*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + aux => work(3*isz+1:) + else if ((4*isz) <= size(work)) then + aux => work(1:) + allocate(ww(isz),tx(isz),ty(isz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_request_,name,i_err=(/3*isz,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + else if ((3*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + allocate(aux(4*isz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_request_,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + else + allocate(ww(isz),tx(isz),ty(isz),& + &aux(4*isz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_request_,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + + endif + + if ((sm%novr == 0).and.(sweeps == 1)) then + ! + ! Shortcut: in this case it's just the same + ! as Block Jacobi. + ! + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + else + +!!$ +!!$ tx(1:nrow_d) = x(1:nrow_d) +!!$ tx(nrow_d+1:isz) = dzero +!!$ +!!$ +!!$ if (sweeps == 1) then +!!$ +!!$ select case(trans_) +!!$ case('N') +!!$ ! +!!$ ! Get the overlap entries of tx (tx == x) +!!$ ! +!!$ if (sm%restr == psb_halo_) then +!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_halo' +!!$ goto 9999 +!!$ end if +!!$ else if (sm%restr /= psb_none_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ case('T','C') +!!$ ! +!!$ ! With transpose, we have to do it here +!!$ ! +!!$ +!!$ select case (sm%prol) +!!$ +!!$ case(psb_none_) +!!$ ! +!!$ ! Do nothing +!!$ +!!$ case(psb_sum_) +!!$ ! +!!$ ! The transpose of sum is halo +!!$ ! +!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_halo' +!!$ goto 9999 +!!$ end if +!!$ +!!$ case(psb_avg_) +!!$ ! +!!$ ! Tricky one: first we have to scale the overlap entries, +!!$ ! which we can do by assignind mode=0, i.e. no communication +!!$ ! (hence only scaling), then we do the halo +!!$ ! +!!$ call psb_ovrl(tx,sm%desc_data,info,& +!!$ & update=psb_avg_,work=aux,mode=0) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_ovrl' +!!$ goto 9999 +!!$ end if +!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_halo' +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_') +!!$ goto 9999 +!!$ end select +!!$ +!!$ +!!$ case default +!!$ info=psb_err_iarg_invalid_i_ +!!$ int_err(1)=6 +!!$ ch_err(2:2)=trans +!!$ goto 9999 +!!$ end select +!!$ +!!$ call sm%sv%apply(done,tx,dzero,ty,sm%desc_data,trans_,aux,info) +!!$ +!!$ if (info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,& +!!$ & a_err='Error in sub_aply Jacobi Sweeps = 1') +!!$ goto 9999 +!!$ endif +!!$ +!!$ select case(trans_) +!!$ case('N') +!!$ +!!$ select case (sm%prol) +!!$ +!!$ case(psb_none_) +!!$ ! +!!$ ! Would work anyway, but since it is supposed to do nothing ... +!!$ ! call psb_ovrl(ty,sm%desc_data,info,& +!!$ ! & update=sm%prol,work=aux) +!!$ +!!$ +!!$ case(psb_sum_,psb_avg_) +!!$ ! +!!$ ! Update the overlap of ty +!!$ ! +!!$ call psb_ovrl(ty,sm%desc_data,info,& +!!$ & update=sm%prol,work=aux) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_ovrl' +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_') +!!$ goto 9999 +!!$ end select +!!$ +!!$ case('T','C') +!!$ ! +!!$ ! With transpose, we have to do it here +!!$ ! +!!$ if (sm%restr == psb_halo_) then +!!$ call psb_ovrl(ty,sm%desc_data,info,& +!!$ & update=psb_sum_,work=aux) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_ovrl' +!!$ goto 9999 +!!$ end if +!!$ else if (sm%restr /= psb_none_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_') +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ info=psb_err_iarg_invalid_i_ +!!$ int_err(1)=6 +!!$ ch_err(2:2)=trans +!!$ goto 9999 +!!$ end select +!!$ +!!$ +!!$ +!!$ else if (sweeps > 1) then +!!$ +!!$ ! +!!$ ! +!!$ ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver +!!$ ! to compute an approximate solution of a linear system. +!!$ ! +!!$ ! +!!$ ty = dzero +!!$ do i=1, sweeps +!!$ select case(trans_) +!!$ case('N') +!!$ ! +!!$ ! Get the overlap entries of tx (tx == x) +!!$ ! +!!$ if (sm%restr == psb_halo_) then +!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_halo' +!!$ goto 9999 +!!$ end if +!!$ else if (sm%restr /= psb_none_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ case('T','C') +!!$ ! +!!$ ! With transpose, we have to do it here +!!$ ! +!!$ +!!$ select case (sm%prol) +!!$ +!!$ case(psb_none_) +!!$ ! +!!$ ! Do nothing +!!$ +!!$ case(psb_sum_) +!!$ ! +!!$ ! The transpose of sum is halo +!!$ ! +!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_halo' +!!$ goto 9999 +!!$ end if +!!$ +!!$ case(psb_avg_) +!!$ ! +!!$ ! Tricky one: first we have to scale the overlap entries, +!!$ ! which we can do by assignind mode=0, i.e. no communication +!!$ ! (hence only scaling), then we do the halo +!!$ ! +!!$ call psb_ovrl(tx,sm%desc_data,info,& +!!$ & update=psb_avg_,work=aux,mode=0) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_ovrl' +!!$ goto 9999 +!!$ end if +!!$ call psb_halo(tx,sm%desc_data,info,work=aux,data=psb_comm_ext_) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_halo' +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_') +!!$ goto 9999 +!!$ end select +!!$ +!!$ +!!$ case default +!!$ info=psb_err_iarg_invalid_i_ +!!$ int_err(1)=6 +!!$ ch_err(2:2)=trans +!!$ goto 9999 +!!$ end select +!!$ ! +!!$ ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the +!!$ ! block diagonal part and the remaining part of the local matrix +!!$ ! and Y(j) is the approximate solution at sweep j. +!!$ ! +!!$ ww(1:n_row) = tx(1:n_row) +!!$ call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,work=aux,trans=trans_) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info) +!!$ +!!$ if (info /= psb_success_) exit +!!$ +!!$ +!!$ select case(trans_) +!!$ case('N') +!!$ +!!$ select case (sm%prol) +!!$ +!!$ case(psb_none_) +!!$ ! +!!$ ! Would work anyway, but since it is supposed to do nothing ... +!!$ ! call psb_ovrl(ty,sm%desc_data,info,& +!!$ ! & update=sm%prol,work=aux) +!!$ +!!$ +!!$ case(psb_sum_,psb_avg_) +!!$ ! +!!$ ! Update the overlap of ty +!!$ ! +!!$ call psb_ovrl(ty,sm%desc_data,info,& +!!$ & update=sm%prol,work=aux) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_ovrl' +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_prol_') +!!$ goto 9999 +!!$ end select +!!$ +!!$ case('T','C') +!!$ ! +!!$ ! With transpose, we have to do it here +!!$ ! +!!$ if (sm%restr == psb_halo_) then +!!$ call psb_ovrl(ty,sm%desc_data,info,& +!!$ & update=psb_sum_,work=aux) +!!$ if(info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ ch_err='psb_ovrl' +!!$ goto 9999 +!!$ end if +!!$ else if (sm%restr /= psb_none_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Invalid mld_sub_restr_') +!!$ goto 9999 +!!$ end if +!!$ +!!$ case default +!!$ info=psb_err_iarg_invalid_i_ +!!$ int_err(1)=6 +!!$ ch_err(2:2)=trans +!!$ goto 9999 +!!$ end select +!!$ end do +!!$ +!!$ if (info /= psb_success_) then +!!$ info=psb_err_internal_error_ +!!$ call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') +!!$ goto 9999 +!!$ end if +!!$ +!!$ +!!$ else +!!$ +!!$ info = psb_err_iarg_neg_ +!!$ call psb_errpush(info,name,& +!!$ & i_err=(/2,sweeps,0,0,0/)) +!!$ goto 9999 +!!$ +!!$ +!!$ end if +!!$ +!!$ ! +!!$ ! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) +!!$ ! +!!$ call psb_geaxpby(alpha,ty,beta,y,desc_data,info) +!!$ + end if + + + if ((6*isz) <= size(work)) then + else if ((4*isz) <= size(work)) then + deallocate(ww,tx,ty) + else if ((3*isz) <= size(work)) then + deallocate(aux) + else + deallocate(ww,aux,tx,ty) + endif + + 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_as_smoother_apply_vect + + subroutine d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) use psb_base_mod type(psb_desc_type), intent(in) :: desc_data @@ -178,7 +615,8 @@ contains trans_ = psb_toupper(trans) select case(trans_) case('N') - case('T','C') + case('T') + case('C') case default call psb_errpush(psb_err_iarg_invalid_i_,name) goto 9999 @@ -585,19 +1023,21 @@ contains end subroutine d_as_smoother_apply - subroutine d_as_smoother_bld(a,desc_a,sm,upd,info,mold) + subroutine d_as_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_as_smoother_type), intent(inout) :: sm - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_as_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold + ! Local variables type(psb_dspmat_type) :: blck, atmp integer :: n_row,n_col, nrow_a, nhalo, novr, data_, nzeros @@ -688,7 +1128,8 @@ contains End if if (info == psb_success_) & - & call sm%sv%build(a,sm%desc_data,upd,info,blck,mold=mold) + & call sm%sv%build(a,sm%desc_data,upd,info,& + & blck,amold=amold,vmold=vmold) nrow_a = a%get_nrows() n_row = sm%desc_data%get_local_rows() @@ -699,9 +1140,16 @@ contains if (info == psb_success_) call blck%csclip(atmp,info,& & jmin=nrow_a+1,rscale=.false.,cscale=.false.) if (info == psb_success_) call psb_rwextd(n_row,sm%nd,info,b=atmp) - if (info == psb_success_) call sm%nd%cscnv(info,& - & type='csr',dupl=psb_dupl_add_) - + + if (info == psb_success_) then + if (present(amold)) then + call sm%nd%cscnv(info,& + & mold=amold,dupl=psb_dupl_add_) + else + call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + end if + end if if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='clip & psb_spcnv csr 4') goto 9999 diff --git a/mlprec/mld_d_diag_solver.f90 b/mlprec/mld_d_diag_solver.f90 index 7b44be86..5c17657a 100644 --- a/mlprec/mld_d_diag_solver.f90 +++ b/mlprec/mld_d_diag_solver.f90 @@ -48,9 +48,11 @@ module mld_d_diag_solver use mld_d_prec_type type, extends(mld_d_base_solver_type) :: mld_d_diag_solver_type - real(psb_dpk_), allocatable :: d(:) + type(psb_d_vect_type), allocatable :: dv + real(psb_dpk_), allocatable :: d(:) contains procedure, pass(sv) :: build => d_diag_solver_bld + procedure, pass(sv) :: apply_v => d_diag_solver_apply_vect procedure, pass(sv) :: apply_a => d_diag_solver_apply procedure, pass(sv) :: free => d_diag_solver_free procedure, pass(sv) :: seti => d_diag_solver_seti @@ -69,16 +71,89 @@ module mld_d_diag_solver contains + subroutine d_diag_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_diag_solver_type), intent(inout) :: sv + type(psb_d_vect_type), intent(inout) :: x + type(psb_d_vect_type), intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_diag_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T','C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/)) + goto 9999 + end if + if (.not.allocated(sv%dv)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + if (sv%dv%get_nrows() < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: D") + goto 9999 + end if + + call y%mlt(alpha,sv%dv,x,beta,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='vect%mlt') + 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_diag_solver_apply_vect + subroutine d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data + type(psb_desc_type), intent(in) :: desc_data class(mld_d_diag_solver_type), intent(in) :: sv - real(psb_dpk_),intent(inout) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info + real(psb_dpk_), intent(inout) :: x(:) + real(psb_dpk_), intent(inout) :: y(:) + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info integer :: n_row,n_col real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) @@ -189,20 +264,21 @@ contains end subroutine d_diag_solver_apply - subroutine d_diag_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_diag_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_diag_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target, optional :: b - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_diag_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables integer :: n_row,n_col, nrow_a, nztota real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) @@ -251,7 +327,20 @@ contains sv%d(i) = done/sv%d(i) end if end do - + allocate(sv%dv,stat=info) + if (info == psb_success_) then + if (present(vmold)) then + allocate(sv%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: sv%dv%v,stat=info) + end if + end if + if (info == psb_success_) then + call sv%dv%bld(sv%d) + else + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate sv%dv') + goto 9999 + end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' end' diff --git a/mlprec/mld_d_id_solver.f90 b/mlprec/mld_d_id_solver.f90 index 64c46400..68493ee1 100644 --- a/mlprec/mld_d_id_solver.f90 +++ b/mlprec/mld_d_id_solver.f90 @@ -50,6 +50,7 @@ module mld_d_id_solver type, extends(mld_d_base_solver_type) :: mld_d_id_solver_type contains procedure, pass(sv) :: build => d_id_solver_bld + procedure, pass(sv) :: apply_v => d_id_solver_apply_vect procedure, pass(sv) :: apply_a => d_id_solver_apply procedure, pass(sv) :: free => d_id_solver_free procedure, pass(sv) :: seti => d_id_solver_seti @@ -68,6 +69,53 @@ module mld_d_id_solver contains + subroutine d_id_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_id_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_id_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + call psb_geaxpby(alpha,x,beta,y,desc_data,info) + + 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_id_solver_apply_vect + + subroutine d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use psb_base_mod type(psb_desc_type), intent(in) :: desc_data @@ -92,7 +140,8 @@ contains trans_ = psb_toupper(trans) select case(trans_) case('N') - case('T','C') + case('T') + case('C') case default call psb_errpush(psb_err_iarg_invalid_i_,name) goto 9999 @@ -113,20 +162,21 @@ contains end subroutine d_id_solver_apply - subroutine d_id_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_id_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_id_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold - type(psb_dspmat_type), intent(in), target, optional :: b + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_id_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables integer :: n_row,n_col, nrow_a, nztota real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) diff --git a/mlprec/mld_d_ilu_solver.f90 b/mlprec/mld_d_ilu_solver.f90 index 1f17bb69..702dcf2e 100644 --- a/mlprec/mld_d_ilu_solver.f90 +++ b/mlprec/mld_d_ilu_solver.f90 @@ -51,11 +51,13 @@ module mld_d_ilu_solver type, extends(mld_d_base_solver_type) :: mld_d_ilu_solver_type type(psb_dspmat_type) :: l, u real(psb_dpk_), allocatable :: d(:) + type(psb_d_vect_type), allocatable :: dv integer :: fact_type, fill_in real(psb_dpk_) :: thresh contains procedure, pass(sv) :: dump => d_ilu_solver_dmp procedure, pass(sv) :: build => d_ilu_solver_bld + procedure, pass(sv) :: apply_v => d_ilu_solver_apply_vect procedure, pass(sv) :: apply_a => d_ilu_solver_apply procedure, pass(sv) :: free => d_ilu_solver_free procedure, pass(sv) :: seti => d_ilu_solver_seti @@ -63,6 +65,7 @@ module mld_d_ilu_solver procedure, pass(sv) :: setr => d_ilu_solver_setr procedure, pass(sv) :: descr => d_ilu_solver_descr procedure, pass(sv) :: sizeof => d_ilu_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_ilu_solver_get_nzeros procedure, pass(sv) :: default => d_ilu_solver_default end type mld_d_ilu_solver_type @@ -71,7 +74,8 @@ module mld_d_ilu_solver & d_ilu_solver_free, d_ilu_solver_seti, & & d_ilu_solver_setc, d_ilu_solver_setr,& & d_ilu_solver_descr, d_ilu_solver_sizeof, & - & d_ilu_solver_default, d_ilu_solver_dmp + & d_ilu_solver_default, d_ilu_solver_dmp, & + & d_ilu_solver_apply_vect, d_ilu_solver_get_nzeros character(len=15), parameter, private :: & @@ -143,6 +147,142 @@ contains end subroutine d_ilu_solver_check + subroutine d_ilu_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_ilu_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + type(psb_d_vect_type) :: wv + real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_ilu_solver_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/)) + goto 9999 + end if + if (y%get_nrows() < 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(sv%dv)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: DV") + goto 9999 + end if + if (sv%dv%get_nrows() < n_row) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner: DV") + 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) + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + endif + if (info == psb_success_) allocate(wv%v,mold=x%v) + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + call wv%bld(n_col) + + + + select case(trans_) + case('N') + call psb_spsm(done,sv%l,x,dzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) + + if (info == psb_success_) call psb_spsm(alpha,sv%u,wv,beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_, work=aux) + + case('T') + call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) + if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case('C') + call psb_spsm(done,sv%u,x,dzero,wv,desc_data,info,& + & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) + if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve') + goto 9999 + endif + + 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_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_ilu_solver_apply_vect + + subroutine d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use psb_base_mod type(psb_desc_type), intent(in) :: desc_data @@ -207,7 +347,14 @@ contains if (info == psb_success_) call psb_spsm(alpha,sv%u,ww,beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - case('T','C') + case('T') + call psb_spsm(done,sv%u,x,dzero,ww,desc_data,info,& + & trans=trans_,scale='L',diag=sv%d,choice=psb_none_,work=aux) + if (info == psb_success_) call psb_spsm(alpha,sv%l,ww,beta,y,desc_data,info,& + & trans=trans_,scale='U',choice=psb_none_,work=aux) + case('C') + ! For real this is the same of T but in the future we will move to + ! a preprocessed version call psb_spsm(done,sv%u,x,dzero,ww,desc_data,info,& & trans=trans_,scale='L',diag=sv%d,choice=psb_none_,work=aux) if (info == psb_success_) call psb_spsm(alpha,sv%l,ww,beta,y,desc_data,info,& @@ -246,20 +393,21 @@ contains end subroutine d_ilu_solver_apply - subroutine d_ilu_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_ilu_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_ilu_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold - type(psb_dspmat_type), intent(in), target, optional :: b + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_ilu_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables integer :: n_row,n_col, nrow_a, nztota real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) @@ -301,11 +449,20 @@ contains endif if (.not.allocated(sv%d)) then allocate(sv%d(n_row),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + endif + if (info == psb_success_) then + allocate(sv%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(sv%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: sv%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 @@ -403,10 +560,11 @@ contains call sv%l%trim() call sv%u%set_asb() call sv%u%trim() + call sv%dv%bld(sv%d) - if (present(mold)) then - call sv%l%cscnv(info,mold=mold) - call sv%u%cscnv(info,mold=mold) + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) end if if (debug_level >= psb_debug_outer_) & @@ -626,6 +784,22 @@ contains return end subroutine d_ilu_solver_descr + function d_ilu_solver_get_nzeros(sv) result(val) + use psb_base_mod + implicit none + ! Arguments + class(mld_d_ilu_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer :: i + + val = 0 + if (allocated(sv%dv)) val = val + sv%dv%get_nrows() + val = val + sv%l%get_nzeros() + val = val + sv%u%get_nzeros() + + return + end function d_ilu_solver_get_nzeros + function d_ilu_solver_sizeof(sv) result(val) use psb_base_mod implicit none diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 24f5b855..e7ffe16c 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -50,16 +50,17 @@ module mld_d_inner_mod interface mld_mlprec_bld - subroutine mld_dmlprec_bld(a,desc_a,prec,info, mold) + subroutine mld_dmlprec_bld(a,desc_a,prec,info, amold, vmold) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, & - & psb_dpk_, psb_d_base_sparse_mat + & psb_dpk_, psb_d_base_sparse_mat, psb_d_base_vect_type use mld_d_prec_type, only : mld_dprec_type implicit none - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dprec_type), intent(inout), target :: prec - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dprec_type), intent(inout), target :: prec + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold !!$ character, intent(in),optional :: upd end subroutine mld_dmlprec_bld end interface mld_mlprec_bld diff --git a/mlprec/mld_d_jac_smoother.f90 b/mlprec/mld_d_jac_smoother.f90 index bca437c6..5b99f280 100644 --- a/mlprec/mld_d_jac_smoother.f90 +++ b/mlprec/mld_d_jac_smoother.f90 @@ -54,26 +54,188 @@ module mld_d_jac_smoother type(psb_dspmat_type) :: nd integer :: nnz_nd_tot contains - procedure, pass(sm) :: build => d_jac_smoother_bld - procedure, pass(sm) :: apply => d_jac_smoother_apply - procedure, pass(sm) :: free => d_jac_smoother_free - procedure, pass(sm) :: seti => d_jac_smoother_seti - procedure, pass(sm) :: setc => d_jac_smoother_setc - procedure, pass(sm) :: setr => d_jac_smoother_setr - procedure, pass(sm) :: descr => d_jac_smoother_descr - procedure, pass(sm) :: sizeof => d_jac_smoother_sizeof + procedure, pass(sm) :: build => d_jac_smoother_bld + procedure, pass(sm) :: apply_v => d_jac_smoother_apply_vect + procedure, pass(sm) :: apply_a => d_jac_smoother_apply + procedure, pass(sm) :: free => d_jac_smoother_free + procedure, pass(sm) :: seti => d_jac_smoother_seti + procedure, pass(sm) :: setc => d_jac_smoother_setc + procedure, pass(sm) :: setr => d_jac_smoother_setr + procedure, pass(sm) :: descr => d_jac_smoother_descr + procedure, pass(sm) :: sizeof => d_jac_smoother_sizeof end type mld_d_jac_smoother_type private :: d_jac_smoother_bld, d_jac_smoother_apply, & & d_jac_smoother_free, d_jac_smoother_seti, & & d_jac_smoother_setc, d_jac_smoother_setr,& - & d_jac_smoother_descr, d_jac_smoother_sizeof + & d_jac_smoother_descr, d_jac_smoother_sizeof, & + & d_jac_smoother_apply_vect contains + + subroutine d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_jac_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + integer :: n_row,n_col + type(psb_d_vect_type) :: tx, ty + real(psb_dpk_), pointer :: ww(:), aux(:) + integer :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='d_jac_smoother_apply' + + call psb_erractionsave(err_act) + + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') + case('T') + case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + if (.not.allocated(sm%sv)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + 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 /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/4*n_col,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/5*n_col,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + endif + + if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then + + call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,& + & name,a_err='Error in sub_aply Jacobi Sweeps = 1') + goto 9999 + endif + + else if (sweeps > 1) then + + ! + ! + ! Apply multiple sweeps of a block-Jacobi solver + ! to compute an approximate solution of a linear system. + ! + ! + allocate(tx%v,mold=x%v,stat=info) + if (info == psb_success_) allocate(ty%v,mold=x%v,stat=info) + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*n_col,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + call tx%bld(x%get_nrows()) + call ty%bld(x%get_nrows()) + + do i=1, sweeps + ! + ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the + ! block diagonal part and the remaining part of the local matrix + ! and Y(j) is the approximate solution at sweep j. + ! + call psb_geaxpby(done,x,dzero,ty,desc_data,info) + call psb_spmm(-done,sm%nd,tx,done,ty,desc_data,info,work=aux,trans=trans_) + + if (info /= psb_success_) exit + + call sm%sv%apply(done,ty,dzero,tx,desc_data,trans_,aux,info) + + if (info /= psb_success_) exit + end do + + if (info == psb_success_) call psb_geaxpby(alpha,tx,beta,y,desc_data,info) + + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='subsolve with Jacobi sweeps > 1') + goto 9999 + end if + + call tx%free(info) + if (info == psb_success_) call ty%free(info) + if (info /= psb_success_) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='final cleanup with Jacobi sweeps > 1') + goto 9999 + end if + + else + + info = psb_err_iarg_neg_ + call psb_errpush(info,name,& + & i_err=(/2,sweeps,0,0,0/)) + goto 9999 + + endif + + + 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_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine d_jac_smoother_apply_vect + subroutine d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) use psb_base_mod type(psb_desc_type), intent(in) :: desc_data @@ -228,19 +390,20 @@ contains end subroutine d_jac_smoother_apply - subroutine d_jac_smoother_bld(a,desc_a,sm,upd,info,mold) + subroutine d_jac_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) use psb_base_mod use mld_d_diag_solver Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_jac_smoother_type), intent(inout) :: sm - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_jac_smoother_type), intent(inout) :: sm + character, intent(in) :: upd + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables integer :: n_row,n_col, nrow_a, nztota, nzeros real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) @@ -268,15 +431,22 @@ contains call a%csclip(sm%nd,info,& & jmin=nrow_a+1,rscale=.false.,cscale=.false.) end select - if (info == psb_success_) call sm%nd%cscnv(info,& - & type='csr',dupl=psb_dupl_add_) + if (info == psb_success_) then + if (present(amold)) then + call sm%nd%cscnv(info,& + & mold=amold,dupl=psb_dupl_add_) + else + call sm%nd%cscnv(info,& + & type='csr',dupl=psb_dupl_add_) + endif + end if if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='clip & psb_spcnv csr 4') goto 9999 end if - call sm%sv%build(a,desc_a,upd,info,mold=mold) + call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,& & a_err='solver build') diff --git a/mlprec/mld_d_prec_mod.f90 b/mlprec/mld_d_prec_mod.f90 index 8594e092..dc07d171 100644 --- a/mlprec/mld_d_prec_mod.f90 +++ b/mlprec/mld_d_prec_mod.f90 @@ -111,16 +111,17 @@ module mld_d_prec_mod end interface interface mld_precbld - subroutine mld_dprecbld(a,desc_a,prec,info,mold) + subroutine mld_dprecbld(a,desc_a,prec,info,amold,vmold) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, & - & psb_dpk_, psb_d_base_sparse_mat + & psb_dpk_, psb_d_base_sparse_mat, psb_d_base_vect_type use mld_d_prec_type, only : mld_dprec_type implicit none - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dprec_type), intent(inout), target :: prec - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dprec_type), intent(inout), target :: prec + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold !!$ character, intent(in),optional :: upd end subroutine mld_dprecbld end interface diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index 000bf1df..f7d9e9b5 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -60,6 +60,7 @@ module mld_d_prec_type use mld_base_prec_type + use psb_base_mod, only : psb_d_vect_type, psb_d_base_vect_type ! ! Type: mld_Tprec_type. ! @@ -191,6 +192,7 @@ module mld_d_prec_type procedure, pass(sv) :: default => d_base_solver_default procedure, pass(sv) :: descr => d_base_solver_descr procedure, pass(sv) :: sizeof => d_base_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_base_solver_get_nzeros end type mld_d_base_solver_type type mld_d_base_smoother_type @@ -199,7 +201,9 @@ module mld_d_prec_type procedure, pass(sm) :: check => d_base_smoother_check procedure, pass(sm) :: dump => d_base_smoother_dmp procedure, pass(sm) :: build => d_base_smoother_bld - procedure, pass(sm) :: apply => d_base_smoother_apply + procedure, pass(sm) :: apply_v => d_base_smoother_apply_vect + procedure, pass(sm) :: apply_a => d_base_smoother_apply + generic, public :: apply => apply_a, apply_v procedure, pass(sm) :: free => d_base_smoother_free procedure, pass(sm) :: seti => d_base_smoother_seti procedure, pass(sm) :: setc => d_base_smoother_setc @@ -208,6 +212,7 @@ module mld_d_prec_type procedure, pass(sm) :: default => d_base_smoother_default procedure, pass(sm) :: descr => d_base_smoother_descr procedure, pass(sm) :: sizeof => d_base_smoother_sizeof + procedure, pass(sm) :: get_nzeros => d_base_smoother_get_nzeros end type mld_d_base_smoother_type type mld_donelev_type @@ -227,6 +232,7 @@ module mld_d_prec_type procedure, pass(lv) :: setr => d_base_onelev_setr procedure, pass(lv) :: setc => d_base_onelev_setc generic, public :: set => seti, setr, setc + procedure, pass(lv) :: get_nzeros => d_base_onelev_get_nzeros end type mld_donelev_type type, extends(psb_dprec_type) :: mld_dprec_type @@ -234,11 +240,13 @@ module mld_d_prec_type real(psb_dpk_) :: op_complexity=-done type(mld_donelev_type), allocatable :: precv(:) contains + procedure, pass(prec) :: d_apply2_vect => mld_d_apply2_vect procedure, pass(prec) :: d_apply2v => mld_d_apply2v procedure, pass(prec) :: d_apply1v => mld_d_apply1v procedure, pass(prec) :: dump => mld_d_dump procedure, pass(prec) :: get_complexity => mld_d_get_compl procedure, pass(prec) :: cmp_complexity => mld_d_cmp_compl + procedure, pass(prec) :: get_nzeros => mld_d_get_nzeros end type mld_dprec_type private :: d_base_solver_bld, d_base_solver_apply, & @@ -246,18 +254,20 @@ module mld_d_prec_type & d_base_solver_setc, d_base_solver_setr, & & d_base_solver_descr, d_base_solver_sizeof, & & d_base_solver_default, d_base_solver_check,& - & d_base_solver_dmp, & + & d_base_solver_dmp, d_base_solver_apply_vect, & & d_base_smoother_bld, d_base_smoother_apply, & & d_base_smoother_free, d_base_smoother_seti, & & d_base_smoother_setc, d_base_smoother_setr,& & d_base_smoother_descr, d_base_smoother_sizeof, & & d_base_smoother_default, d_base_smoother_check, & - & d_base_smoother_dmp, & + & d_base_smoother_dmp, d_base_smoother_apply_vect, & & d_base_onelev_seti, d_base_onelev_setc, & & d_base_onelev_setr, d_base_onelev_check, & & d_base_onelev_default, d_base_onelev_dump, & - & d_base_onelev_descr, mld_d_dump, & - & mld_d_get_compl, mld_d_cmp_compl + & d_base_onelev_descr, mld_d_dump, & + & mld_d_get_compl, mld_d_cmp_compl,& + & mld_d_get_nzeros, d_base_onelev_get_nzeros, & + & d_base_smoother_get_nzeros, d_base_solver_get_nzeros ! @@ -282,6 +292,18 @@ module mld_d_prec_type end interface interface mld_precaply + subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work) + use psb_base_mod, only : psb_dspmat_type, psb_desc_type, & + & psb_dpk_, psb_d_vect_type + import mld_dprec_type + type(psb_desc_type),intent(in) :: desc_data + type(mld_dprec_type), intent(inout) :: prec + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + end subroutine mld_dprecaply_vect subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ import mld_dprec_type @@ -309,6 +331,48 @@ contains ! Function returning the size of the mld_prec_type data structure ! + function d_base_solver_get_nzeros(sv) result(val) + implicit none + class(mld_d_base_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + end function d_base_solver_get_nzeros + + function d_base_smoother_get_nzeros(sm) result(val) + implicit none + class(mld_d_base_smoother_type), intent(in) :: sm + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(sm%sv)) & + & val = sm%sv%get_nzeros() + end function d_base_smoother_get_nzeros + + function d_base_onelev_get_nzeros(lv) result(val) + implicit none + class(mld_donelev_type), intent(in) :: lv + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(lv%sm)) & + & val = lv%sm%get_nzeros() + end function d_base_onelev_get_nzeros + + function mld_d_get_nzeros(prec) result(val) + implicit none + class(mld_dprec_type), intent(in) :: prec + integer(psb_long_int_k_) :: val + integer :: i + val = 0 + if (allocated(prec%precv)) then + do i=1, size(prec%precv) + val = val + prec%precv(i)%get_nzeros() + end do + end if + end function mld_d_get_nzeros + + function mld_dprec_sizeof(prec) result(val) implicit none type(mld_dprec_type), intent(in) :: prec @@ -451,6 +515,10 @@ contains endif call p%precv(1)%sm%descr(info,iout=iout_) if (nlev == 1) then + if (p%precv(1)%parms%sweeps > 1) then + write(iout_,*) ' Number of sweeps : ',& + & p%precv(1)%parms%sweeps + end if write(iout_,*) return end if @@ -688,6 +756,47 @@ contains end subroutine d_base_smoother_apply + subroutine d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& + & trans,sweeps,work,info) + use psb_base_mod + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_base_smoother_type), intent(inout) :: sm + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='d_base_smoother_apply' + + call psb_erractionsave(err_act) + info = psb_success_ + if (allocated(sm%sv)) then + call sm%sv%apply(alpha,x,beta,y,desc_data,trans,work,info) + else + info = 1121 + endif + if (info /= psb_success_) then + call psb_errpush(info,name) + 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_base_smoother_apply_vect + subroutine d_base_smoother_check(sm,info) use psb_base_mod @@ -830,19 +939,20 @@ contains return end subroutine d_base_smoother_setr - subroutine d_base_smoother_bld(a,desc_a,sm,upd,info,mold) + subroutine d_base_smoother_bld(a,desc_a,sm,upd,info,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a + type(psb_dspmat_type), intent(in), target :: a Type(psb_desc_type), Intent(in) :: desc_a class(mld_d_base_smoother_type), intent(inout) :: sm character, intent(in) :: upd integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold Integer :: err_act character(len=20) :: name='d_base_smoother_bld' @@ -850,7 +960,7 @@ contains info = psb_success_ if (allocated(sm%sv)) then - call sm%sv%build(a,desc_a,upd,info,mold=mold) + call sm%sv%build(a,desc_a,upd,info,amold=amold,vmold=vmold) else info = 1121 call psb_errpush(info,name) @@ -989,7 +1099,6 @@ contains end subroutine d_base_smoother_default - subroutine d_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use psb_base_mod type(psb_desc_type), intent(in) :: desc_data @@ -1023,17 +1132,16 @@ contains end subroutine d_base_solver_apply - subroutine d_base_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) use psb_base_mod - type(psb_desc_type), intent(in) :: desc_data - class(mld_d_base_solver_type), intent(in) :: sv - type(psb_d_vect_type),intent(inout) :: x - type(psb_d_vect_type),intent(inout) :: y - real(psb_dpk_),intent(in) :: alpha,beta - character(len=1),intent(in) :: trans - real(psb_dpk_),target, intent(inout) :: work(:) - integer, intent(out) :: info + type(psb_desc_type), intent(in) :: desc_data + class(mld_d_base_solver_type), intent(inout) :: sv + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + real(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_dpk_),target, intent(inout) :: work(:) + integer, intent(out) :: info Integer :: err_act character(len=20) :: name='d_base_solver_apply' @@ -1057,20 +1165,21 @@ contains end subroutine d_base_solver_apply_vect - subroutine d_base_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_base_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_base_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target, optional :: b - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_base_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold Integer :: err_act character(len=20) :: name='d_base_solver_bld' @@ -1287,6 +1396,43 @@ contains return end subroutine d_base_solver_default + + subroutine mld_d_apply2_vect(prec,x,y,desc_data,info,trans,work) + use psb_base_mod + type(psb_desc_type),intent(in) :: desc_data + class(mld_dprec_type), intent(inout) :: prec + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + 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' + + call psb_erractionsave(err_act) + + select type(prec) + type is (mld_dprec_type) + call mld_precaply(prec,x,y,desc_data,info,trans,work) + class default + info = psb_err_missing_override_method_ + call psb_errpush(info,name) + 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 mld_d_apply2_vect + subroutine mld_d_apply2v(prec,x,y,desc_data,info,trans,work) use psb_base_mod diff --git a/mlprec/mld_d_slu_solver.f90 b/mlprec/mld_d_slu_solver.f90 index 45032a07..fdfd004d 100644 --- a/mlprec/mld_d_slu_solver.f90 +++ b/mlprec/mld_d_slu_solver.f90 @@ -185,20 +185,21 @@ contains end subroutine d_slu_solver_apply - subroutine d_slu_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_slu_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_slu_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold - type(psb_dspmat_type), intent(in), target, optional :: b + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_slu_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables type(psb_dspmat_type) :: atmp type(psb_d_csr_sparse_mat) :: acsr diff --git a/mlprec/mld_d_sludist_solver.f90 b/mlprec/mld_d_sludist_solver.f90 index cae8bcd4..663a7bc6 100644 --- a/mlprec/mld_d_sludist_solver.f90 +++ b/mlprec/mld_d_sludist_solver.f90 @@ -185,20 +185,21 @@ contains end subroutine d_sludist_solver_apply - subroutine d_sludist_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_sludist_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_sludist_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold - type(psb_dspmat_type), intent(in), target, optional :: b + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_sludist_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables type(psb_dspmat_type) :: atmp type(psb_d_csr_sparse_mat) :: acsr diff --git a/mlprec/mld_d_umf_solver.f90 b/mlprec/mld_d_umf_solver.f90 index 26d98ae3..db2e835b 100644 --- a/mlprec/mld_d_umf_solver.f90 +++ b/mlprec/mld_d_umf_solver.f90 @@ -185,20 +185,21 @@ contains end subroutine d_umf_solver_apply - subroutine d_umf_solver_bld(a,desc_a,sv,upd,info,b,mold) + subroutine d_umf_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold) use psb_base_mod Implicit None ! Arguments - type(psb_dspmat_type), intent(in), target :: a - Type(psb_desc_type), Intent(in) :: desc_a - class(mld_d_umf_solver_type), intent(inout) :: sv - character, intent(in) :: upd - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold - type(psb_dspmat_type), intent(in), target, optional :: b + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_d_umf_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target, optional :: b + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold ! Local variables type(psb_dspmat_type) :: atmp type(psb_d_csc_sparse_mat) :: acsc diff --git a/mlprec/mld_dmlprec_bld.f90 b/mlprec/mld_dmlprec_bld.f90 index 964aab7a..5d618e67 100644 --- a/mlprec/mld_dmlprec_bld.f90 +++ b/mlprec/mld_dmlprec_bld.f90 @@ -63,7 +63,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_dmlprec_bld(a,desc_a,p,info,mold) +subroutine mld_dmlprec_bld(a,desc_a,p,info,amold,vmold) use psb_base_mod use mld_d_inner_mod, mld_protect_name => mld_dmlprec_bld @@ -72,11 +72,12 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,mold) Implicit None ! Arguments - type(psb_dspmat_type),intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dprec_type),intent(inout),target :: p - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type),intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dprec_type),intent(inout),target :: p + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold !!$ character, intent(in), optional :: upd ! Local Variables @@ -309,10 +310,10 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info,mold) ! Test version for beginning of OO stuff. ! call p%precv(i)%sm%build(p%precv(i)%base_a,p%precv(i)%base_desc,& - & 'F',info,mold=mold) + & 'F',info,amold=amold,vmold=vmold) - if ((info == psb_success_).and.(i>1).and.(present(mold))) then - call psb_map_cscnv(p%precv(i)%map,info,mold=mold) + if ((info == psb_success_).and.(i>1).and.(present(amold))) then + call psb_map_cscnv(p%precv(i)%map,info,mold=amold) end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& diff --git a/mlprec/mld_dprecaply.f90 b/mlprec/mld_dprecaply.f90 index 3329ee0b..f691c227 100644 --- a/mlprec/mld_dprecaply.f90 +++ b/mlprec/mld_dprecaply.f90 @@ -263,3 +263,105 @@ subroutine mld_dprecaply1(prec,x,desc_data,info,trans) end if return end subroutine mld_dprecaply1 + + + +subroutine mld_dprecaply_vect(prec,x,y,desc_data,info,trans,work) + + use psb_base_mod + use mld_d_inner_mod + + implicit none + + ! Arguments + type(psb_desc_type),intent(in) :: desc_data + type(mld_dprec_type), intent(inout) :: prec + type(psb_d_vect_type),intent(inout) :: x + type(psb_d_vect_type),intent(inout) :: y + integer, intent(out) :: info + character(len=1), optional :: trans + real(psb_dpk_),intent(inout), optional, target :: work(:) + + ! Local variables + character :: trans_ + real(psb_dpk_), pointer :: work_(:) + integer :: ictxt,np,me,err_act,iwsz + character(len=20) :: name + + name='mld_dprecaply' + info = psb_success_ + call psb_erractionsave(err_act) + + 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 (present(work)) then + work_ => work + else + iwsz = max(1,4*psb_cd_get_local_cols(desc_data)) + allocate(work_(iwsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_request_,name,i_err=(/iwsz,0,0,0,0/),& + &a_err='real(psb_dpk_)') + goto 9999 + end if + + end if + + if (.not.(allocated(prec%precv))) then + !! Error 1: should call mld_dprecbld + info=3112 + call psb_errpush(info,name) + goto 9999 + end if + if (size(prec%precv) >1) then + ! + ! Number of levels > 1: apply the multilevel preconditioner + ! +!!$ call mld_mlprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info) + info = psb_err_missing_override_method_ + if(info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply') + goto 9999 + end if + + else if (size(prec%precv) == 1) then + ! + ! Number of levels = 1: apply the base preconditioner + ! + call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,& + & prec%precv(1)%parms%sweeps, work_,info) + else + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='Invalid size of precv',& + & i_Err=(/size(prec%precv),0,0,0,0/)) + goto 9999 + endif + + ! If the original distribution has an overlap we should fix that. +!!$ call psb_halo(y,desc_data,info,data=psb_comm_mov_) + + + if (present(work)) then + else + deallocate(work_) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_dprecaply_vect diff --git a/mlprec/mld_dprecbld.f90 b/mlprec/mld_dprecbld.f90 index 32d468de..68562e7f 100644 --- a/mlprec/mld_dprecbld.f90 +++ b/mlprec/mld_dprecbld.f90 @@ -58,7 +58,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_dprecbld(a,desc_a,p,info,mold) +subroutine mld_dprecbld(a,desc_a,p,info,amold,vmold) use psb_base_mod use mld_d_inner_mod @@ -67,11 +67,12 @@ subroutine mld_dprecbld(a,desc_a,p,info,mold) Implicit None ! Arguments - type(psb_dspmat_type),intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dprec_type),intent(inout),target :: p - integer, intent(out) :: info - class(psb_d_base_sparse_mat), intent(in), optional :: mold + type(psb_dspmat_type),intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dprec_type),intent(inout), target :: p + integer, intent(out) :: info + class(psb_d_base_sparse_mat), intent(in), optional :: amold + class(psb_d_base_vect_type), intent(in), optional :: vmold !!$ character, intent(in), optional :: upd ! Local Variables @@ -171,7 +172,7 @@ subroutine mld_dprecbld(a,desc_a,p,info,mold) goto 9999 endif - call p%precv(1)%sm%build(a,desc_a,upd_,info,mold=mold) + call p%precv(1)%sm%build(a,desc_a,upd_,info,amold=amold,vmold=vmold) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& & a_err='One level preconditioner build.') @@ -185,7 +186,7 @@ subroutine mld_dprecbld(a,desc_a,p,info,mold) ! ! Build the multilevel preconditioner ! - call mld_mlprec_bld(a,desc_a,p,info,mold=mold) + call mld_mlprec_bld(a,desc_a,p,info,amold=amold,vmold=vmold) if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,&