From 048453938bdc0fcec8e145548e835c986923b272 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 3 Oct 2006 15:35:14 +0000 Subject: [PATCH] Included first development version of multilevel stuff. --- src/modules/psb_const.fh | 2 +- src/modules/psb_prec_mod.f90 | 7 +- src/modules/psb_prec_type.f90 | 174 +++++-- src/modules/psb_serial_mod.f90 | 6 +- src/prec/psb_dbaseprc_aply.f90 | 85 ++-- src/prec/psb_dbaseprc_bld.f90 | 15 +- src/prec/psb_dbjac_aply.f90 | 43 +- src/prec/psb_dbldaggrmat.f90 | 41 +- src/prec/psb_ddiagsc_bld.f90 | 5 + src/prec/psb_dgenaggrmap.f90 | 2 +- src/prec/psb_dmlprc_aply.f90 | 834 +++++++++++++++++++++----------- src/prec/psb_dmlprc_bld.f90 | 14 +- src/prec/psb_dprc_aply.f90 | 12 +- src/prec/psb_dprecbld.f90 | 52 +- src/prec/psb_dprecset.f90 | 209 ++++---- src/prec/psb_zbaseprc_aply.f90 | 67 ++- src/prec/psb_zbaseprc_bld.f90 | 15 +- src/prec/psb_zbjac_aply.f90 | 41 +- src/prec/psb_zbldaggrmat.f90 | 21 +- src/prec/psb_zmlprc_aply.f90 | 845 ++++++++++++++++++++++----------- src/prec/psb_zmlprc_bld.f90 | 6 +- src/prec/psb_zprc_aply.f90 | 12 +- src/prec/psb_zprecbld.f90 | 40 +- src/serial/dp/check_dim.f | 2 +- src/serial/f77/smmp.f | 19 +- src/serial/jad/djdnrmi.f | 2 +- src/serial/psb_cest.f90 | 4 +- src/serial/psb_dcsprt.f90 | 6 +- src/serial/psb_dsymbmm.f90 | 38 +- src/serial/psb_zsymbmm.f90 | 33 +- src/tools/psb_cdalv.f90 | 1 - test/Fileread/df_sample.f90 | 1 - 32 files changed, 1653 insertions(+), 1001 deletions(-) diff --git a/src/modules/psb_const.fh b/src/modules/psb_const.fh index 97030853..8d7f56d0 100644 --- a/src/modules/psb_const.fh +++ b/src/modules/psb_const.fh @@ -124,4 +124,4 @@ real(kind(1.d0)), parameter :: epstol=1.d-32 character, parameter :: psb_all_='A', psb_topdef_=' ' - character(len=5) :: psb_fidef_='CSR' + character(len=5) :: psb_fidef_='CSR' diff --git a/src/modules/psb_prec_mod.f90 b/src/modules/psb_prec_mod.f90 index 9d7d9d06..c7b5ee66 100644 --- a/src/modules/psb_prec_mod.f90 +++ b/src/modules/psb_prec_mod.f90 @@ -39,7 +39,7 @@ module psb_prec_mod use psb_prec_type implicit none type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(in), target :: desc_a type(psb_dprec_type), intent(inout) :: prec integer, intent(out) :: info character, intent(in),optional :: upd @@ -49,7 +49,7 @@ module psb_prec_mod use psb_prec_type implicit none type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(in), target :: desc_a type(psb_zprec_type), intent(inout) :: prec integer, intent(out) :: info character, intent(in),optional :: upd @@ -57,7 +57,7 @@ module psb_prec_mod end interface interface psb_precset - subroutine psb_dprecset(prec,ptype,info,iv,rs,rv) + subroutine psb_dprecset(prec,ptype,info,iv,rs,rv,ilev,nlev) use psb_serial_mod use psb_descriptor_type use psb_prec_type @@ -66,6 +66,7 @@ module psb_prec_mod character(len=*), intent(in) :: ptype integer, intent(out) :: info integer, optional, intent(in) :: iv(:) + integer, optional, intent(in) :: nlev,ilev real(kind(1.d0)), optional, intent(in) :: rs real(kind(1.d0)), optional, intent(in) :: rv(:) end subroutine psb_dprecset diff --git a/src/modules/psb_prec_type.f90 b/src/modules/psb_prec_type.f90 index 9c50a674..c8842ef3 100644 --- a/src/modules/psb_prec_type.f90 +++ b/src/modules/psb_prec_type.f90 @@ -86,16 +86,60 @@ module psb_prec_type type(psb_dspmat_type), pointer :: av(:) => null() ! real(kind(1.d0)), pointer :: d(:) => null() - type(psb_desc_type), pointer :: desc_data => null() ! + type(psb_desc_type), pointer :: desc_data => null(), desc_ac=>null()! ! integer, pointer :: iprcparm(:) => null() ! real(kind(1.d0)), pointer :: dprcparm(:) => null() ! integer, pointer :: perm(:) => null(), invperm(:) => null() integer, pointer :: mlia(:) => null(), nlaggr(:) => null() ! - type(psb_dspmat_type), pointer :: aorig => null() ! + type(psb_dspmat_type), pointer :: base_a => null() ! + type(psb_desc_type), pointer :: base_desc => null() ! real(kind(1.d0)), pointer :: dorig(:) => null() ! end type psb_dbaseprc_type + + ! + ! Multilevel preconditioning + ! + ! To each level I there corresponds a matrix A(I) and a preconditioner K(I) + ! + ! A notational difference: in the DD reference above the preconditioner for + ! a given level K(I) is written out as a sum over the subdomains + ! + ! SUM_k(R_k^T A_k R_k) + ! + ! whereas in this code the sum is implicit in the parallelization, + ! i.e. each process takes care of one subdomain, and for each level we have + ! as many subdomains as there are processes (except for the coarsest level where + ! we might have a replicated index space). Thus the sum apparently disappears + ! from our code, but only apparently, because it is implicit in the call + ! to psb_baseprc_aply. + ! + ! A bit of description of the baseprecv(:) data structure: + ! 1. Number of levels = NLEV = size(baseprecv(:)) + ! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level. + ! Includes: + ! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners + ! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners + ! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps + ! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV + ! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors + ! (ilev-1) ---> (ilev) + ! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors + ! (ilev) ---> (ilev-1) + ! Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe + ! + ! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV + ! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix + ! of the current level, i.e.: if ILEV=1 then A + ! else the aggregated matrix av(ac_); so we have + ! a unified treatment of residuals. Need this to + ! avoid passing explicitly matrix A to the + ! outer prec. routine + ! 5. baseprecv(ilev)%mlia The aggregation map from (ilev-1)-->(ilev) + ! if no smoother, it is used instead of sm_pr_ + ! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs. + ! type psb_dprec_type type(psb_dbaseprc_type), pointer :: baseprecv(:) => null() ! contain type of preconditioning to be performed @@ -105,14 +149,15 @@ module psb_prec_type type psb_zbaseprc_type type(psb_zspmat_type), pointer :: av(:) => null() ! - complex(kind(1.d0)), pointer :: d(:) => null() - type(psb_desc_type), pointer :: desc_data => null() ! + complex(kind(1.d0)), pointer :: d(:) => null() + type(psb_desc_type), pointer :: desc_data => null() , desc_ac=>null()! integer, pointer :: iprcparm(:) => null() ! real(kind(1.d0)), pointer :: dprcparm(:) => null() ! integer, pointer :: perm(:) => null(), invperm(:) => null() integer, pointer :: mlia(:) => null(), nlaggr(:) => null() ! - type(psb_zspmat_type), pointer :: aorig => null() ! - complex(kind(1.d0)), pointer :: dorig(:) => null() ! + type(psb_zspmat_type), pointer :: base_a => null() ! + type(psb_desc_type), pointer :: base_desc => null() ! + complex(kind(1.d0)), pointer :: dorig(:) => null() ! end type psb_zbaseprc_type @@ -181,6 +226,7 @@ contains subroutine psb_file_prec_descr(iout,p) integer, intent(in) :: iout type(psb_dprec_type), intent(in) :: p + integer :: ilev write(iout,*) 'Preconditioner description' if (associated(p%baseprecv)) then @@ -206,40 +252,44 @@ contains end select end if if (size(p%baseprecv)>=2) then - if (.not.associated(p%baseprecv(2)%iprcparm)) then - write(iout,*) 'Inconsistent MLPREC part!' - return - endif - write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_)) - if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then - write(iout,*) 'Multilevel aggregation: ', & - & aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_)) - write(iout,*) 'Smoother: ', & - & smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_)) - if (p%baseprecv(2)%iprcparm(smth_kind_) /= no_smth_) then - write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_) - write(iout,*) 'Smoothing position: ',& - & smooth_names(p%baseprecv(2)%iprcparm(smth_pos_)) + do ilev = 2, size(p%baseprecv) + if (.not.associated(p%baseprecv(ilev)%iprcparm)) then + write(iout,*) 'Inconsistent MLPREC part!' + return + endif + + write(iout,*) 'Multilevel: Level No', ilev + write(iout,*) 'Multilevel type: ',& + & ml_names(p%baseprecv(ilev)%iprcparm(ml_type_)) + if (p%baseprecv(ilev)%iprcparm(ml_type_)>no_ml_) then + write(iout,*) 'Multilevel aggregation: ', & + & aggr_names(p%baseprecv(ilev)%iprcparm(aggr_alg_)) + write(iout,*) 'Smoother: ', & + & smooth_kinds(p%baseprecv(ilev)%iprcparm(smth_kind_)) + if (p%baseprecv(ilev)%iprcparm(smth_kind_) /= no_smth_) then + write(iout,*) 'Smoothing omega: ', p%baseprecv(ilev)%dprcparm(smooth_omega_) + write(iout,*) 'Smoothing position: ',& + & smooth_names(p%baseprecv(ilev)%iprcparm(smth_pos_)) + end if + write(iout,*) 'Coarse matrix: ',& + & matrix_names(p%baseprecv(ilev)%iprcparm(coarse_mat_)) + write(iout,*) 'Aggregation sizes: ', & + & sum( p%baseprecv(ilev)%nlaggr(:)),' : ',p%baseprecv(ilev)%nlaggr(:) + write(iout,*) 'Factorization type: ',& + & fact_names(p%baseprecv(ilev)%iprcparm(f_type_)) + select case(p%baseprecv(ilev)%iprcparm(f_type_)) + case(f_ilu_n_) + write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(ilu_fill_in_) + case(f_ilu_e_) + write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(fact_eps_) + case(f_slu_,f_umf_) + case default + write(iout,*) 'Should never get here!' + end select + write(iout,*) 'Number of Jacobi sweeps: ', & + & (p%baseprecv(ilev)%iprcparm(jac_sweeps_)) end if - write(iout,*) 'Coarse matrix: ',& - & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_)) - write(iout,*) 'Aggregation sizes: ', & - & sum( p%baseprecv(2)%nlaggr(:)),' : ',p%baseprecv(2)%nlaggr(:) - write(iout,*) 'Factorization type: ',& - & fact_names(p%baseprecv(2)%iprcparm(f_type_)) - select case(p%baseprecv(2)%iprcparm(f_type_)) - case(f_ilu_n_) - write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(ilu_fill_in_) - case(f_ilu_e_) - write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_) - case(f_slu_,f_umf_) - case default - write(iout,*) 'Should never get here!' - end select - write(iout,*) 'Number of Jacobi sweeps: ', & - & (p%baseprecv(2)%iprcparm(jac_sweeps_)) - - end if + end do end if else @@ -355,13 +405,14 @@ contains & aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_)) write(iout,*) 'Smoother: ', & & smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_)) - write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_) if (p%baseprecv(2)%iprcparm(smth_kind_) /= no_smth_) then + write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_) write(iout,*) 'Smoothing position: ',& & smooth_names(p%baseprecv(2)%iprcparm(smth_pos_)) - write(iout,*) 'Coarse matrix: ',& - & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_)) end if + + write(iout,*) 'Coarse matrix: ',& + & matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_)) write(iout,*) 'Aggregation sizes: ', & & sum( p%baseprecv(2)%nlaggr(:)),' : ',p%baseprecv(2)%nlaggr(:) write(iout,*) 'Factorization type: ',& @@ -631,12 +682,23 @@ contains end if deallocate(p%desc_data) endif + if (associated(p%desc_ac)) then + if (associated(p%desc_ac%matrix_data)) then + call psb_cdfree(p%desc_ac,info) + end if + deallocate(p%desc_ac) + endif + if (associated(p%dprcparm)) then deallocate(p%dprcparm,stat=info) end if - if (associated(p%aorig)) then + if (associated(p%base_a)) then + ! This is a pointer to something else, must not free it here. + nullify(p%base_a) + endif + if (associated(p%base_desc)) then ! This is a pointer to something else, must not free it here. - nullify(p%aorig) + nullify(p%base_desc) endif if (associated(p%dorig)) then deallocate(p%dorig,stat=info) @@ -676,7 +738,7 @@ contains type(psb_dbaseprc_type), intent(inout) :: p nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,& - & p%nlaggr,p%aorig,p%dorig,p%desc_data) + & p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac) end subroutine psb_nullify_dbaseprec @@ -712,12 +774,22 @@ contains end if deallocate(p%desc_data) endif + if (associated(p%desc_ac)) then + if (associated(p%desc_ac%matrix_data)) then + call psb_cdfree(p%desc_ac,info) + end if + deallocate(p%desc_ac) + endif if (associated(p%dprcparm)) then deallocate(p%dprcparm,stat=info) end if - if (associated(p%aorig)) then + if (associated(p%base_a)) then ! This is a pointer to something else, must not free it here. - nullify(p%aorig) + nullify(p%base_a) + endif + if (associated(p%base_desc)) then + ! This is a pointer to something else, must not free it here. + nullify(p%base_desc) endif if (associated(p%dorig)) then deallocate(p%dorig,stat=info) @@ -741,11 +813,11 @@ contains if (associated(p%iprcparm)) then if (p%iprcparm(f_type_)==f_slu_) then -!!$ call psb_zslu_free(p%iprcparm(slu_ptr_),info) + call psb_zslu_free(p%iprcparm(slu_ptr_),info) end if if (p%iprcparm(f_type_)==f_umf_) then -!!$ call psb_zumf_free(p%iprcparm(umf_symptr_),& -!!$ & p%iprcparm(umf_numptr_),info) + call psb_zumf_free(p%iprcparm(umf_symptr_),& + & p%iprcparm(umf_numptr_),info) end if deallocate(p%iprcparm,stat=info) end if @@ -757,7 +829,7 @@ contains type(psb_zbaseprc_type), intent(inout) :: p nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,& - & p%nlaggr,p%aorig,p%dorig,p%desc_data) + & p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data,p%desc_ac) end subroutine psb_nullify_zbaseprec diff --git a/src/modules/psb_serial_mod.f90 b/src/modules/psb_serial_mod.f90 index b58f5a88..58f55c0a 100644 --- a/src/modules/psb_serial_mod.f90 +++ b/src/modules/psb_serial_mod.f90 @@ -359,13 +359,15 @@ module psb_serial_mod interface psb_symbmm - subroutine psb_dsymbmm(a,b,c) + subroutine psb_dsymbmm(a,b,c,info) use psb_spmat_type type(psb_dspmat_type) :: a,b,c + integer :: info end subroutine psb_dsymbmm - subroutine psb_zsymbmm(a,b,c) + subroutine psb_zsymbmm(a,b,c,info) use psb_spmat_type type(psb_zspmat_type) :: a,b,c + integer :: info end subroutine psb_zsymbmm end interface diff --git a/src/prec/psb_dbaseprc_aply.f90 b/src/prec/psb_dbaseprc_aply.f90 index 8c3c09aa..cb16b786 100644 --- a/src/prec/psb_dbaseprc_aply.f90 +++ b/src/prec/psb_dbaseprc_aply.f90 @@ -34,9 +34,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) +subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! - ! Compute Y <- beta*Y + K^-1 X + ! Compute Y <- beta*Y + alpha*K^-1 X ! where K is a a basic preconditioner stored in prec ! @@ -52,7 +52,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) type(psb_desc_type),intent(in) :: desc_data type(psb_dbaseprc_type), intent(in) :: prec real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -68,17 +68,17 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) character(len=20) :: name, ch_err interface psb_bjac_aply - subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) - use psb_descriptor_type - use psb_prec_type - type(psb_desc_type), intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbjac_aply + subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type), intent(in) :: desc_data + type(psb_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(inout) :: x(:), y(:) + real(kind(0.d0)),intent(in) :: alpha,beta + character(len=1) :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dbjac_aply end interface name='psb_dbaseprc_aply' @@ -105,33 +105,35 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) case(noprec_) - n_row=desc_data%matrix_data(psb_n_row_) - if (beta==dzero) then - y(1:n_row) = x(1:n_row) - else if (beta==done) then - y(1:n_row) = x(1:n_row) + y(1:n_row) - else if (beta==-done) then - y(1:n_row) = x(1:n_row) - y(1:n_row) - else - y(1:n_row) = x(1:n_row) + beta*y(1:n_row) - end if + call psb_geaxpby(alpha,x,beta,y,desc_data,info) case(diagsc_) + + if (size(work) >= size(x)) then + ww => work + else + allocate(ww(size(x)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + end if n_row=desc_data%matrix_data(psb_n_row_) - if (beta==dzero) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - else if (beta==done) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row) - else if (beta==-done) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row) - else - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row) + ww(1:n_row) = x(1:n_row)*prec%d(1:n_row) + 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 case(bja_) - call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) + call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 ch_err='psb_bjac_aply' @@ -142,7 +144,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) if (prec%iprcparm(n_ovr_)==0) then ! shortcut: this fixes performance for RAS(0) == BJA - call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) + call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 ch_err='psb_bjacaply' @@ -214,7 +216,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) end if endif - call psb_bjac_aply(prec,tx,dzero,ty,prec%desc_data,trans,aux,info) + call psb_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info) if(info.ne.0) then info=4010 ch_err='psb_bjac_aply' @@ -250,18 +252,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) & prec%iprcparm(prol_) end select - if (beta == dzero) then - y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_)) - else if (beta == done) then - y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - else if (beta == -done) then - y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - else - y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) if ((6*isz) <= size(work)) then diff --git a/src/prec/psb_dbaseprc_bld.f90 b/src/prec/psb_dbaseprc_bld.f90 index d15da26c..735edc93 100644 --- a/src/prec/psb_dbaseprc_bld.f90 +++ b/src/prec/psb_dbaseprc_bld.f90 @@ -49,7 +49,7 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) Implicit None type(psb_dspmat_type), target :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(in), target :: desc_a type(psb_dbaseprc_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd @@ -118,7 +118,7 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) integer :: int_err(5) character :: iupd - logical, parameter :: debug=.false. + logical, parameter :: debug=.false. integer,parameter :: iroot=0,iout=60,ilout=40 character(len=20) :: name, ch_err @@ -169,7 +169,13 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) select case(p%iprcparm(p_type_)) case (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 (diagsc_) @@ -256,7 +262,8 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) end select - + p%base_a => a + p%base_desc => desc_a call psb_erractionrestore(err_act) return diff --git a/src/prec/psb_dbjac_aply.f90 b/src/prec/psb_dbjac_aply.f90 index b196a35a..c6e4e0e0 100644 --- a/src/prec/psb_dbjac_aply.f90 +++ b/src/prec/psb_dbjac_aply.f90 @@ -34,9 +34,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) +subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! - ! Compute Y <- beta*Y + K^-1 X + ! Compute Y <- beta*Y + alpha*K^-1 X ! where K is a a Block Jacobi preconditioner stored in prec ! Note that desc_data may or may not be the same as prec%desc_data, ! but since both are INTENT(IN) this should be legal. @@ -54,7 +54,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) type(psb_desc_type), intent(in) :: desc_data type(psb_dbaseprc_type), intent(in) :: prec real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -125,7 +125,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) & trans='N',unit=diagl,choice=psb_none_,work=aux) if(info /=0) goto 9999 ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(done,prec%av(u_pr_),ww,beta,y,desc_data,info,& + call psb_spsm(alpha,prec%av(u_pr_),ww,beta,y,desc_data,info,& & trans='N',unit=diagu,choice=psb_none_, work=aux) if(info /=0) goto 9999 @@ -134,7 +134,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) & trans=trans,unit=diagu,choice=psb_none_, work=aux) if(info /=0) goto 9999 ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(done,prec%av(l_pr_),ww,beta,y,desc_data,info,& + call psb_spsm(alpha,prec%av(l_pr_),ww,beta,y,desc_data,info,& & trans=trans,unit=diagl,choice=psb_none_,work=aux) if(info /=0) goto 9999 @@ -152,16 +152,8 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) end select if(info /=0) goto 9999 + call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - if (beta == dzero) then - y(1:n_row) = ww(1:n_row) - else if (beta==done) then - y(1:n_row) = ww(1:n_row) + y(1:n_row) - else if (beta==-done) then - y(1:n_row) = ww(1:n_row) - y(1:n_row) - else - y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) - endif case (f_umf_) @@ -174,15 +166,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) if(info /=0) goto 9999 - if (beta == dzero) then - y(1:n_row) = ww(1:n_row) - else if (beta==dzero) then - y(1:n_row) = ww(1:n_row) + y(1:n_row) - else if (beta==-dzero) then - y(1:n_row) = ww(1:n_row) - y(1:n_row) - else - y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) - endif + call psb_geaxpby(alpha,ww,beta,y,desc_data,info) case default write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_) @@ -252,16 +236,9 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) end do end select - - if (beta == dzero) then - y(1:n_row) = tx(1:n_row) - else if (beta==done) then - y(1:n_row) = tx(1:n_row) + y(1:n_row) - else if (beta==-done) then - y(1:n_row) = tx(1:n_row) - y(1:n_row) - else - y(1:n_row) = tx(1:n_row) + beta*y(1:n_row) - endif + + call psb_geaxpby(alpha,tx,beta,y,desc_data,info) + deallocate(tx,ty) diff --git a/src/prec/psb_dbldaggrmat.f90 b/src/prec/psb_dbldaggrmat.f90 index c5ad0211..5786b27c 100644 --- a/src/prec/psb_dbldaggrmat.f90 +++ b/src/prec/psb_dbldaggrmat.f90 @@ -75,13 +75,14 @@ subroutine psb_dbldaggrmat(a,desc_a,ac,p,desc_p,info) if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.') case(smth_omg_,smth_biz_) - + if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix') call smooth_aggregate(info) if(info /= 0) then call psb_errpush(4010,name,a_err='smooth_aggregate') goto 9999 end if + if (aggr_dump) call psb_csprt(90+me,ac,head='% Smooth aggregate.') case default call psb_errpush(4010,name,a_err=name) goto 9999 @@ -117,6 +118,7 @@ contains integer :: ictxt, nrow, nglob, ncol, ntaggr, nzbg, ip, ndx,& & naggr, np, me, nzt,irs,jl,nzl,nlr,& & icomm,naggrm1, i, j, k, err_act + name='raw_aggregate' if(psb_get_errstatus().ne.0) return info=0 @@ -228,11 +230,6 @@ contains enddo end if - call psb_fixcoo(b,info) - if(info /= 0) then - call psb_errpush(4010,name,a_err='fixcoo') - goto 9999 - end if irs = b%infoa(psb_nnz_) call psb_sp_reall(b,irs,info) @@ -523,8 +520,6 @@ contains endif - if (test_dump) call & - & psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob) call psb_ipcoo2csr(am4,info) @@ -542,7 +537,7 @@ contains ! ! WARNING: the cycles below assume that AM3 does have ! its diagonal elements stored explicitly!!! - ! Should we swicth to something safer? + ! Should we switch to something safer? ! call psb_spscal(am3,p%dorig,info) if(info /= 0) goto 9999 @@ -604,12 +599,15 @@ contains am3%aspk(j) = done - omega*am3%aspk(j) endif end do + call psb_ipcoo2csr(am3,info) else write(0,*) 'Missing implementation of I sum' call psb_errpush(4010,name) goto 9999 end if + if (test_dump) call & + & psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob) if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,& & ivc=desc_a%loc_to_glob) if (debug) write(0,*) me,'Done gather, going for SYMBMM 1' @@ -620,7 +618,15 @@ contains ! Doing it this way means to consider diag(Ai) ! ! - call psb_symbmm(am3,am4,am1) + call psb_symbmm(am3,am4,am1,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 1') + goto 9999 + end if + am1%aspk(:) = 0.d0 + if (test_dump) & + & call psb_csprt(50+me,am1,head='% (I-wDA)Pt symbmm ') + call psb_numbmm(am3,am4,am1) if (debug) write(0,*) me,'Done NUMBMM 1' @@ -667,7 +673,12 @@ contains if (test_dump) & & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) - call psb_symbmm(a,am1,am3) + call psb_symbmm(a,am1,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 2') + goto 9999 + end if + call psb_numbmm(a,am1,am3) if (debug) write(0,*) me,'Done NUMBMM 2' @@ -724,7 +735,12 @@ contains endif if (debug) write(0,*) me,'starting symbmm 3' - call psb_symbmm(am2,am3,b) + call psb_symbmm(am2,am3,b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm 3') + goto 9999 + end if + if (debug) write(0,*) me,'starting numbmm 3' call psb_numbmm(am2,am3,b) if (debug) write(0,*) me,'Done NUMBMM 3' @@ -1018,6 +1034,7 @@ contains deallocate(nzbr,idisp) end select + call psb_ipcoo2csr(bg,info) call psb_erractionrestore(err_act) return diff --git a/src/prec/psb_ddiagsc_bld.f90 b/src/prec/psb_ddiagsc_bld.f90 index 4ad4f15c..014a885e 100644 --- a/src/prec/psb_ddiagsc_bld.f90 +++ b/src/prec/psb_ddiagsc_bld.f90 @@ -98,6 +98,11 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if + call psb_cdcpy(desc_a,p%desc_Data,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_cdcpy') + goto 9999 + end if if (debug) write(ilout+me,*) 'VDIAG ',n_row do i=1,n_row diff --git a/src/prec/psb_dgenaggrmap.f90 b/src/prec/psb_dgenaggrmap.f90 index b9e6fa12..8c03c645 100644 --- a/src/prec/psb_dgenaggrmap.f90 +++ b/src/prec/psb_dgenaggrmap.f90 @@ -92,7 +92,7 @@ subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Very simple minded strategy. ! naggr = 0 - nlp = 0 + nlp = 0 do icnt = 0 do i=1, nr diff --git a/src/prec/psb_dmlprc_aply.f90 b/src/prec/psb_dmlprc_aply.f90 index 9fbf1a14..cede1a0c 100644 --- a/src/prec/psb_dmlprc_aply.f90 +++ b/src/prec/psb_dmlprc_aply.f90 @@ -1,4 +1,4 @@ -!!$ +!!$ !!$ !!$ MD2P4 !!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS @@ -34,11 +34,54 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) +subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) ! - ! Compute Y <- beta*Y + K^-1 X - ! where K is a multilevel (actually 2-level) preconditioner stored in prec + ! Compute Y <- beta*Y + alpha*K^-1 X + ! where K is a multilevel preconditioner stored in baseprecv ! + ! cfr.: Smith, Biorstad & Gropp + ! Domain Decomposition + ! Cambridge Univ. Press + ! + ! To each level I there corresponds a matrix A(I) and a preconditioner K(I) + ! + ! A notational difference: in the DD reference above the preconditioner for + ! a given level K(I) is written out as a sum over the subdomains + ! + ! SUM_k(R_k^T A_k R_k) + ! + ! whereas in this code the sum is implicit in the parallelization, + ! i.e. each process takes care of one subdomain, and for each level we have + ! as many subdomains as there are processes (except for the coarsest level where + ! we might have a replicated index space). Thus the sum apparently disappears + ! from our code, but only apparently, because it is implicit in the call + ! to psb_baseprc_aply. + ! + ! A bit of description of the baseprecv(:) data structure: + ! 1. Number of levels = NLEV = size(baseprecv(:)) + ! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level. + ! Includes: + ! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners + ! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners + ! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps + ! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV + ! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors + ! (ilev-1) ---> (ilev) + ! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors + ! (ilev) ---> (ilev-1) + ! Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe + ! + ! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV + ! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix + ! of the current level, i.e.: if ILEV=1 then A + ! else the aggregated matrix av(ac_); so we have + ! a unified treatment of residuals. Need this to + ! avoid passing explicitly matrix A to the + ! outer prec. routine + ! 5. baseprecv(ilev)%mlia The aggregation map from (ilev-1)-->(ilev) + ! if no smoother, it is used instead of sm_pr_ + ! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs. + ! use psb_serial_mod use psb_descriptor_type @@ -52,41 +95,49 @@ subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) type(psb_desc_type),intent(in) :: desc_data type(psb_dbaseprc_type), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(inout) :: x(:), y(:) character :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info + ! Local variables integer :: n_row,n_col - real(kind(1.d0)), allocatable :: tx(:),ty(:),t2l(:),w2l(:),& - & x2l(:),b2l(:),tz(:),tty(:) + real(kind(1.d0)), allocatable :: tx(:),ty(:),x2l(:),y2l(:),& + & b2l(:),tty(:) character ::diagl, diagu integer :: ictxt,np,me,i, isz, nrg,nr2l,err_act, iptype, int_err(5) real(kind(1.d0)) :: omega real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime logical, parameter :: debug=.false., debugprt=.false. - integer :: ismth + integer :: ismth, nlev, ilev external mpi_wtime character(len=20) :: name, ch_err + type psb_mlprec_wrk_type + real(kind(1.d0)), pointer :: tx(:)=>null(),ty(:)=>null(),& + & x2l(:)=>null(),y2l(:)=>null(),& + & b2l(:)=>null(),tty(:)=>null() + end type psb_mlprec_wrk_type + type(psb_mlprec_wrk_type), pointer :: mlprec_wrk(:) + interface psb_baseprc_aply - subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) + subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_descriptor_type use psb_prec_type type(psb_desc_type),intent(in) :: desc_data type(psb_dbaseprc_type), intent(in) :: prec real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info end subroutine psb_dbaseprc_aply end interface - name='psb_dmlprc_aply' + name='psb_mlprc_aply' info = 0 call psb_erractionsave(err_act) @@ -94,426 +145,625 @@ subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) ictxt=desc_data%matrix_data(psb_ctxt_) call psb_info(ictxt, me, np) - omega=baseprecv(2)%dprcparm(smooth_omega_) - ismth=baseprecv(2)%iprcparm(smth_kind_) + nlev = size(baseprecv) + allocate(mlprec_wrk(nlev),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + select case(baseprecv(2)%iprcparm(ml_type_)) + case(no_ml_) ! Should not really get here. - write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????' + call psb_errpush(4010,name,a_err='no_ml_ in mlprc_aply?') + goto 9999 + case(add_ml_prec_) - ! - ! Additive multilevel + ! - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - call psb_baseprc_aply(baseprecv(1),x,beta,y,desc_data,trans,work,info) - if(info /=0) goto 9999 - - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - t2l(:) = dzero - w2l(:) = dzero + ! Additive is very simple. + ! 1. X(1) = Xext + ! 2. DO ILEV=2,NLEV + ! X(ILEV) = AV(PR_SM_T_)*X(ILEV-1) + ! 3. Y(ILEV) = (K(ILEV)**(-1))*X(ILEV) + ! 4. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = AV(PR_SM_)*Y(ILEV+1) + ! 5. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + - if (ismth /= no_smth_) then - ! - ! Smoothed aggregation - ! - allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),& - & tz(max(n_row,n_col)),stat=info) + call psb_baseprc_aply(alpha,baseprecv(1),x,beta,y,& + & baseprecv(1)%base_desc,trans,work,info) + if(info /=0) goto 9999 + allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y))) + mlprec_wrk(1)%x2l(:) = x(:) + + + do ilev = 2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + allocate(mlprec_wrk(ilev)%x2l(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%tx(max(n_row,n_col)),& + & mlprec_wrk(ilev)%ty(max(n_row,n_col)), stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero - ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero - tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero + mlprec_wrk(ilev)%x2l(:) = dzero + mlprec_wrk(ilev)%y2l(:) = dzero + mlprec_wrk(ilev)%tx(1:n_row) = mlprec_wrk(ilev-1)%x2l(1:n_row) + mlprec_wrk(ilev)%tx(n_row+1:max(n_row,n_col)) = dzero + mlprec_wrk(ilev)%ty(:) = dzero + ismth=baseprecv(ilev)%iprcparm(smth_kind_) - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + + + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,& + & info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = dzero + end if + + call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l,& + & dzero,mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero - end if + ! + ! Raw aggregation, may take shortcut + ! + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%x2l(i) + end do - call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info) - if(info /=0) goto 9999 + end if - else - ! - ! Raw aggregation, may take shortcut - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i) - end do + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif - end if + call psb_baseprc_aply(done,baseprecv(ilev),& + & mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,& + & baseprecv(ilev)%desc_data, 'N',work,info) - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif + enddo - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) + do ilev =nlev,2,-1 + ismth=baseprecv(ilev)%iprcparm(smth_kind_) + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) - if (ismth /= no_smth_) then + if (ismth /= no_smth_) then - call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info) - if(info /=0) goto 9999 - ! - ! Finally add back into Y. - ! - call psb_geaxpby(done,ty,done,y,desc_data,info) - if(info /=0) goto 9999 - deallocate(tx,ty,tz) + call psb_csmm(done,baseprecv(ilev)%av(sm_pr_),mlprec_wrk(ilev)%y2l,& + & done,mlprec_wrk(ilev-1)%y2l,info) + if(info /=0) goto 9999 - else + else - do i=1, desc_data%matrix_data(psb_n_row_) - y(i) = y(i) + t2l(baseprecv(2)%mlia(i)) - enddo + do i=1, n_row + mlprec_wrk(ilev-1)%y2l(i) = mlprec_wrk(ilev-1)%y2l(i) + & + & mlprec_wrk(ilev)%y2l(baseprecv(ilev)%mlia(i)) + enddo - end if + end if + end do - if (debugprt) write(0,*)' Y2: ',Y(:) + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,done,y,baseprecv(1)%base_desc,info) + if(info /=0) goto 9999 - deallocate(t2l,w2l) case(mult_ml_prec_) ! ! Multiplicative multilevel ! Pre/post smoothing versions. + ! select case(baseprecv(2)%iprcparm(smth_pos_)) case(post_smooth_) + + ! + ! Post smoothing. + ! 1. X(1) = Xext + ! 2. DO ILEV=2, NLEV :: X(ILEV) = AV(PR_SM_T_,ILEV)*X(ILEV-1) + ! 3. Y(NLEV) = (K(NLEV)**(-1))*X(NLEV) + ! 4. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = AV(PR_SM_,ILEV+1)*Y(ILEV+1) + ! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV)) + ! + ! 5. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + ! + ! Also: post smoothing is not spelled out in detail in DD. + ! + ! + - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if + n_col = desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - t2l(:) = dzero - w2l(:) = dzero + allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), & + & mlprec_wrk(1)%tx(nr2l), stat=info) + mlprec_wrk(1)%x2l(:) = dzero + mlprec_wrk(1)%y2l(:) = dzero + mlprec_wrk(1)%tx(:) = dzero - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==dzero - ! + call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,& + & baseprecv(1)%base_desc,info) + call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,& + & baseprecv(1)%base_desc,info) + + do ilev=2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + ismth = baseprecv(ilev)%iprcparm(smth_kind_) + + allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%x2l(nr2l), stat=info) - if (debug) write(0,*)' mult_ml_apply omega ',omega - if (debugprt) write(0,*)' mult_ml_apply X: ',X(:) - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - if(info /=0) then - if (debug) write(0,*)' From axpby1 ',size(x),size(tx),n_row,n_col,nr2l,nrg - call psb_errpush(4010,name,a_err='axpby post_smooth 1') - goto 9999 - endif - if (ismth /= no_smth_) then - ! - ! Smoothed aggregation - ! - allocate(tz(max(n_row,n_col)),stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) + mlprec_wrk(ilev)%x2l(:) = dzero + mlprec_wrk(ilev)%y2l(:) = dzero + mlprec_wrk(ilev)%tx(:) = dzero + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + call psb_halo(mlprec_wrk(ilev-1)%x2l,& + & baseprecv(ilev-1)%base_desc,info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = dzero + end if + + call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l, & + & dzero,mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero + ! + ! Raw aggregation, may take shortcut + ! + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%x2l(i) + end do end if - call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info) + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif + call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,& + & baseprecv(ilev)%base_desc,info) if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcut - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if + enddo - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,& - &'N',work,info) + call psb_baseprc_aply(done,baseprecv(nlev),mlprec_wrk(nlev)%x2l, & + & dzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info) + if(info /=0) goto 9999 - if (ismth /= no_smth_) then - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info) - if(info /=0) goto 9999 - ! - ! Finally add back into Y. - ! - deallocate(tz) - else - ty(:) = dzero - do i=1, desc_data%matrix_data(psb_n_row_) - ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i)) - enddo + + do ilev=nlev-1, 1, -1 + ismth = baseprecv(ilev+1)%iprcparm(smth_kind_) + if (ismth /= no_smth_) then + if (ismth == smth_omg_) & + & call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,& + & info,work=work) + call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,& + & dzero,mlprec_wrk(ilev)%y2l,info) + if(info /=0) goto 9999 - end if - deallocate(t2l,w2l) + else + n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_) + mlprec_wrk(ilev)%y2l(:) = dzero + do i=1, n_row + mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & + & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) + enddo - call psb_spmm(-done,baseprecv(2)%aorig,ty,done,tx,desc_data,info,work=work) - if(info /=0) goto 9999 + end if - call psb_baseprc_aply(baseprecv(1),tx,done,ty,desc_data,trans,& - & work,info) - if(info /=0) goto 9999 + call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) - call psb_geaxpby(done,ty,beta,y,desc_data,info) - if(info /=0) goto 9999 + if(info /=0) goto 9999 + + call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,& + & done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info) + + if(info /=0) goto 9999 + + enddo - deallocate(tx,ty) + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info) + if(info /=0) goto 9999 case(pre_smooth_) - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) + + ! + ! Pre smoothing. + ! 1. X(1) = Xext + ! 2. Y(1) = (K(1)**(-1))*X(1) + ! 3. TX(1) = X(1) - A(1)*Y(1) + ! 4. DO ILEV=2, NLEV + ! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1) + ! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV) + ! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV)) + ! 5. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1) + ! 6. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + ! + ! + + n_col = desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + + allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), & + & mlprec_wrk(1)%tx(nr2l), stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') - goto 9999 + goto 9999 end if - t2l(:) = dzero - w2l(:) = dzero + mlprec_wrk(1)%y2l(:) = dzero - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - if(info /=0) goto 9999 - call psb_baseprc_aply(baseprecv(1),x,dzero,tty,desc_data,& + mlprec_wrk(1)%x2l(:) = x + + call psb_baseprc_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,& + & dzero,mlprec_wrk(1)%y2l,& + & baseprecv(1)%base_desc,& & trans,work,info) + if(info /=0) goto 9999 - call psb_spmm(-done,baseprecv(2)%aorig,tty,done,tx,desc_data,info,work=work) + mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l + + call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,& + & done,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work) if(info /=0) goto 9999 - if (ismth /= no_smth_) then - allocate(tz(max(n_row,n_col)),stat=info) + do ilev = 2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + ismth = baseprecv(ilev)%iprcparm(smth_kind_) + allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%x2l(nr2l), stat=info) + + if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if + mlprec_wrk(ilev)%x2l(:) = dzero + mlprec_wrk(ilev)%y2l(:) = dzero + mlprec_wrk(ilev)%tx(:) = dzero + - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) + if (ismth /= no_smth_) then + ! + !Smoothed Aggregation + ! + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + + call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,& + & info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%tx(n_row+1:max(n_row,n_col)) = dzero + end if + + call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%tx,dzero,& + & mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero + ! + ! Raw aggregation, may take shortcuts + ! + mlprec_wrk(ilev)%x2l = dzero + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%tx(i) + end do end if - call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info) + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif + + + call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,& + & dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info) + if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcuts - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if + if(ilev < nlev) then + mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l + call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) + if(info /=0) goto 9999 + endif - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif + enddo - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - if(info /=0) goto 9999 + do ilev = nlev-1, 1, -1 - if (ismth /= no_smth_) then + ismth=baseprecv(ilev+1)%iprcparm(smth_kind_) - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info) - if(info /=0) goto 9999 + if (ismth /= no_smth_) then - call psb_geaxpby(done,ty,done,tty,desc_data,info) - if(info /=0) goto 9999 + if (ismth == smth_omg_) & + & call psb_halo(mlprec_wrk(ilev+1)%y2l,& + & baseprecv(ilev+1)%desc_data,info,work=work) + call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,& + & done,mlprec_wrk(ilev)%y2l,info) - deallocate(tz) - else + if(info /=0) goto 9999 - do i=1, desc_data%matrix_data(psb_n_row_) - tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) - enddo + else - end if + n_row = baseprecv(ilev+1)%base_desc%matrix_data(psb_n_row_) + do i=1, n_row + mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & + & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) + enddo + + end if + + enddo + + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& + & baseprecv(1)%base_desc,info) - call psb_geaxpby(done,tty,beta,y,desc_data,info) if(info /=0) goto 9999 - deallocate(t2l,w2l,tx,ty,tty) case(smooth_both_) - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) + ! + ! Symmetrized smoothing. + ! 1. X(1) = Xext + ! 2. Y(1) = (K(1)**(-1))*X(1) + ! 3. TX(1) = X(1) - A(1)*Y(1) + ! 4. DO ILEV=2, NLEV + ! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1) + ! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV) + ! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV)) + ! 5. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1) + ! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV)) + ! 6. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + ! + ! + n_col = desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + + allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), & + & mlprec_wrk(1)%ty(nr2l), mlprec_wrk(1)%tx(nr2l), stat=info) + + mlprec_wrk(1)%x2l(:) = dzero + mlprec_wrk(1)%y2l(:) = dzero + mlprec_wrk(1)%tx(:) = dzero + mlprec_wrk(1)%ty(:) = dzero + + if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') - goto 9999 + goto 9999 end if - t2l(:) = dzero - w2l(:) = dzero - tx(:) = dzero - ty(:) = dzero - tty(:) = dzero + call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,& + & baseprecv(1)%base_desc,info) + call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,& + & baseprecv(1)%base_desc,info) - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - call psb_geaxpby(done,y,dzero,ty,desc_data,info) - if(info /=0) goto 9999 + call psb_baseprc_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,& + & dzero,mlprec_wrk(1)%y2l,& + & baseprecv(1)%base_desc,& + & trans,work,info) - call psb_baseprc_aply(baseprecv(1),tx,dzero,tty,desc_data,trans,work,info) if(info /=0) goto 9999 - call psb_spmm(-done,baseprecv(2)%aorig,tty,done,tx,desc_data,info,work=work) + mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l + + call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,& + & done,mlprec_wrk(1)%ty,baseprecv(1)%base_desc,info,work=work) if(info /=0) goto 9999 - if (ismth /= no_smth_) then - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,baseprecv(1)%desc_data,info,work=work) + do ilev = 2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + ismth=baseprecv(ilev)%iprcparm(smth_kind_) + allocate(mlprec_wrk(ilev)%ty(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%x2l(nr2l), stat=info) + + mlprec_wrk(ilev)%x2l(:) = dzero + mlprec_wrk(ilev)%y2l(:) = dzero + mlprec_wrk(ilev)%tx(:) = dzero + mlprec_wrk(ilev)%ty(:) = dzero + + + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + + if (ismth /= no_smth_) then + ! + !Smoothed Aggregation + ! + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + + call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,& + & info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%ty(n_row+1:max(n_row,n_col)) = dzero + end if + + call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%ty,dzero,& + & mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero + ! + ! Raw aggregation, may take shortcuts + ! + mlprec_wrk(ilev)%x2l = dzero + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%ty(i) + end do end if - call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info) + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif + + call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,& + & baseprecv(ilev)%base_desc,info) if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcuts - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if + call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,& + & dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info) - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif + if(info /=0) goto 9999 + if(ilev < nlev) then + mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l + call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & done,mlprec_wrk(ilev)%ty,baseprecv(ilev)%base_desc,info,work=work) + if(info /=0) goto 9999 + endif - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - if(info /=0) goto 9999 + enddo - if (ismth /= no_smth_) then - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info) - if(info /=0) goto 9999 + do ilev=nlev-1, 1, -1 - call psb_geaxpby(done,ty,done,tty,desc_data,info) - if(info /=0) goto 9999 + ismth=baseprecv(ilev+1)%iprcparm(smth_kind_) + if (ismth /= no_smth_) then + if (ismth == smth_omg_) & + & call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,& + & info,work=work) + call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,& + & done,mlprec_wrk(ilev)%y2l,info) + if(info /=0) goto 9999 - else + else + n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_) + do i=1, n_row + mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & + & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) + enddo - do i=1, desc_data%matrix_data(psb_n_row_) - tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) - enddo + end if - end if + call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) - call psb_geaxpby(done,x,dzero,tx,desc_data,info) - if(info /=0) goto 9999 + if(info /=0) goto 9999 - call psb_spmm(-done,baseprecv(2)%aorig,tty,done,tx,desc_data,info,work=work) - if(info /=0) goto 9999 - call psb_baseprc_aply(baseprecv(1),tx,done,tty,desc_data,'N',work,info) + call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,& + & done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info) + if(info /=0) goto 9999 - call psb_geaxpby(done,tty,beta,y,desc_data,info) + enddo + + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& + & baseprecv(1)%base_desc,info) - deallocate(t2l,w2l,tx,ty,tty) + if(info /=0) goto 9999 case default - write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_) + call psb_errpush(4013,name,a_err='wrong smooth_pos',& + & i_Err=(/baseprecv(2)%iprcparm(smth_pos_),0,0,0,0/)) + goto 9999 end select case default - write(0,*) me, 'Wrong mltype into PRC_APLY ',& - & baseprecv(2)%iprcparm(ml_type_) + call psb_errpush(4013,name,a_err='wrong mltype',& + & i_Err=(/baseprecv(2)%iprcparm(ml_type_),0,0,0,0/)) + goto 9999 + end select + + call mlprec_wrk_free(mlprec_wrk) + deallocate(mlprec_wrk) + call psb_erractionrestore(err_act) return @@ -526,4 +776,24 @@ subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) end if return +contains + subroutine mlprec_wrk_free(wrk) + type(psb_mlprec_wrk_type) :: wrk(:) + ! This will not be needed when we have allocatables, as + ! it is sufficient to deallocate the container, and + ! the compiler is supposed to recursively deallocate the + ! various components. + integer i + + do i=1, size(wrk) + if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx) + if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty) + if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l) + if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l) + if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l) + if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty) + end do + end subroutine mlprec_wrk_free + end subroutine psb_dmlprc_aply + diff --git a/src/prec/psb_dmlprc_bld.f90 b/src/prec/psb_dmlprc_bld.f90 index f887a7bc..339ba95a 100644 --- a/src/prec/psb_dmlprc_bld.f90 +++ b/src/prec/psb_dmlprc_bld.f90 @@ -135,8 +135,6 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) - p%aorig => a - nullify(p%d) @@ -167,7 +165,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) call psb_baseprc_bld(ac,desc_p,p,info) - if (debug) write(0,*) 'Out from basaeprcbld',info + if (debug) write(0,*) 'Out from baseprcbld',info if(info /= 0) then info=4010 ch_err='psb_baseprc_bld' @@ -180,13 +178,13 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) ! We have used a separate ac because: ! 1. We want to reuse the same routines psb_ilu_bld etc. ! 2. We do NOT want to pass an argument twice to them - ! p%av(ac_) and p - ! Hence a separate AC and a TRANSFER function. + ! p%av(ac_) and p, as this would violate the Fortran standard + ! Hence a separate AC and a TRANSFER function at the end. ! call psb_sp_transfer(ac,p%av(ac_),info) - - call psb_cdfree(desc_p,info) - deallocate(desc_p) + p%base_a => p%av(ac_) + p%desc_ac => desc_p + nullify(desc_p) call psb_erractionrestore(err_act) return diff --git a/src/prec/psb_dprc_aply.f90 b/src/prec/psb_dprc_aply.f90 index 12751344..38a36a1d 100644 --- a/src/prec/psb_dprc_aply.f90 +++ b/src/prec/psb_dprc_aply.f90 @@ -61,13 +61,13 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) character(len=20) :: name interface psb_baseprc_aply - subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) + subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_descriptor_type use psb_prec_type type(psb_desc_type),intent(in) :: desc_data type(psb_dbaseprc_type), intent(in) :: prec real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -75,12 +75,12 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) end interface interface psb_mlprc_aply - subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) + subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) use psb_descriptor_type use psb_prec_type type(psb_desc_type),intent(in) :: desc_data type(psb_dbaseprc_type), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: beta + real(kind(0.d0)),intent(in) :: alpha,beta real(kind(0.d0)),intent(inout) :: x(:), y(:) character :: trans real(kind(0.d0)),target :: work(:) @@ -117,14 +117,14 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) end if if (size(prec%baseprecv) >1) then if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y) - call psb_mlprc_aply(prec%baseprecv,x,dzero,y,desc_data,trans_,work_,info) + call psb_mlprc_aply(done,prec%baseprecv,x,dzero,y,desc_data,trans_,work_,info) if(info /= 0) then call psb_errpush(4010,name,a_err='psb_dmlprc_aply') goto 9999 end if else if (size(prec%baseprecv) == 1) then - call psb_baseprc_aply(prec%baseprecv(1),x,dzero,y,desc_data,trans_, work_,info) + call psb_baseprc_aply(done,prec%baseprecv(1),x,dzero,y,desc_data,trans_, work_,info) else write(0,*) 'Inconsistent preconditioner: size of baseprecv???' endif diff --git a/src/prec/psb_dprecbld.f90 b/src/prec/psb_dprecbld.f90 index 37f7dbbd..548237ec 100644 --- a/src/prec/psb_dprecbld.f90 +++ b/src/prec/psb_dprecbld.f90 @@ -48,19 +48,19 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) use psb_penv_mod Implicit None - type(psb_dspmat_type), target :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dprec_type),intent(inout) :: p - integer, intent(out) :: info - character, intent(in), optional :: upd + type(psb_dspmat_type), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(psb_dprec_type),intent(inout) :: p + integer, intent(out) :: info + character, intent(in), optional :: upd interface psb_baseprc_bld subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) Use psb_spmat_type use psb_descriptor_type use psb_prec_type - type(psb_dspmat_type) :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), target :: a + type(psb_desc_type), intent(in), target :: desc_a type(psb_dbaseprc_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd @@ -87,7 +87,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) integer :: int_err(5) character :: iupd - logical, parameter :: debug=.false. + logical, parameter :: debug=.false. integer,parameter :: iroot=0,iout=60,ilout=40 character(len=20) :: name, ch_err @@ -149,29 +149,8 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) endif if (size(p%baseprecv) > 1) then - call init_baseprc_av(p%baseprecv(2),info) - if (info /= 0) then - info=4010 - ch_err='allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - call psb_mlprc_bld(a,desc_a,p%baseprecv(2),info) - - if(info /= 0) then - info=4010 - ch_err='psb_mlprc_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ! - ! Note: this here has not been tried out. We probably need - ! a different baseprc field %desc_ac, in case we try RAS on lev. 2 of - ! a 3-level prec. - ! - do i=3, size(p%baseprecv) + do i=2, size(p%baseprecv) call init_baseprc_av(p%baseprecv(i),info) if (info /= 0) then info=4010 @@ -180,8 +159,17 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) goto 9999 endif - call psb_mlprc_bld(p%baseprecv(i-1)%av(ac_),p%baseprecv(i-1)%desc_data,& + call psb_mlprc_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,& & p%baseprecv(i),info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name) + goto 9999 + endif + if (debug) then + write(0,*) 'Return from ',i-1,' call to mlprcbld ',info + endif + end do endif @@ -206,7 +194,7 @@ contains ! Have not decided what to do yet end if allocate(p%av(max_avsz),stat=info) - if (info /= 0) return +!!$ if (info /= 0) return do k=1,size(p%av) call psb_nullify_sp(p%av(k)) end do diff --git a/src/prec/psb_dprecset.f90 b/src/prec/psb_dprecset.f90 index 3085da69..19504a34 100644 --- a/src/prec/psb_dprecset.f90 +++ b/src/prec/psb_dprecset.f90 @@ -34,33 +34,57 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dprecset(p,ptype,info,iv,rs,rv) +subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev) use psb_serial_mod use psb_descriptor_type use psb_prec_type implicit none - type(psb_dprec_type), intent(inout) :: p character(len=*), intent(in) :: ptype integer, intent(out) :: info integer, optional, intent(in) :: iv(:) + integer, optional, intent(in) :: nlev,ilev real(kind(1.d0)), optional, intent(in) :: rs real(kind(1.d0)), optional, intent(in) :: rv(:) type(psb_dbaseprc_type), pointer :: bpv(:)=>null() character(len=len(ptype)) :: typeup - integer :: isz, err + integer :: isz, err, nlev_, ilev_, i info = 0 + if (present(ilev)) then + ilev_ = max(1, ilev) + else + ilev_ = 1 + end if + if (present(nlev)) then + if (associated(p%baseprecv)) then + write(0,*) 'Warning: NLEV is ignored when P is already associated' + end if + nlev_ = max(1, nlev) + else + nlev_ = 1 + end if if (.not.associated(p%baseprecv)) then - allocate(p%baseprecv(1),stat=err) - call psb_nullify_baseprec(p%baseprecv(1)) + allocate(p%baseprecv(nlev_),stat=err) + do i=1, nlev_ + call psb_nullify_baseprec(p%baseprecv(i)) + end do + else + nlev_ = size(p%baseprecv) endif - if (.not.associated(p%baseprecv(1)%iprcparm)) then - allocate(p%baseprecv(1)%iprcparm(ifpsz),stat=err) + if ((ilev_<1).or.(ilev_ > nlev_)) then + write(0,*) 'PRECSET ERRROR: ilev out of bounds' + info = -1 + return + endif + + if (.not.associated(p%baseprecv(ilev_)%iprcparm)) then + allocate(p%baseprecv(ilev_)%iprcparm(ifpsz),& + & p%baseprecv(ilev_)%dprcparm(dfpsz),stat=err) if (err/=0) then write(0,*)'Precset Memory Failure',err endif @@ -68,124 +92,103 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv) select case(toupper(ptype(1:len_trim(ptype)))) case ('NONE','NOPREC') - p%baseprecv(1)%iprcparm(:) = 0 - p%baseprecv(1)%iprcparm(p_type_) = noprec_ - p%baseprecv(1)%iprcparm(f_type_) = f_none_ - p%baseprecv(1)%iprcparm(restr_) = psb_none_ - p%baseprecv(1)%iprcparm(prol_) = psb_none_ - p%baseprecv(1)%iprcparm(iren_) = 0 - p%baseprecv(1)%iprcparm(n_ovr_) = 0 - p%baseprecv(1)%iprcparm(jac_sweeps_) = 1 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_ + p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_ + p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(iren_) = 0 + p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0 + p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 case ('DIAG','DIAGSC') - p%baseprecv(1)%iprcparm(:) = 0 - p%baseprecv(1)%iprcparm(p_type_) = diagsc_ - p%baseprecv(1)%iprcparm(f_type_) = f_none_ - p%baseprecv(1)%iprcparm(restr_) = psb_none_ - p%baseprecv(1)%iprcparm(prol_) = psb_none_ - p%baseprecv(1)%iprcparm(iren_) = 0 - p%baseprecv(1)%iprcparm(n_ovr_) = 0 - p%baseprecv(1)%iprcparm(jac_sweeps_) = 1 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(p_type_) = diagsc_ + p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_ + p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(iren_) = 0 + p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0 + p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 case ('BJA','ILU') - p%baseprecv(1)%iprcparm(:) = 0 - p%baseprecv(1)%iprcparm(p_type_) = bja_ - p%baseprecv(1)%iprcparm(f_type_) = f_ilu_n_ - p%baseprecv(1)%iprcparm(restr_) = psb_none_ - p%baseprecv(1)%iprcparm(prol_) = psb_none_ - p%baseprecv(1)%iprcparm(iren_) = 0 - p%baseprecv(1)%iprcparm(n_ovr_) = 0 - p%baseprecv(1)%iprcparm(ilu_fill_in_) = 0 - p%baseprecv(1)%iprcparm(jac_sweeps_) = 1 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(p_type_) = bja_ + p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_ + p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(iren_) = 0 + p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0 + p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0 + p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 case ('ASM','AS') - p%baseprecv(1)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 ! Defaults first - p%baseprecv(1)%iprcparm(p_type_) = asm_ - p%baseprecv(1)%iprcparm(f_type_) = f_ilu_n_ - p%baseprecv(1)%iprcparm(restr_) = psb_halo_ - p%baseprecv(1)%iprcparm(prol_) = psb_none_ - p%baseprecv(1)%iprcparm(iren_) = 0 - p%baseprecv(1)%iprcparm(n_ovr_) = 1 - p%baseprecv(1)%iprcparm(ilu_fill_in_) = 0 - p%baseprecv(1)%iprcparm(jac_sweeps_) = 1 + p%baseprecv(ilev_)%iprcparm(p_type_) = asm_ + p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_ + p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_ + p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(iren_) = 0 + p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1 + p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0 + p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 if (present(iv)) then isz = size(iv) - if (isz >= 1) p%baseprecv(1)%iprcparm(n_ovr_) = iv(1) - if (isz >= 2) p%baseprecv(1)%iprcparm(restr_) = iv(2) - if (isz >= 3) p%baseprecv(1)%iprcparm(prol_) = iv(3) - if (isz >= 4) p%baseprecv(1)%iprcparm(f_type_) = iv(4) + if (isz >= 1) p%baseprecv(ilev_)%iprcparm(n_ovr_) = iv(1) + if (isz >= 2) p%baseprecv(ilev_)%iprcparm(restr_) = iv(2) + if (isz >= 3) p%baseprecv(ilev_)%iprcparm(prol_) = iv(3) + if (isz >= 4) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(4) ! Do not consider renum for the time being. -!!$ if (isz >= 5) p%baseprecv(1)%iprcparm(iren_) = iv(5) +!!$ if (isz >= 5) p%baseprecv(ilev_)%iprcparm(iren_) = iv(5) end if case ('ML', '2L', '2LEV') - select case (size(p%baseprecv)) - case(1) - ! Reallocate - allocate(bpv(2),stat=err) - if (err/=0) then - write(0,*)'Precset Memory Failure 2l:1',err - endif - bpv(1) = p%baseprecv(1) - call psb_nullify_baseprec(bpv(2)) - deallocate(p%baseprecv) - p%baseprecv => bpv - nullify(bpv) - - case(2) - ! Do nothing - - case default - ! Error - - end select - - allocate(p%baseprecv(2)%iprcparm(ifpsz),stat=err) - if (err/=0) then - write(0,*)'Precset Memory Failure 2l:2',err - endif - allocate(p%baseprecv(2)%dprcparm(dfpsz),stat=err) - if (err/=0) then - write(0,*)'Precset Memory Failure 2l:3',err - endif - - p%baseprecv(2)%iprcparm(:) = 0 - - p%baseprecv(2)%iprcparm(p_type_) = bja_ - p%baseprecv(2)%iprcparm(restr_) = psb_none_ - p%baseprecv(2)%iprcparm(prol_) = psb_none_ - p%baseprecv(2)%iprcparm(iren_) = 0 - p%baseprecv(2)%iprcparm(n_ovr_) = 0 - p%baseprecv(2)%iprcparm(ml_type_) = mult_ml_prec_ - p%baseprecv(2)%iprcparm(aggr_alg_) = loc_aggr_ - p%baseprecv(2)%iprcparm(smth_kind_) = smth_omg_ - p%baseprecv(2)%iprcparm(coarse_mat_) = mat_distr_ - p%baseprecv(2)%iprcparm(smth_pos_) = post_smooth_ - p%baseprecv(2)%iprcparm(glb_smth_) = 1 - p%baseprecv(2)%iprcparm(om_choice_) = lib_choice_ - p%baseprecv(2)%iprcparm(f_type_) = f_ilu_n_ - p%baseprecv(2)%iprcparm(ilu_fill_in_) = 0 - p%baseprecv(2)%dprcparm(smooth_omega_) = 4.d0/3.d0 - p%baseprecv(2)%iprcparm(jac_sweeps_) = 1 +!!$ allocate(p%baseprecv(ilev_)%iprcparm(ifpsz),stat=err) +!!$ if (err/=0) then +!!$ write(0,*)'Precset Memory Failure 2l:2',err +!!$ endif +!!$ allocate(p%baseprecv(ilev_)%dprcparm(dfpsz),stat=err) +!!$ if (err/=0) then +!!$ write(0,*)'Precset Memory Failure 2l:3',err +!!$ endif + + p%baseprecv(ilev_)%iprcparm(:) = 0 + + p%baseprecv(ilev_)%iprcparm(p_type_) = bja_ + p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_ + p%baseprecv(ilev_)%iprcparm(iren_) = 0 + p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0 + p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_ + p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_ + p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_ + p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_ + p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_ + p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1 + p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_ + p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_ + p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0 + p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0 + p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1 if (present(iv)) then isz = size(iv) - if (isz >= 1) p%baseprecv(2)%iprcparm(ml_type_) = iv(1) - if (isz >= 2) p%baseprecv(2)%iprcparm(aggr_alg_) = iv(2) - if (isz >= 3) p%baseprecv(2)%iprcparm(smth_kind_) = iv(3) - if (isz >= 4) p%baseprecv(2)%iprcparm(coarse_mat_) = iv(4) - if (isz >= 5) p%baseprecv(2)%iprcparm(smth_pos_) = iv(5) - if (isz >= 6) p%baseprecv(2)%iprcparm(f_type_) = iv(6) - if (isz >= 7) p%baseprecv(2)%iprcparm(jac_sweeps_) = iv(7) + if (isz >= 1) p%baseprecv(ilev_)%iprcparm(ml_type_) = iv(1) + if (isz >= 2) p%baseprecv(ilev_)%iprcparm(aggr_alg_) = iv(2) + if (isz >= 3) p%baseprecv(ilev_)%iprcparm(smth_kind_) = iv(3) + if (isz >= 4) p%baseprecv(ilev_)%iprcparm(coarse_mat_) = iv(4) + if (isz >= 5) p%baseprecv(ilev_)%iprcparm(smth_pos_) = iv(5) + if (isz >= 6) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(6) + if (isz >= 7) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(7) end if if (present(rs)) then - p%baseprecv(2)%iprcparm(om_choice_) = user_choice_ - p%baseprecv(2)%dprcparm(smooth_omega_) = rs + p%baseprecv(ilev_)%iprcparm(om_choice_) = user_choice_ + p%baseprecv(ilev_)%dprcparm(smooth_omega_) = rs end if diff --git a/src/prec/psb_zbaseprc_aply.f90 b/src/prec/psb_zbaseprc_aply.f90 index 15d79655..35f117c6 100644 --- a/src/prec/psb_zbaseprc_aply.f90 +++ b/src/prec/psb_zbaseprc_aply.f90 @@ -34,9 +34,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) +subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! - ! Compute Y <- beta*Y + K^-1 X + ! Compute Y <- beta*Y + alpha*K^-1 X ! where K is a a basic preconditioner stored in prec ! @@ -52,7 +52,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) type(psb_desc_type),intent(in) :: desc_data type(psb_zbaseprc_type), intent(in) :: prec complex(kind(0.d0)),intent(inout) :: x(:), y(:) - complex(kind(0.d0)),intent(in) :: beta + complex(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans complex(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -68,13 +68,13 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) character(len=20) :: name, ch_err interface psb_bjac_aply - subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) + subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_descriptor_type use psb_prec_type type(psb_desc_type), intent(in) :: desc_data type(psb_zbaseprc_type), intent(in) :: prec complex(kind(0.d0)),intent(inout) :: x(:), y(:) - complex(kind(0.d0)),intent(in) :: beta + complex(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans complex(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -105,33 +105,35 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) case(noprec_) - n_row=desc_data%matrix_data(psb_n_row_) - if (beta==zzero) then - y(1:n_row) = x(1:n_row) - else if (beta==zone) then - y(1:n_row) = x(1:n_row) + y(1:n_row) - else if (beta==-zone) then - y(1:n_row) = x(1:n_row) - y(1:n_row) - else - y(1:n_row) = x(1:n_row) + beta*y(1:n_row) - end if + call psb_geaxpby(alpha,x,beta,y,desc_data,info) case(diagsc_) + + if (size(work) >= size(x)) then + ww => work + else + allocate(ww(size(x)),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + end if n_row=desc_data%matrix_data(psb_n_row_) - if (beta==zzero) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - else if (beta==zone) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row) - else if (beta==-zone) then - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row) - else - y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row) + ww(1:n_row) = x(1:n_row)*prec%d(1:n_row) + 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 case(bja_) - call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) + call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 ch_err='psb_bjac_aply' @@ -142,7 +144,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) if (prec%iprcparm(n_ovr_)==0) then ! shortcut: this fixes performance for RAS(0) == BJA - call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) + call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 ch_err='psb_bjacaply' @@ -214,7 +216,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) end if endif - call psb_bjac_aply(prec,tx,zzero,ty,prec%desc_data,trans,aux,info) + call psb_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,trans,aux,info) if(info.ne.0) then info=4010 ch_err='psb_bjac_aply' @@ -250,18 +252,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) & prec%iprcparm(prol_) end select - if (beta == zzero) then - y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_)) - else if (beta == zone) then - y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - else if (beta == -zone) then - y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - else - y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + & - & ty(1:desc_data%matrix_data(psb_n_row_)) - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) if ((6*isz) <= size(work)) then diff --git a/src/prec/psb_zbaseprc_bld.f90 b/src/prec/psb_zbaseprc_bld.f90 index 0ef2bb3b..bc3c2874 100644 --- a/src/prec/psb_zbaseprc_bld.f90 +++ b/src/prec/psb_zbaseprc_bld.f90 @@ -49,7 +49,7 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) Implicit None type(psb_zspmat_type), target :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(in), target :: desc_a type(psb_zbaseprc_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd @@ -88,7 +88,7 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) use psb_const_mod implicit none - type(psb_zspmat_type), intent(in) :: a + type(psb_zspmat_type), intent(inout) :: a type(psb_desc_type), intent(in) :: desc_a type(psb_zbaseprc_type), intent(inout) :: p integer, intent(out) :: info @@ -169,7 +169,13 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) select case(p%iprcparm(p_type_)) case (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 (diagsc_) @@ -256,7 +262,8 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) end select - + p%base_a => a + p%base_desc => desc_a call psb_erractionrestore(err_act) return diff --git a/src/prec/psb_zbjac_aply.f90 b/src/prec/psb_zbjac_aply.f90 index ef8e2cbe..9442a7a2 100644 --- a/src/prec/psb_zbjac_aply.f90 +++ b/src/prec/psb_zbjac_aply.f90 @@ -34,9 +34,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) +subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! - ! Compute Y <- beta*Y + K^-1 X + ! Compute Y <- beta*Y + alpha*K^-1 X ! where K is a a Block Jacobi preconditioner stored in prec ! Note that desc_data may or may not be the same as prec%desc_data, ! but since both are INTENT(IN) this should be legal. @@ -54,7 +54,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) type(psb_desc_type), intent(in) :: desc_data type(psb_zbaseprc_type), intent(in) :: prec complex(kind(0.d0)),intent(inout) :: x(:), y(:) - complex(kind(0.d0)),intent(in) :: beta + complex(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans complex(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -125,7 +125,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) & trans='N',unit=diagl,choice=psb_none_,work=aux) if(info /=0) goto 9999 ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(zone,prec%av(u_pr_),ww,beta,y,desc_data,info,& + call psb_spsm(alpha,prec%av(u_pr_),ww,beta,y,desc_data,info,& & trans='N',unit=diagu,choice=psb_none_, work=aux) if(info /=0) goto 9999 @@ -134,7 +134,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) & trans=trans,unit=diagu,choice=psb_none_, work=aux) if(info /=0) goto 9999 ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row) - call psb_spsm(zone,prec%av(l_pr_),ww,beta,y,desc_data,info,& + call psb_spsm(alpha,prec%av(l_pr_),ww,beta,y,desc_data,info,& & trans=trans,unit=diagl,choice=psb_none_,work=aux) if(info /=0) goto 9999 @@ -152,16 +152,8 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) end select if(info /=0) goto 9999 + call psb_geaxpby(alpha,ww,beta,y,desc_data,info) - if (beta == zzero) then - y(1:n_row) = ww(1:n_row) - else if (beta==zone) then - y(1:n_row) = ww(1:n_row) + y(1:n_row) - else if (beta==-zone) then - y(1:n_row) = ww(1:n_row) - y(1:n_row) - else - y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) - endif case (f_umf_) @@ -174,15 +166,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) if(info /=0) goto 9999 - if (beta == zzero) then - y(1:n_row) = ww(1:n_row) - else if (beta==zone) then - y(1:n_row) = ww(1:n_row) + y(1:n_row) - else if (beta==-zone) then - y(1:n_row) = ww(1:n_row) - y(1:n_row) - else - y(1:n_row) = ww(1:n_row) + beta*y(1:n_row) - endif + call psb_geaxpby(alpha,ww,beta,y,desc_data,info) case default write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_) @@ -253,15 +237,8 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info) end select - if (beta == zzero) then - y(1:n_row) = tx(1:n_row) - else if (beta==zone) then - y(1:n_row) = tx(1:n_row) + y(1:n_row) - else if (beta==-zone) then - y(1:n_row) = tx(1:n_row) - y(1:n_row) - else - y(1:n_row) = tx(1:n_row) + beta*y(1:n_row) - endif + call psb_geaxpby(alpha,tx,beta,y,desc_data,info) + deallocate(tx,ty) diff --git a/src/prec/psb_zbldaggrmat.f90 b/src/prec/psb_zbldaggrmat.f90 index d7d1defe..fb57d93a 100644 --- a/src/prec/psb_zbldaggrmat.f90 +++ b/src/prec/psb_zbldaggrmat.f90 @@ -620,7 +620,12 @@ contains ! Doing it this way means to consider diag(Ai) ! ! - call psb_symbmm(am3,am4,am1) + call psb_symbmm(am3,am4,am1,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm') + goto 9999 + end if + call psb_numbmm(am3,am4,am1) if (debug) write(0,*) me,'Done NUMBMM 1' @@ -667,7 +672,12 @@ contains if (test_dump) & & call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob) - call psb_symbmm(a,am1,am3) + call psb_symbmm(a,am1,am3,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm') + goto 9999 + end if + call psb_numbmm(a,am1,am3) if (debug) write(0,*) me,'Done NUMBMM 2' @@ -724,7 +734,12 @@ contains endif if (debug) write(0,*) me,'starting symbmm 3' - call psb_symbmm(am2,am3,b) + call psb_symbmm(am2,am3,b,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='symbmm') + goto 9999 + end if + if (debug) write(0,*) me,'starting numbmm 3' call psb_numbmm(am2,am3,b) if (debug) write(0,*) me,'Done NUMBMM 3' diff --git a/src/prec/psb_zmlprc_aply.f90 b/src/prec/psb_zmlprc_aply.f90 index 6ef7fbee..216fd720 100644 --- a/src/prec/psb_zmlprc_aply.f90 +++ b/src/prec/psb_zmlprc_aply.f90 @@ -34,11 +34,54 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) +subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) ! - ! Compute Y <- beta*Y + K^-1 X - ! where K is a multilevel (actually 2-level) preconditioner stored in prec + ! Compute Y <- beta*Y + alpha*K^-1 X + ! where K is a multilevel preconditioner stored in baseprecv ! + ! cfr.: Smith, Biorstad & Gropp + ! Domain Decomposition + ! Cambridge Univ. Press + ! + ! To each level I there corresponds a matrix A(I) and a preconditioner K(I) + ! + ! A notational difference: in the DD reference above the preconditioner for + ! a given level K(I) is written out as a sum over the subdomains + ! + ! SUM_k(R_k^T A_k R_k) + ! + ! whereas in this code the sum is implicit in the parallelization, + ! i.e. each process takes care of one subdomain, and for each level we have + ! as many subdomains as there are processes (except for the coarsest level where + ! we might have a replicated index space). Thus the sum apparently disappears + ! from our code, but only apparently, because it is implicit in the call + ! to psb_baseprc_aply. + ! + ! A bit of description of the baseprecv(:) data structure: + ! 1. Number of levels = NLEV = size(baseprecv(:)) + ! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level. + ! Includes: + ! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners + ! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners + ! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps + ! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV + ! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors + ! (ilev-1) ---> (ilev) + ! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors + ! (ilev) ---> (ilev-1) + ! Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe + ! + ! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV + ! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix + ! of the current level, i.e.: if ILEV=1 then A + ! else the aggregated matrix av(ac_); so we have + ! a unified treatment of residuals. Need this to + ! avoid passing explicitly matrix A to the + ! outer prec. routine + ! 5. baseprecv(ilev)%mlia The aggregation map from (ilev-1)-->(ilev) + ! if no smoother, it is used instead of sm_pr_ + ! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs. + ! use psb_serial_mod use psb_descriptor_type @@ -52,7 +95,7 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) type(psb_desc_type),intent(in) :: desc_data type(psb_zbaseprc_type), intent(in) :: baseprecv(:) - complex(kind(1.d0)),intent(in) :: beta + complex(kind(1.d0)),intent(in) :: alpha,beta complex(kind(1.d0)),intent(inout) :: x(:), y(:) character :: trans complex(kind(1.d0)),target :: work(:) @@ -68,25 +111,32 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) real(kind(1.d0)) :: omega real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime logical, parameter :: debug=.false., debugprt=.false. - integer :: ismth + integer :: ismth, nlev, ilev external mpi_wtime character(len=20) :: name, ch_err + type psb_mlprec_wrk_type + complex(kind(1.d0)), pointer :: tx(:)=>null(),ty(:)=>null(),& + & x2l(:)=>null(),y2l(:)=>null(),& + & b2l(:)=>null(),tty(:)=>null() + end type psb_mlprec_wrk_type + type(psb_mlprec_wrk_type), pointer :: mlprec_wrk(:) + interface psb_baseprc_aply - subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) - use psb_descriptor_type - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_zbaseprc_type), intent(in) :: prec - complex(kind(1.d0)),intent(inout) :: x(:), y(:) - complex(kind(1.d0)),intent(in) :: beta - character(len=1) :: trans - complex(kind(1.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_zbaseprc_aply + subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + use psb_descriptor_type + use psb_prec_type + type(psb_desc_type),intent(in) :: desc_data + type(psb_zbaseprc_type), intent(in) :: prec + complex(kind(1.d0)),intent(inout) :: x(:), y(:) + complex(kind(1.d0)),intent(in) :: alpha,beta + character(len=1) :: trans + complex(kind(1.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_zbaseprc_aply end interface - name='psb_zmlprc_aply' + name='psb_mlprc_aply' info = 0 call psb_erractionsave(err_act) @@ -94,426 +144,625 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) ictxt=desc_data%matrix_data(psb_ctxt_) call psb_info(ictxt, me, np) - omega=baseprecv(2)%dprcparm(smooth_omega_) - ismth=baseprecv(2)%iprcparm(smth_kind_) + nlev = size(baseprecv) + allocate(mlprec_wrk(nlev),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + select case(baseprecv(2)%iprcparm(ml_type_)) + case(no_ml_) ! Should not really get here. - write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????' + call psb_errpush(4010,name,a_err='no_ml_ in mlprc_aply?') + goto 9999 + case(add_ml_prec_) - ! - ! Additive multilevel + ! - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - call psb_baseprc_aply(baseprecv(1),x,beta,y,desc_data,trans,work,info) - if(info /=0) goto 9999 - - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - t2l(:) = zzero - w2l(:) = zzero + ! Additive is very simple. + ! 1. X(1) = Xext + ! 2. DO ILEV=2,NLEV + ! X(ILEV) = AV(PR_SM_T_)*X(ILEV-1) + ! 3. Y(ILEV) = (K(ILEV)**(-1))*X(ILEV) + ! 4. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = AV(PR_SM_)*Y(ILEV+1) + ! 5. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + - if (ismth /= no_smth_) then - ! - ! Smoothed aggregation - ! - allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),& - & tz(max(n_row,n_col)),stat=info) + call psb_baseprc_aply(alpha,baseprecv(1),x,beta,y,& + & baseprecv(1)%base_desc,trans,work,info) + if(info /=0) goto 9999 + allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y))) + mlprec_wrk(1)%x2l(:) = x(:) + + + do ilev = 2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + allocate(mlprec_wrk(ilev)%x2l(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%tx(max(n_row,n_col)),& + & mlprec_wrk(ilev)%ty(max(n_row,n_col)), stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_)) - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero - ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero - tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero + mlprec_wrk(ilev)%x2l(:) = zzero + mlprec_wrk(ilev)%y2l(:) = zzero + mlprec_wrk(ilev)%tx(1:n_row) = mlprec_wrk(ilev-1)%x2l(1:n_row) + mlprec_wrk(ilev)%tx(n_row+1:max(n_row,n_col)) = zzero + mlprec_wrk(ilev)%ty(:) = zzero + + ismth=baseprecv(ilev)%iprcparm(smth_kind_) + + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,& + & info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = zzero + end if - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) + call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l,& + & zzero,mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero - end if + ! + ! Raw aggregation, may take shortcut + ! + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%x2l(i) + end do - call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info) - if(info /=0) goto 9999 + end if - else - ! - ! Raw aggregation, may take shortcut - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i) - end do + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif - end if + call psb_baseprc_aply(zone,baseprecv(ilev),& + & mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,& + & baseprecv(ilev)%desc_data, 'N',work,info) - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif + enddo - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) + do ilev =nlev,2,-1 + ismth=baseprecv(ilev)%iprcparm(smth_kind_) + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) - if (ismth /= no_smth_) then + if (ismth /= no_smth_) then - call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info) - if(info /=0) goto 9999 - ! - ! Finally add back into Y. - ! - call psb_geaxpby(zone,ty,zone,y,desc_data,info) - if(info /=0) goto 9999 - deallocate(tx,ty,tz) + call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_),mlprec_wrk(ilev)%y2l,& + & zone,mlprec_wrk(ilev-1)%y2l,info) + if(info /=0) goto 9999 - else + else - do i=1, desc_data%matrix_data(psb_n_row_) - y(i) = y(i) + t2l(baseprecv(2)%mlia(i)) - enddo + do i=1, n_row + mlprec_wrk(ilev-1)%y2l(i) = mlprec_wrk(ilev-1)%y2l(i) + & + & mlprec_wrk(ilev)%y2l(baseprecv(ilev)%mlia(i)) + enddo - end if + end if + end do - if (debugprt) write(0,*)' Y2: ',Y(:) + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,zone,y,baseprecv(1)%base_desc,info) + if(info /=0) goto 9999 - deallocate(t2l,w2l) case(mult_ml_prec_) ! ! Multiplicative multilevel ! Pre/post smoothing versions. + ! select case(baseprecv(2)%iprcparm(smth_pos_)) case(post_smooth_) + + ! + ! Post smoothing. + ! 1. X(1) = Xext + ! 2. DO ILEV=2, NLEV :: X(ILEV) = AV(PR_SM_T_,ILEV)*X(ILEV-1) + ! 3. Y(NLEV) = (K(NLEV)**(-1))*X(NLEV) + ! 4. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = AV(PR_SM_,ILEV+1)*Y(ILEV+1) + ! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV)) + ! + ! 5. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + ! + ! Also: post smoothing is not spelled out in detail in DD. + ! + ! + - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if + n_col = desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - t2l(:) = zzero - w2l(:) = zzero + allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), & + & mlprec_wrk(1)%tx(nr2l), stat=info) + mlprec_wrk(1)%x2l(:) = zzero + mlprec_wrk(1)%y2l(:) = zzero + mlprec_wrk(1)%tx(:) = zzero - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zzero - ! + call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& + & baseprecv(1)%base_desc,info) + call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& + & baseprecv(1)%base_desc,info) + + do ilev=2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + ismth = baseprecv(ilev)%iprcparm(smth_kind_) + + allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%x2l(nr2l), stat=info) - if (debug) write(0,*)' mult_ml_apply omega ',omega - if (debugprt) write(0,*)' mult_ml_apply X: ',X(:) - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - if(info /=0) then - if (debug) write(0,*)' From axpby1 ',size(x),size(tx),n_row,n_col,nr2l,nrg - call psb_errpush(4010,name,a_err='axpby post_smooth 1') - goto 9999 - endif - if (ismth /= no_smth_) then - ! - ! Smoothed aggregation - ! - allocate(tz(max(n_row,n_col)),stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) + mlprec_wrk(ilev)%x2l(:) = zzero + mlprec_wrk(ilev)%y2l(:) = zzero + mlprec_wrk(ilev)%tx(:) = zzero + if (ismth /= no_smth_) then + ! + ! Smoothed aggregation + ! + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + call psb_halo(mlprec_wrk(ilev-1)%x2l,& + & baseprecv(ilev-1)%base_desc,info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = zzero + end if + + call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l, & + & zzero,mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero + ! + ! Raw aggregation, may take shortcut + ! + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%x2l(i) + end do end if - call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info) + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif + call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& + & baseprecv(ilev)%base_desc,info) if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcut - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if + enddo - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,& - &'N',work,info) + call psb_baseprc_aply(zone,baseprecv(nlev),mlprec_wrk(nlev)%x2l, & + & zzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info) + if(info /=0) goto 9999 - if (ismth /= no_smth_) then - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info) - if(info /=0) goto 9999 - ! - ! Finally add back into Y. - ! - deallocate(tz) - else - ty(:) = zzero - do i=1, desc_data%matrix_data(psb_n_row_) - ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i)) - enddo - end if - deallocate(t2l,w2l) + do ilev=nlev-1, 1, -1 + ismth = baseprecv(ilev+1)%iprcparm(smth_kind_) + if (ismth /= no_smth_) then + if (ismth == smth_omg_) & + & call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,& + & info,work=work) + call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,& + & zzero,mlprec_wrk(ilev)%y2l,info) + if(info /=0) goto 9999 - call psb_spmm(-zone,baseprecv(2)%aorig,ty,zone,tx,desc_data,info,work=work) - if(info /=0) goto 9999 + else + n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_) + mlprec_wrk(ilev)%y2l(:) = zzero + do i=1, n_row + mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & + & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) + enddo - call psb_baseprc_aply(baseprecv(1),tx,zone,ty,desc_data,trans,& - & work,info) - if(info /=0) goto 9999 + end if - call psb_geaxpby(zone,ty,beta,y,desc_data,info) - if(info /=0) goto 9999 + call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) - deallocate(tx,ty) + if(info /=0) goto 9999 + call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,& + & zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info) + + if(info /=0) goto 9999 + + enddo + + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info) + + if(info /=0) goto 9999 case(pre_smooth_) - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) + + ! + ! Pre smoothing. + ! 1. X(1) = Xext + ! 2. Y(1) = (K(1)**(-1))*X(1) + ! 3. TX(1) = X(1) - A(1)*Y(1) + ! 4. DO ILEV=2, NLEV + ! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1) + ! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV) + ! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV)) + ! 5. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1) + ! 6. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + ! + ! + + n_col = desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + + allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), & + & mlprec_wrk(1)%tx(nr2l), stat=info) if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - t2l(:) = zzero - w2l(:) = zzero + mlprec_wrk(1)%y2l(:) = zzero - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - if(info /=0) goto 9999 - call psb_baseprc_aply(baseprecv(1),x,zzero,tty,desc_data,& + mlprec_wrk(1)%x2l(:) = x + + call psb_baseprc_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,& + & zzero,mlprec_wrk(1)%y2l,& + & baseprecv(1)%base_desc,& & trans,work,info) + if(info /=0) goto 9999 - call psb_spmm(-zone,baseprecv(2)%aorig,tty,zone,tx,desc_data,info,work=work) + mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l + + call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,& + & zone,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work) if(info /=0) goto 9999 - if (ismth /= no_smth_) then - allocate(tz(max(n_row,n_col)),stat=info) + do ilev = 2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + ismth = baseprecv(ilev)%iprcparm(smth_kind_) + allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%x2l(nr2l), stat=info) + + if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if + mlprec_wrk(ilev)%x2l(:) = zzero + mlprec_wrk(ilev)%y2l(:) = zzero + mlprec_wrk(ilev)%tx(:) = zzero + - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,desc_data,info,work=work) + if (ismth /= no_smth_) then + ! + !Smoothed Aggregation + ! + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + + call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,& + & info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%tx(n_row+1:max(n_row,n_col)) = zzero + end if + + call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%tx,zzero,& + & mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero + ! + ! Raw aggregation, may take shortcuts + ! + mlprec_wrk(ilev)%x2l = zzero + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%tx(i) + end do end if - call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info) + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif + + + call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,& + & zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info) + if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcuts - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if + if(ilev < nlev) then + mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l + call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) + if(info /=0) goto 9999 + endif - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif + enddo - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - if(info /=0) goto 9999 + do ilev = nlev-1, 1, -1 - if (ismth /= no_smth_) then + ismth=baseprecv(ilev+1)%iprcparm(smth_kind_) - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info) - if(info /=0) goto 9999 + if (ismth /= no_smth_) then - call psb_geaxpby(zone,ty,zone,tty,desc_data,info) - if(info /=0) goto 9999 + if (ismth == smth_omg_) & + & call psb_halo(mlprec_wrk(ilev+1)%y2l,& + & baseprecv(ilev+1)%desc_data,info,work=work) + call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,& + & zone,mlprec_wrk(ilev)%y2l,info) - deallocate(tz) - else + if(info /=0) goto 9999 - do i=1, desc_data%matrix_data(psb_n_row_) - tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) - enddo + else - end if + n_row = baseprecv(ilev+1)%base_desc%matrix_data(psb_n_row_) + do i=1, n_row + mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & + & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) + enddo + + end if + + enddo + + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& + & baseprecv(1)%base_desc,info) - call psb_geaxpby(zone,tty,beta,y,desc_data,info) if(info /=0) goto 9999 - deallocate(t2l,w2l,tx,ty,tty) case(smooth_both_) - t1 = mpi_wtime() - n_row = desc_data%matrix_data(psb_n_row_) - n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_) - nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_) - nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_) - allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info) + ! + ! Symmetrized smoothing. + ! 1. X(1) = Xext + ! 2. Y(1) = (K(1)**(-1))*X(1) + ! 3. TX(1) = X(1) - A(1)*Y(1) + ! 4. DO ILEV=2, NLEV + ! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1) + ! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV) + ! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV)) + ! 5. DO ILEV=NLEV-1,1,-1 + ! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1) + ! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV)) + ! 6. Yext = beta*Yext + Y(1) + ! + ! Note: level numbering reversed wrt ref. DD, i.e. + ! 1..NLEV <=> (j) <-> 0 + ! + ! + n_col = desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_) + + allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), & + & mlprec_wrk(1)%ty(nr2l), mlprec_wrk(1)%tx(nr2l), stat=info) + + mlprec_wrk(1)%x2l(:) = zzero + mlprec_wrk(1)%y2l(:) = zzero + mlprec_wrk(1)%tx(:) = zzero + mlprec_wrk(1)%ty(:) = zzero + + if (info /= 0) then call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - t2l(:) = zzero - w2l(:) = zzero - tx(:) = zzero - ty(:) = zzero - tty(:) = zzero + call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,& + & baseprecv(1)%base_desc,info) + call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,& + & baseprecv(1)%base_desc,info) - ! - ! Need temp copies to handle Y<- betaY + K^-1 X - ! One of the temp copies is not strictly needed when beta==zero - ! - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - call psb_geaxpby(zone,y,zzero,ty,desc_data,info) - if(info /=0) goto 9999 + call psb_baseprc_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,& + & zzero,mlprec_wrk(1)%y2l,& + & baseprecv(1)%base_desc,& + & trans,work,info) - call psb_baseprc_aply(baseprecv(1),tx,zzero,tty,desc_data,trans,work,info) if(info /=0) goto 9999 - call psb_spmm(-zone,baseprecv(2)%aorig,tty,zone,tx,desc_data,info,work=work) + mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l + + call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,& + & zone,mlprec_wrk(1)%ty,baseprecv(1)%base_desc,info,work=work) if(info /=0) goto 9999 - if (ismth /= no_smth_) then - if (baseprecv(2)%iprcparm(glb_smth_) >0) then - call psb_halo(tx,baseprecv(1)%desc_data,info,work=work) + do ilev = 2, nlev + n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_) + n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_) + nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_) + nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_) + ismth=baseprecv(ilev)%iprcparm(smth_kind_) + allocate(mlprec_wrk(ilev)%ty(nr2l),mlprec_wrk(ilev)%y2l(nr2l),& + & mlprec_wrk(ilev)%x2l(nr2l), stat=info) + + mlprec_wrk(ilev)%x2l(:) = zzero + mlprec_wrk(ilev)%y2l(:) = zzero + mlprec_wrk(ilev)%tx(:) = zzero + mlprec_wrk(ilev)%ty(:) = zzero + + + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + + if (ismth /= no_smth_) then + ! + !Smoothed Aggregation + ! + if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then + + call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,& + & info,work=work) + if(info /=0) goto 9999 + else + mlprec_wrk(ilev-1)%ty(n_row+1:max(n_row,n_col)) = zzero + end if + + call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%ty,zzero,& + & mlprec_wrk(ilev)%x2l,info) if(info /=0) goto 9999 + else - tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero + ! + ! Raw aggregation, may take shortcuts + ! + mlprec_wrk(ilev)%x2l = zzero + do i=1,n_row + mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = & + & mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + & + & mlprec_wrk(ilev-1)%ty(i) + end do end if - call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info) + if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then + call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg)) + else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then + write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& + & baseprecv(ilev)%iprcparm(coarse_mat_) + endif + + call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,& + & baseprecv(ilev)%base_desc,info) if(info /=0) goto 9999 - else - ! - ! Raw aggregation, may take shortcuts - ! - do i=1,desc_data%matrix_data(psb_n_row_) - t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i) - end do - end if + call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,& + & zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info) - if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then - call psb_sum(ictxt,t2l(1:nrg)) - else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then - write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',& - & baseprecv(2)%iprcparm(coarse_mat_) - endif + if(info /=0) goto 9999 + if(ilev < nlev) then + mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l + call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & zone,mlprec_wrk(ilev)%ty,baseprecv(ilev)%base_desc,info,work=work) + if(info /=0) goto 9999 + endif - t6 = mpi_wtime() - w2l=t2l - call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,& - & 'N',work,info) - if(info /=0) goto 9999 + enddo - if (ismth /= no_smth_) then - if (ismth == smth_omg_) & - & call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work) - call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info) - if(info /=0) goto 9999 + do ilev=nlev-1, 1, -1 - call psb_geaxpby(zone,ty,zone,tty,desc_data,info) - if(info /=0) goto 9999 + ismth=baseprecv(ilev+1)%iprcparm(smth_kind_) + if (ismth /= no_smth_) then + if (ismth == smth_omg_) & + & call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,& + & info,work=work) + call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,& + & zone,mlprec_wrk(ilev)%y2l,info) + if(info /=0) goto 9999 - else + else + n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_) + do i=1, n_row + mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + & + & mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i)) + enddo - do i=1, desc_data%matrix_data(psb_n_row_) - tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i)) - enddo + end if - end if - - call psb_geaxpby(zone,x,zzero,tx,desc_data,info) - if(info /=0) goto 9999 + call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,& + & zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work) - call psb_spmm(-zone,baseprecv(2)%aorig,tty,zone,tx,desc_data,info,work=work) - if(info /=0) goto 9999 - call psb_baseprc_aply(baseprecv(1),tx,zone,tty,desc_data,'N',work,info) + if(info /=0) goto 9999 + call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,& + & zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info) - call psb_geaxpby(zone,tty,beta,y,desc_data,info) - - deallocate(t2l,w2l,tx,ty,tty) + if(info /=0) goto 9999 + + enddo + + call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,& + & baseprecv(1)%base_desc,info) + + if(info /=0) goto 9999 case default - write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_) + call psb_errpush(4013,name,a_err='wrong smooth_pos',& + & i_Err=(/baseprecv(2)%iprcparm(smth_pos_),0,0,0,0/)) + goto 9999 end select case default - write(0,*) me, 'Wrong mltype into PRC_APLY ',& - & baseprecv(2)%iprcparm(ml_type_) + call psb_errpush(4013,name,a_err='wrong mltype',& + & i_Err=(/baseprecv(2)%iprcparm(ml_type_),0,0,0,0/)) + goto 9999 + end select + + call mlprec_wrk_free(mlprec_wrk) + deallocate(mlprec_wrk) + call psb_erractionrestore(err_act) return @@ -526,4 +775,24 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) end if return +contains + subroutine mlprec_wrk_free(wrk) + type(psb_mlprec_wrk_type) :: wrk(:) + ! This will not be needed when we have allocatables, as + ! it is sufficient to deallocate the container, and + ! the compiler is supposed to recursively deallocate the + ! various components. + integer i + + do i=1, size(wrk) + if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx) + if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty) + if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l) + if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l) + if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l) + if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty) + end do + end subroutine mlprec_wrk_free + end subroutine psb_zmlprc_aply + diff --git a/src/prec/psb_zmlprc_bld.f90 b/src/prec/psb_zmlprc_bld.f90 index 2ae4cf73..aa496839 100644 --- a/src/prec/psb_zmlprc_bld.f90 +++ b/src/prec/psb_zmlprc_bld.f90 @@ -135,8 +135,6 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info) - p%aorig => a - nullify(p%d) @@ -180,8 +178,8 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info) ! We have used a separate ac because: ! 1. We want to reuse the same routines psb_ilu_bld etc. ! 2. We do NOT want to pass an argument twice to them - ! p%av(ac_) and p - ! Hence a separate AC and a TRANSFER function. + ! p%av(ac_) and p, as this would violate the Fortran standard + ! Hence a separate AC and a TRANSFER function at the end. ! call psb_sp_transfer(ac,p%av(ac_),info) diff --git a/src/prec/psb_zprc_aply.f90 b/src/prec/psb_zprc_aply.f90 index 261168a1..5fd87cc8 100644 --- a/src/prec/psb_zprc_aply.f90 +++ b/src/prec/psb_zprc_aply.f90 @@ -61,13 +61,13 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) character(len=20) :: name interface psb_baseprc_aply - subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) + subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_descriptor_type use psb_prec_type type(psb_desc_type),intent(in) :: desc_data type(psb_zbaseprc_type), intent(in) :: prec complex(kind(0.d0)),intent(inout) :: x(:), y(:) - complex(kind(0.d0)),intent(in) :: beta + complex(kind(0.d0)),intent(in) :: alpha,beta character(len=1) :: trans complex(kind(0.d0)),target :: work(:) integer, intent(out) :: info @@ -75,12 +75,12 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) end interface interface psb_mlprc_aply - subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info) + subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) use psb_descriptor_type use psb_prec_type type(psb_desc_type),intent(in) :: desc_data type(psb_zbaseprc_type), intent(in) :: baseprecv(:) - complex(kind(0.d0)),intent(in) :: beta + complex(kind(0.d0)),intent(in) :: alpha,beta complex(kind(0.d0)),intent(inout) :: x(:), y(:) character :: trans complex(kind(0.d0)),target :: work(:) @@ -117,14 +117,14 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) end if if (size(prec%baseprecv) >1) then if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y) - call psb_mlprc_aply(prec%baseprecv,x,zzero,y,desc_data,trans_,work_,info) + call psb_mlprc_aply(zone,prec%baseprecv,x,zzero,y,desc_data,trans_,work_,info) if(info /= 0) then call psb_errpush(4010,name,a_err='psb_zmlprc_aply') goto 9999 end if else if (size(prec%baseprecv) == 1) then - call psb_baseprc_aply(prec%baseprecv(1),x,zzero,y,desc_data,trans_, work_,info) + call psb_baseprc_aply(zone,prec%baseprecv(1),x,zzero,y,desc_data,trans_, work_,info) else write(0,*) 'Inconsistent preconditioner: size of baseprecv???' endif diff --git a/src/prec/psb_zprecbld.f90 b/src/prec/psb_zprecbld.f90 index 2b17bcc7..8447a6d9 100644 --- a/src/prec/psb_zprecbld.f90 +++ b/src/prec/psb_zprecbld.f90 @@ -49,7 +49,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) Implicit None type(psb_zspmat_type), target :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(in), target :: desc_a type(psb_zprec_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd @@ -60,7 +60,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) use psb_descriptor_type use psb_prec_type type(psb_zspmat_type), target :: a - type(psb_desc_type), intent(in) :: desc_a + type(psb_desc_type), intent(in), target :: desc_a type(psb_zbaseprc_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd @@ -149,29 +149,8 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) endif if (size(p%baseprecv) > 1) then - call init_baseprc_av(p%baseprecv(2),info) - if (info /= 0) then - info=4010 - ch_err='allocate' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - - call psb_mlprc_bld(a,desc_a,p%baseprecv(2),info) - if(info /= 0) then - info=4010 - ch_err='psb_mlprc_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - ! - ! Note: this here has not been tried out. We probably need - ! a different baseprc field %desc_ac, in case we try RAS on lev. 2 of - ! a 3-level prec. - ! - do i=3, size(p%baseprecv) + do i=2, size(p%baseprecv) call init_baseprc_av(p%baseprecv(i),info) if (info /= 0) then info=4010 @@ -180,10 +159,19 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) goto 9999 endif - call psb_mlprc_bld(p%baseprecv(i-1)%av(ac_),p%baseprecv(i-1)%desc_data,& + call psb_mlprc_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,& & p%baseprecv(i),info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name) + goto 9999 + endif + if (debug) then + write(0,*) 'Return from ',i-1,' call to mlprcbld ',info + endif + end do - + endif call psb_erractionrestore(err_act) diff --git a/src/serial/dp/check_dim.f b/src/serial/dp/check_dim.f index 707d23a2..51b569d4 100644 --- a/src/serial/dp/check_dim.f +++ b/src/serial/dp/check_dim.f @@ -51,7 +51,7 @@ C Local scalars NNZ = NZ - LIMIT = INT(DIM_BLOCK*PSB_PERCENT_) +c$$$ LIMIT = INT(DIM_BLOCK*PSB_PERCENT_) DO BLOCK = 1, NG DIM_BLOCK = IA(1,BLOCK+1)-IA(1,BLOCK) diff --git a/src/serial/f77/smmp.f b/src/serial/f77/smmp.f index 8b42405f..8665e980 100644 --- a/src/serial/f77/smmp.f +++ b/src/serial/f77/smmp.f @@ -24,10 +24,21 @@ c * index(*) integer, pointer :: ic(:),jc(:) integer :: nze, info + integer, save :: iunit=11 c c symbolic matrix multiply c=a*b c -c$$$ write(0,*) 'SYMBMM: ',n,m,l,ib(m+1)-1,jb(ib(m+1)-1) +c$$$ open(iunit) +c$$$ write(iunit,*) 'SYMBMM: ',n,m,l +c$$$ write(iunit,*) 'SYMBMM: A:' +c$$$ do i=1,n +c$$$ write(iunit,*) 'Row:',i,' : ',ja(ia(i):ia(i+1)-1) +c$$$ enddo +c$$$ write(iunit,*) 'SYMBMM: B:' +c$$$ do i=1,m +c$$$ write(iunit,*) 'Row:',i,' : ',jb(ib(i):ib(i+1)-1) +c$$$ enddo + if (size(ic) < n+1) then call psb_realloc(n+1,ic,info) endif @@ -71,7 +82,7 @@ c b = d + ... endif do 20 k=ib(j),ib(j+1)-1 if ((jb(k)<1).or.(jb(k)>maxlmn)) then - write(0,*) 'Problem in SYMBMM 1:',j,k,jb(k),maxlmn + write(0,*) 'Problem in SYMBMM 1:',j,k,jb(k),maxlmn else if(index(jb(k)).eq.0) then index(jb(k))=istart @@ -95,6 +106,7 @@ c endif call psb_realloc(nze,jc,info) end if + do 40 j= ic(i),ic(i+1)-1 if (diagc.eq.1 .and. istart.eq.i) then istart = index(istart) @@ -105,8 +117,11 @@ c index(jc(j))=0 40 continue call isr(length,jc(ic(i))) +c$$$ write(iunit,*) length,' : ',jc(ic(i):ic(i)+length-1) index(i) = 0 50 continue +c$$$ close(iunit) +c$$$ iunit = iunit + 1 c$$$ write(0,*) 'SYMBMM: on exit',ic(n+1)-1,jc(ic(n+1)-1) return end diff --git a/src/serial/jad/djdnrmi.f b/src/serial/jad/djdnrmi.f index 8cec70aa..dd260731 100644 --- a/src/serial/jad/djdnrmi.f +++ b/src/serial/jad/djdnrmi.f @@ -49,7 +49,7 @@ C .. External routines .. PNG = IA(1) PIA = IA(2) PJA = IA(3) - + djdnrmi = 0.d0 IF (DESCRA(1:1).EQ.'G') THEN DJDNRMI = DJADNR(TRANS,M,N,IA(PNG), + A,JA,IA(PJA),IA(PIA), diff --git a/src/serial/psb_cest.f90 b/src/serial/psb_cest.f90 index b7d9eb8d..ed01deae 100644 --- a/src/serial/psb_cest.f90 +++ b/src/serial/psb_cest.f90 @@ -66,7 +66,7 @@ subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, iup, info) lar = nnz else info = 136 - call psb_errpush(info,name,a_err=afmt) + call psb_errpush(info,name,a_err=toupper(afmt)) goto 9999 endif @@ -86,7 +86,7 @@ subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, iup, info) lar = nnz else info = 136 - call psb_errpush(info,name,a_err=afmt) + call psb_errpush(info,name,a_err=toupper(afmt)) goto 9999 endif diff --git a/src/serial/psb_dcsprt.f90 b/src/serial/psb_dcsprt.f90 index a1ac3c24..10dc1598 100644 --- a/src/serial/psb_dcsprt.f90 +++ b/src/serial/psb_dcsprt.f90 @@ -70,11 +70,12 @@ subroutine psb_dcsprt(iout,a,iv,eirs,eics,head,ivr,ivc) else ics = 0 endif - + open(iout) if (present(head)) then write(iout,'(a)') '%%MatrixMarket matrix coordinate real general' write(iout,'(a,a)') '% ',head - write(iout,'(a)') '%' + write(iout,'(a)') '%' + write(iout,'(a,a)') '% ',toupper(a%fida) endif select case(toupper(a%fida)) @@ -181,4 +182,5 @@ subroutine psb_dcsprt(iout,a,iv,eirs,eics,head,ivr,ivc) case default write(0,*) 'Feeling lazy today, format not implemented: "',a%fida,'"' end select + close(iout) end subroutine psb_dcsprt diff --git a/src/serial/psb_dsymbmm.f90 b/src/serial/psb_dsymbmm.f90 index 67e1a092..3db814dd 100644 --- a/src/serial/psb_dsymbmm.f90 +++ b/src/serial/psb_dsymbmm.f90 @@ -38,8 +38,9 @@ ! rewritten in Fortran 95 making use of our sparse matrix facilities. ! -subroutine psb_dsymbmm(a,b,c) +subroutine psb_dsymbmm(a,b,c,info) use psb_spmat_type + use psb_string_mod implicit none type(psb_dspmat_type) :: a,b,c @@ -52,9 +53,8 @@ subroutine psb_dsymbmm(a,b,c) integer n,m,l, ia(*), ja(*), diaga, ib(*), jb(*), diagb,& & diagc, index(*) integer, pointer :: ic(:),jc(:) - end subroutine symbmm + end subroutine symbmm end interface - interface psb_sp_getrow subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw) use psb_spmat_type @@ -69,6 +69,28 @@ subroutine psb_dsymbmm(a,b,c) end subroutine psb_dspgetrow end interface + character(len=20) :: name, ch_err + integer :: err_act + name='psb_symbmm' + call psb_erractionsave(err_act) + + select case(toupper(a%fida(1:3))) + case ('CSR') + case default + info=135 + ch_err=a%fida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select + select case(toupper(b%fida(1:3))) + case ('CSR') + case default + info=136 + ch_err=b%fida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select + if (b%m /= a%k) then write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k endif @@ -78,7 +100,6 @@ subroutine psb_dsymbmm(a,b,c) endif nze = max(a%m+1,2*a%m) call psb_sp_reall(c,nze,info) -!!$ write(0,*) 'SYMBMM90 ',size(c%pl),size(c%pr) ! ! Note: we need to test whether there is a performance impact ! in not using the original Douglas & Bank code. @@ -97,6 +118,15 @@ subroutine psb_dsymbmm(a,b,c) c%fida='CSR' c%descra='GUN' deallocate(itemp) + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if return contains diff --git a/src/serial/psb_zsymbmm.f90 b/src/serial/psb_zsymbmm.f90 index 4b7d23f7..2609b97b 100644 --- a/src/serial/psb_zsymbmm.f90 +++ b/src/serial/psb_zsymbmm.f90 @@ -38,8 +38,9 @@ ! rewritten in Fortran 95 making use of our sparse matrix facilities. ! -subroutine psb_zsymbmm(a,b,c) +subroutine psb_zsymbmm(a,b,c,info) use psb_spmat_type + use psb_string_mod implicit none type(psb_zspmat_type) :: a,b,c @@ -68,6 +69,27 @@ subroutine psb_zsymbmm(a,b,c) integer, intent(out) :: info end subroutine psb_zspgetrow end interface + character(len=20) :: name, ch_err + integer :: err_act + name='psb_symbmm' + call psb_erractionsave(err_act) + + select case(toupper(a%fida(1:3))) + case ('CSR') + case default + info=135 + ch_err=a%fida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select + select case(toupper(b%fida(1:3))) + case ('CSR') + case default + info=136 + ch_err=b%fida(1:3) + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end select if (b%m /= a%k) then write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k @@ -97,6 +119,15 @@ subroutine psb_zsymbmm(a,b,c) c%fida='CSR' c%descra='GUN' deallocate(itemp) + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.act_abort) then + call psb_error() + return + end if return contains diff --git a/src/tools/psb_cdalv.f90 b/src/tools/psb_cdalv.f90 index 2f0f825e..207cf615 100644 --- a/src/tools/psb_cdalv.f90 +++ b/src/tools/psb_cdalv.f90 @@ -235,7 +235,6 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag) endif ! estimate local cols number -!!$ loc_col = int((psb_colrow_+1.d0)*loc_row)+1 loc_col = min(2*loc_row,m) allocate(desc_a%loc_to_glob(loc_col),& &desc_a%lprm(1),stat=info) diff --git a/test/Fileread/df_sample.f90 b/test/Fileread/df_sample.f90 index af813da9..0c75dd5f 100644 --- a/test/Fileread/df_sample.f90 +++ b/test/Fileread/df_sample.f90 @@ -245,7 +245,6 @@ program df_sample call psb_precset(pre,'asm',info,iv=(/novr,halo_,none_/)) case(rash_) call psb_precset(pre,'asm',info,iv=(/novr,nohalo_,none_/)) - case default call psb_precset(pre,'ilu',info) end select