From 580605c087b2cacc440e31e6754face8a95ad977 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 8 Mar 2006 17:11:02 +0000 Subject: [PATCH] Completed reorganization of preconditioner build routines. Now it should even work with multilevel L>2, even though we don't have a test program right now. Started modifying the prc_aply set. --- src/methd/psb_dbicg.f90 | 4 +- src/methd/psb_dcg.f90 | 2 +- src/methd/psb_dcgs.f90 | 4 +- src/methd/psb_dcgstab.f90 | 6 +- src/methd/psb_dcgstabl.f90 | 6 +- src/methd/psb_dgmresr.f90 | 4 +- src/modules/psb_prec_mod.f90 | 25 +-- src/modules/psb_prec_type.f90 | 58 +++---- src/prec/Makefile | 2 +- src/prec/psb_dbaseprc_bld.f90 | 308 ++++++++++++++++++++++++++++++++++ src/prec/psb_dbldaggrmat.f90 | 2 +- src/prec/psb_dilu_bld.f90 | 10 +- src/prec/psb_dmlprc_bld.f90 | 134 +++++---------- src/prec/psb_dprec.f90 | 186 ++++++++++---------- src/prec/psb_dprecbld.f90 | 304 +++++++++------------------------ src/prec/psb_dprecset.f90 | 2 +- src/prec/psb_dslu_bld.f90 | 8 +- src/prec/psb_dsp_renum.f90 | 8 +- src/prec/psb_dumf_bld.f90 | 8 +- 19 files changed, 593 insertions(+), 488 deletions(-) create mode 100644 src/prec/psb_dbaseprc_bld.f90 diff --git a/src/methd/psb_dbicg.f90 b/src/methd/psb_dbicg.f90 index 5c5a5614..b83f3543 100644 --- a/src/methd/psb_dbicg.f90 +++ b/src/methd/psb_dbicg.f90 @@ -270,8 +270,8 @@ subroutine psb_dbicg(a,prec,b,x,eps,desc_a,info,& itx = itx + 1 if (debug) write(*,*) 'iteration: ',itx - call psb_prcaply(prec,r,z,desc_a,info,work=aux) - call psb_prcaply(prec,rt,zt,desc_a,info,trans='t',work=aux) + call psb_prc_aply(prec,r,z,desc_a,info,work=aux) + call psb_prc_aply(prec,rt,zt,desc_a,info,trans='t',work=aux) rho_old = rho rho = psb_dot(rt,z,desc_a,info) diff --git a/src/methd/psb_dcg.f90 b/src/methd/psb_dcg.f90 index 9357fede..f45a2ba1 100644 --- a/src/methd/psb_dcg.f90 +++ b/src/methd/psb_dcg.f90 @@ -215,7 +215,7 @@ Subroutine psb_dcg(a,prec,b,x,eps,desc_a,info,& it = it + 1 itx = itx + 1 - Call psb_prcaply(prec,r,z,desc_a,info,work=aux) + Call psb_prc_aply(prec,r,z,desc_a,info,work=aux) rho_old = rho rho = psb_dot(r,z,desc_a,info) diff --git a/src/methd/psb_dcgs.f90 b/src/methd/psb_dcgs.f90 index df4e9cff..04bcc311 100644 --- a/src/methd/psb_dcgs.f90 +++ b/src/methd/psb_dcgs.f90 @@ -274,7 +274,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& End If - Call psb_prcaply(prec,p,f,desc_a,info,work=aux) + Call psb_prc_aply(prec,p,f,desc_a,info,work=aux) Call psb_spmm(one,a,f,zero,v,desc_a,info,& & work=aux) @@ -292,7 +292,7 @@ Subroutine psb_dcgs(a,prec,b,x,eps,desc_a,info,& Call psb_axpby(one,uv,zero,s,desc_a,info) Call psb_axpby(one,q,one,s,desc_a,info) - Call psb_prcaply(prec,s,z,desc_a,info,work=aux) + Call psb_prc_aply(prec,s,z,desc_a,info,work=aux) Call psb_axpby(alpha,z,one,x,desc_a,info) diff --git a/src/methd/psb_dcgstab.f90 b/src/methd/psb_dcgstab.f90 index fd53031a..dd227259 100644 --- a/src/methd/psb_dcgstab.f90 +++ b/src/methd/psb_dcgstab.f90 @@ -300,7 +300,7 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,& Call psb_axpby(one,r,beta,p,desc_a,info) End If - Call psb_prcaply(prec,p,f,desc_a,info,work=aux) + Call psb_prc_aply(prec,p,f,desc_a,info,work=aux) Call psb_spmm(one,a,f,zero,v,desc_a,info,& & work=aux) @@ -323,9 +323,9 @@ Subroutine psb_dcgstab(a,prec,b,x,eps,desc_a,info,& goto 9999 end if - Call psb_prcaply(prec,s,z,desc_a,info,work=aux) + Call psb_prc_aply(prec,s,z,desc_a,info,work=aux) if(info.ne.0) then - call psb_errpush(4010,name,a_err='psb_prcaply') + call psb_errpush(4010,name,a_err='psb_prc_aply') goto 9999 end if diff --git a/src/methd/psb_dcgstabl.f90 b/src/methd/psb_dcgstabl.f90 index bb2d0fc1..678a7aa7 100644 --- a/src/methd/psb_dcgstabl.f90 +++ b/src/methd/psb_dcgstabl.f90 @@ -237,7 +237,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& Call psb_axpby(one,b,zero,r,desc_a,info) Call psb_spmm(-one,a,x,one,r,desc_a,info,work=aux) - call psb_prcaply(prec,r,desc_a,info) + call psb_prc_aply(prec,r,desc_a,info) Call psb_axpby(one,r,zero,rt0,desc_a,info) Call psb_axpby(one,r,zero,rh(:,0),desc_a,info) @@ -301,7 +301,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& If (debug) Write(0,*) 'bicg part: ',rh(1,0),beta Call psb_spmm(one,a,uh(:,j),zero,uh(:,j+1),desc_a,info,work=aux) - call psb_prcaply(prec,uh(:,j+1),desc_a,info) + call psb_prc_aply(prec,uh(:,j+1),desc_a,info) gamma(j) = psb_dot(uh(:,j+1),rt0,desc_a,info) If (gamma(j)==zero) Then @@ -315,7 +315,7 @@ Subroutine psb_dcgstabl(a,prec,b,x,eps,desc_a,info,& Call psb_axpby(alpha,uh(:,0),one,x,desc_a,info) Call psb_spmm(one,a,rh(:,j),zero,rh(:,j+1),desc_a,info,work=aux) - call psb_prcaply(prec,rh(:,j+1),desc_a,info) + call psb_prc_aply(prec,rh(:,j+1),desc_a,info) Enddo diff --git a/src/methd/psb_dgmresr.f90 b/src/methd/psb_dgmresr.f90 index a0515828..031a8a18 100644 --- a/src/methd/psb_dgmresr.f90 +++ b/src/methd/psb_dgmresr.f90 @@ -234,7 +234,7 @@ Subroutine psb_dgmresr(a,prec,b,x,eps,desc_a,info,& End If Call psb_spmm(-one,a,x,one,v(:,1),desc_a,info,work=aux) - call psb_prcaply(prec,v(:,1),desc_a,info) + call psb_prc_aply(prec,v(:,1),desc_a,info) rs(1) = psb_nrm2(v(:,1),desc_a,info) if (info.ne.0) Then info=4011 @@ -279,7 +279,7 @@ Subroutine psb_dgmresr(a,prec,b,x,eps,desc_a,info,& itx = itx + 1 Call psb_spmm(one,a,v(:,i),zero,w,desc_a,info,work=aux) - call psb_prcaply(prec,w,desc_a,info) + call psb_prc_aply(prec,w,desc_a,info) do k = 1, i h(k,i) = psb_dot(v(:,k),w,desc_a,info) diff --git a/src/modules/psb_prec_mod.f90 b/src/modules/psb_prec_mod.f90 index 23af25b1..a837c322 100644 --- a/src/modules/psb_prec_mod.f90 +++ b/src/modules/psb_prec_mod.f90 @@ -40,7 +40,7 @@ module psb_prec_mod type(psb_dspmat_type), intent(in), target :: a type(psb_desc_type), intent(in) :: desc_a type(psb_dspmat_type), intent(out), target :: ac - type(psb_dbase_prec), intent(inout) :: p + type(psb_dbaseprc_type), intent(inout) :: p type(psb_desc_type), intent(inout) :: desc_p integer, intent(out) :: info end subroutine psb_dbldaggrmat @@ -100,19 +100,6 @@ end interface end subroutine psb_dprecfree end interface - interface psb_cslu - subroutine psb_dcslu(a,desc_data,p,upd,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbase_prec), intent(inout) :: p - character, intent(in) :: upd - end subroutine psb_dcslu - end interface - interface psb_asmatbld Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_serial_mod @@ -129,8 +116,8 @@ end interface end Subroutine psb_dasmatbld end interface - interface psb_prcaply - subroutine psb_dprecaply(prec,x,y,desc_data,info,trans,work) + interface psb_prc_aply + subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans,work) use psb_serial_mod use psb_descriptor_type use psb_prec_type @@ -140,8 +127,8 @@ end interface integer, intent(out) :: info character(len=1), optional :: trans real(kind(0.d0)),intent(inout), optional, target :: work(:) - end subroutine psb_dprecaply - subroutine psb_dprecaply1(prec,x,desc_data,info,trans) + end subroutine psb_dprc_aply + subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) use psb_serial_mod use psb_descriptor_type use psb_prec_type @@ -150,7 +137,7 @@ end interface real(kind(0.d0)),intent(inout) :: x(:) integer, intent(out) :: info character(len=1), optional :: trans - end subroutine psb_dprecaply1 + end subroutine psb_dprc_aply1 end interface diff --git a/src/modules/psb_prec_type.f90 b/src/modules/psb_prec_type.f90 index 128b2d29..1f9ff203 100644 --- a/src/modules/psb_prec_type.f90 +++ b/src/modules/psb_prec_type.f90 @@ -37,7 +37,7 @@ module psb_prec_type use psb_spmat_type use psb_descriptor_type - + integer, parameter :: min_prec_=0, noprec_=0, diagsc_=1, bja_=2,& & asm_=3, ras_=5, ash_=4, rash_=6, ras2lv_=7, ras2lvm_=8,& & lv2mras_=9, lv2smth_=10, lv2lsm_=11, sl2sm_=12, superlu_=13,& @@ -82,7 +82,7 @@ module psb_prec_type integer, parameter :: smth_avsz=6, max_avsz=smth_avsz - type psb_dbase_prec + type psb_dbaseprc_type type(psb_dspmat_type), pointer :: av(:) => null() ! real(kind(1.d0)), pointer :: d(:) => null() @@ -93,16 +93,16 @@ module psb_prec_type integer, pointer :: mlia(:) => null(), nlaggr(:) => null() ! type(psb_dspmat_type), pointer :: aorig => null() ! real(kind(1.d0)), pointer :: dorig(:) => null() ! - - end type psb_dbase_prec - + + end type psb_dbaseprc_type + type psb_dprec_type - type(psb_dbase_prec), pointer :: baseprecv(:) => null() + type(psb_dbaseprc_type), pointer :: baseprecv(:) => null() ! contain type of preconditioning to be performed integer :: prec, base_prec - end type psb_dprec_type + end type psb_dprec_type + - character(len=15), parameter, private :: & & smooth_names(1:3)=(/'Pre-smoothing ','Post-smoothing',& & 'Smooth both '/) @@ -139,7 +139,7 @@ module psb_prec_type interface psb_prec_short_descr module procedure psb_prec_short_descr end interface - + contains subroutine psb_file_prec_descr(iout,p) @@ -282,14 +282,14 @@ contains function is_legal_base_prec(ip) integer, intent(in) :: ip logical :: is_legal_base_prec - + is_legal_base_prec = ((ip>=noprec_).and.(ip<=rash_)) return end function is_legal_base_prec function is_legal_n_ovr(ip) integer, intent(in) :: ip logical :: is_legal_n_ovr - + is_legal_n_ovr = (ip >=0) return end function is_legal_n_ovr @@ -303,89 +303,89 @@ contains function is_legal_jac_sweeps(ip) integer, intent(in) :: ip logical :: is_legal_jac_sweeps - + is_legal_jac_sweeps = (ip >= 1) return end function is_legal_jac_sweeps function is_legal_prolong(ip) integer, intent(in) :: ip logical :: is_legal_prolong - + is_legal_prolong = ((ip>=none_).and.(ip<=square_root_)) return end function is_legal_prolong function is_legal_restrict(ip) integer, intent(in) :: ip logical :: is_legal_restrict - + is_legal_restrict = ((ip==nohalo_).or.(ip==halo_)) return end function is_legal_restrict function is_legal_ml_type(ip) integer, intent(in) :: ip logical :: is_legal_ml_type - + is_legal_ml_type = ((ip>=no_ml_).and.(ip<=max_ml_)) return end function is_legal_ml_type function is_legal_ml_aggr_kind(ip) integer, intent(in) :: ip logical :: is_legal_ml_aggr_kind - + is_legal_ml_aggr_kind = ((ip>=loc_aggr_).and.(ip<=max_aggr_)) return end function is_legal_ml_aggr_kind function is_legal_ml_smooth_pos(ip) integer, intent(in) :: ip logical :: is_legal_ml_smooth_pos - + is_legal_ml_smooth_pos = ((ip>=pre_smooth_).and.(ip<=max_smooth_)) return end function is_legal_ml_smooth_pos function is_legal_ml_smth_kind(ip) integer, intent(in) :: ip logical :: is_legal_ml_smth_kind - + is_legal_ml_smth_kind = ((ip>=no_smth_).and.(ip<=smth_biz_)) return end function is_legal_ml_smth_kind function is_legal_ml_coarse_mat(ip) integer, intent(in) :: ip logical :: is_legal_ml_coarse_mat - + is_legal_ml_coarse_mat = ((ip>=mat_distr_).and.(ip<=mat_repl_)) return end function is_legal_ml_coarse_mat function is_legal_ml_fact(ip) integer, intent(in) :: ip logical :: is_legal_ml_fact - + is_legal_ml_fact = ((ip>=f_ilu_n_).and.(ip<=f_umf_)) return end function is_legal_ml_fact function is_legal_ml_lev(ip) integer, intent(in) :: ip logical :: is_legal_ml_lev - + is_legal_ml_lev = (ip>=0) return end function is_legal_ml_lev function is_legal_omega(ip) real(kind(1.d0)), intent(in) :: ip logical :: is_legal_omega - + is_legal_omega = ((ip>=0.0d0).and.(ip<=2.0d0)) return end function is_legal_omega function is_legal_ml_eps(ip) real(kind(1.d0)), intent(in) :: ip logical :: is_legal_ml_eps - + is_legal_ml_eps = (ip>=0.0d0) return end function is_legal_ml_eps - + subroutine psb_icheck_def(ip,name,id,is_legal) integer, intent(inout) :: ip integer, intent(in) :: id @@ -424,7 +424,7 @@ contains use psb_serial_mod use psb_descriptor_type use psb_tools_mod - type(psb_dbase_prec), intent(inout) :: p + type(psb_dbaseprc_type), intent(inout) :: p integer, intent(out) :: info integer :: i @@ -454,7 +454,7 @@ contains endif if (associated(p%dprcparm)) then deallocate(p%dprcparm,stat=info) - end if + end if if (associated(p%aorig)) then ! This is a pointer to something else, must not free it here. nullify(p%aorig) @@ -494,11 +494,11 @@ contains subroutine psb_nullify_baseprec(p) use psb_descriptor_type - type(psb_dbase_prec), intent(inout) :: p - + 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) - + end subroutine psb_nullify_baseprec function pr_to_str(iprec) diff --git a/src/prec/Makefile b/src/prec/Makefile index 24e79fd8..28f5d0a8 100644 --- a/src/prec/Makefile +++ b/src/prec/Makefile @@ -7,7 +7,7 @@ MPFOBJS=psb_dilu_bld.o psb_dbldaggrmat.o F90OBJS=psb_dasmatbld.o psb_dslu_bld.o psb_dumf_bld.o psb_dilu_fct.o\ psb_dmlprc_bld.o psb_dsp_renum.o\ psb_dprec.o psb_dprecbld.o gps.o psb_dprecfree.o psb_dprecset.o \ - psb_dgenaggrmap.o $(MPFOBJS) + psb_dbaseprc_bld.o psb_dgenaggrmap.o $(MPFOBJS) COBJS=psb_slu_impl.o psb_umf_impl.o INCDIRS=-I. -I.. -I$(LIBDIR) diff --git a/src/prec/psb_dbaseprc_bld.f90 b/src/prec/psb_dbaseprc_bld.f90 new file mode 100644 index 00000000..14911b43 --- /dev/null +++ b/src/prec/psb_dbaseprc_bld.f90 @@ -0,0 +1,308 @@ +!!$ +!!$ +!!$ MPcube: Multilevel Parallel Preconditioners Package +!!$ for +!!$ Parallel Sparse BLAS v2.0 +!!$ (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari University of Rome Tor Vergata +!!$ Daniela Di Serafino II University of Naples +!!$ Pasqua D'Ambra ICAR-CNR +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MPCUBE group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MPCUBE GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) + + use psb_serial_mod + Use psb_spmat_type + use psb_descriptor_type + use psb_prec_type + use psb_tools_mod + use psb_comm_mod + use psb_const_mod + use psb_psblas_mod + use psb_error_mod + Implicit None + + type(psb_dspmat_type), target :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type),intent(inout) :: p + integer, intent(out) :: info + character, intent(in), optional :: upd + + interface psb_ilu_bld + subroutine psb_dilu_bld(a,desc_data,p,upd,info) + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + integer, intent(out) :: info + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type),intent(in) :: desc_data + type(psb_dbaseprc_type), intent(inout) :: p + character, intent(in) :: upd + end subroutine psb_dilu_bld + end interface + + interface psb_slu_bld + subroutine psb_dslu_bld(a,desc_a,p,info) + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_const_mod + implicit none + + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type), intent(inout) :: p + integer, intent(out) :: info + end subroutine psb_dslu_bld + end interface + + interface psb_umf_bld + subroutine psb_dumf_bld(a,desc_a,p,info) + use psb_serial_mod + use psb_descriptor_type + use psb_prec_type + use psb_const_mod + implicit none + + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type), intent(inout) :: p + integer, intent(out) :: info + end subroutine psb_dumf_bld + end interface + + ! Local scalars + Integer :: err, nnzero, n_row, n_col,I,j,k,icontxt,& + & me,mycol,nprow,npcol,mglob,lw, mtype, nrg, nzg, err_act + real(kind(1.d0)) :: temp, real_err(5) + real(kind(1.d0)),pointer :: gd(:), work(:) + integer :: int_err(5) + character :: iupd + + logical, parameter :: debug=.false. + integer,parameter :: iroot=0,iout=60,ilout=40 + character(len=20) :: name, ch_err + + if(psb_get_errstatus().ne.0) return + info=0 + err=0 + call psb_erractionsave(err_act) + name = 'psb_baseprc_bld' + + if (debug) write(0,*) 'Entering baseprc_bld' + info = 0 + int_err(1) = 0 + icontxt = desc_a%matrix_data(psb_ctxt_) + n_row = desc_a%matrix_data(psb_n_row_) + n_col = desc_a%matrix_data(psb_n_col_) + mglob = desc_a%matrix_data(psb_m_) + if (debug) write(0,*) 'Preconditioner Blacs_gridinfo' + call blacs_gridinfo(icontxt, nprow, npcol, me, mycol) + + if (present(upd)) then + if (debug) write(0,*) 'UPD ', upd + if ((UPD.eq.'F').or.(UPD.eq.'T')) then + IUPD=UPD + else + IUPD='F' + endif + else + IUPD='F' + endif + + ! + ! Should add check to ensure all procs have the same... + ! + ! ALso should define symbolic names for the preconditioners. + ! + + call psb_check_def(p%iprcparm(p_type_),'base_prec',& + & diagsc_,is_legal_base_prec) + + allocate(p%desc_data,stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + call psb_nullify_desc(p%desc_data) + + select case(p%iprcparm(p_type_)) + case (noprec_) + ! Do nothing. + + + case (diagsc_) + + if (debug) write(0,*) 'Precond: Diagonal scaling' + ! diagonal scaling + + call psb_realloc(n_col,p%d,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_realloc') + goto 9999 + end if + + call psb_csrws(p%d,a,info,trans='N') + if(info /= 0) then + info=4010 + ch_err='psb_csrws' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug) write(ilout+me,*) 'VDIAG ',n_row + do i=1,n_row + if (p%d(i).eq.0.0d0) then + p%d(i)=1.d0 + else + p%d(i) = 1.d0/p%d(i) + endif + + if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i), p%d(i) + if (p%d(i).lt.0.d0) then + write(0,*) me,'Negative RWS? ',i,p%d(i) + endif + end do + if (a%pl(1) /= 0) then + allocate(work(n_row),stat=info) + if (info /= 0) then + info=4000 + call psb_errpush(info,name) + goto 9999 + end if + call psb_gelp('n',a%pl,p%d,desc_a,info) + if(info /= 0) then + info=4010 + ch_err='psb_dgelp' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + deallocate(work) + endif + + if (debug) then + allocate(gd(mglob),stat=info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='Allocate') + goto 9999 + end if + + call psb_gather(gd, p%d, desc_a, info, iroot=iroot) + if(info /= 0) then + info=4010 + ch_err='psb_dgatherm' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (me.eq.iroot) then + write(iout+nprow,*) 'VDIAG CHECK ',mglob + do i=1,mglob + write(iout+nprow,*) i,gd(i) + enddo + endif + deallocate(gd) + endif + if (debug) write(*,*) 'Preconditioner DIAG computed OK' + + + case (bja_,asm_) + + call psb_check_def(p%iprcparm(n_ovr_),'overlap',& + & 0,is_legal_n_ovr) + call psb_check_def(p%iprcparm(restr_),'restriction',& + & psb_halo_,is_legal_restrict) + call psb_check_def(p%iprcparm(prol_),'prolongator',& + & psb_none_,is_legal_prolong) + call psb_check_def(p%iprcparm(iren_),'renumbering',& + & renum_none_,is_legal_renum) + call psb_check_def(p%iprcparm(f_type_),'fact',& + & f_ilu_n_,is_legal_ml_fact) + + if (debug) write(0,*)me, ': Calling PSB_ILU_BLD' + + + select case(p%iprcparm(f_type_)) + + case(f_ilu_n_,f_ilu_e_) + call psb_ilu_bld(a,desc_a,p,iupd,info) + if(debug) write(0,*)me,': out of psb_ilu_bld' + if(info /= 0) then + info=4010 + ch_err='psb_ilu_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(f_slu_) + + if(debug) write(0,*)me,': calling slu_bld' + call psb_slu_bld(a,desc_a,p,info) + if(info /= 0) then + info=4010 + ch_err='slu_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(f_umf_) + if(debug) write(0,*)me,': calling umf_bld' + call psb_umf_bld(a,desc_a,p,info) + if(info /= 0) then + info=4010 + ch_err='umf_bld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(f_none_) + write(0,*) 'Fact=None in BASEPRC_BLD Bja/ASM??' + + case default + write(0,*) 'Unknown factor type in baseprc_bld bja/asm: ',& + &p%iprcparm(f_type_) + end select + + end select + + + + 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 + +end subroutine psb_dbaseprc_bld + diff --git a/src/prec/psb_dbldaggrmat.f90 b/src/prec/psb_dbldaggrmat.f90 index eb52c3ac..49484ca7 100644 --- a/src/prec/psb_dbldaggrmat.f90 +++ b/src/prec/psb_dbldaggrmat.f90 @@ -44,7 +44,7 @@ subroutine psb_dbldaggrmat(a,desc_a,ac,p,desc_p,info) implicit none type(psb_dspmat_type), intent(in), target :: a - type(psb_dbase_prec), intent(inout) :: p + type(psb_dbaseprc_type), intent(inout) :: p type(psb_dspmat_type), intent(out), target :: ac type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(inout) :: desc_p diff --git a/src/prec/psb_dilu_bld.f90 b/src/prec/psb_dilu_bld.f90 index c369d538..bf20e340 100644 --- a/src/prec/psb_dilu_bld.f90 +++ b/src/prec/psb_dilu_bld.f90 @@ -63,7 +63,7 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info) integer, intent(out) :: info ! .. array Arguments .. type(psb_dspmat_type), intent(in), target :: a - type(psb_dbase_prec), intent(inout) :: p + type(psb_dbaseprc_type), intent(inout) :: p type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd @@ -118,10 +118,10 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info) implicit none ! .. array Arguments .. - type(psb_dspmat_type), intent(in) :: a,blck - type(psb_dspmat_type), intent(inout) :: atmp - type(psb_dbase_prec), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(in) :: a,blck + type(psb_dspmat_type), intent(inout) :: atmp + type(psb_dbaseprc_type), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info end subroutine psb_dsp_renum end interface diff --git a/src/prec/psb_dmlprc_bld.f90 b/src/prec/psb_dmlprc_bld.f90 index 53073dd4..c3f900a8 100644 --- a/src/prec/psb_dmlprc_bld.f90 +++ b/src/prec/psb_dmlprc_bld.f90 @@ -45,7 +45,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) type(psb_dspmat_type), intent(in), target :: a type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p + type(psb_dbaseprc_type), intent(inout) :: p integer, intent(out) :: info type(psb_desc_type), pointer :: desc_p @@ -55,15 +55,17 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) logical, parameter :: debug=.false. type(psb_dspmat_type) :: ac - interface psb_ilu_fct - subroutine psb_dilu_fct(a,l,u,d,info,blck) - use psb_spmat_type - integer, intent(out) :: info - type(psb_dspmat_type),intent(in) :: a - type(psb_dspmat_type),intent(inout) :: l,u - type(psb_dspmat_type),intent(in), optional, target :: blck - real(kind(1.d0)), intent(inout) :: d(:) - end subroutine psb_dilu_fct + 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), target :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type),intent(inout) :: p + integer, intent(out) :: info + character, intent(in), optional :: upd + end subroutine psb_dbaseprc_bld end interface interface psb_genaggrmap @@ -85,7 +87,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) use psb_descriptor_type use psb_spmat_type type(psb_dspmat_type), intent(in), target :: a - type(psb_dbase_prec), intent(inout) :: p + type(psb_dbaseprc_type), intent(inout) :: p type(psb_dspmat_type), intent(out),target :: ac type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(inout) :: desc_p @@ -93,49 +95,6 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) end subroutine psb_dbldaggrmat end interface - interface psb_ilu_bld - subroutine psb_dilu_bld(a,desc_data,p,upd,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbase_prec), intent(inout) :: p - character, intent(in) :: upd - end subroutine psb_dilu_bld - end interface - - interface psb_slu_bld - subroutine psb_dslu_bld(a,desc_a,p,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_const_mod - implicit none - - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info - end subroutine psb_dslu_bld - end interface - - interface psb_umf_bld - subroutine psb_dumf_bld(a,desc_a,p,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_const_mod - implicit none - - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info - end subroutine psb_dumf_bld - end interface - integer :: icontxt, nprow, npcol, me, mycol name='psb_mlprec_bld' @@ -144,6 +103,37 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) call psb_erractionsave(err_act) call psb_nullify_sp(ac) + + if (.not.associated(p%iprcparm)) then + info = 2222 + call psb_errpush(info,name) + goto 9999 + endif + call psb_check_def(p%iprcparm(ml_type_),'Multilevel type',& + & mult_ml_prec_,is_legal_ml_type) + call psb_check_def(p%iprcparm(aggr_alg_),'aggregation',& + & loc_aggr_,is_legal_ml_aggr_kind) + call psb_check_def(p%iprcparm(smth_kind_),'Smoother kind',& + & smth_omg_,is_legal_ml_smth_kind) + call psb_check_def(p%iprcparm(coarse_mat_),'Coarse matrix',& + & mat_distr_,is_legal_ml_coarse_mat) + call psb_check_def(p%iprcparm(smth_pos_),'smooth_pos',& + & pre_smooth_,is_legal_ml_smooth_pos) + + + nullify(p%desc_data) + select case(p%iprcparm(f_type_)) + case(f_ilu_n_) + call psb_check_def(p%iprcparm(ilu_fill_in_),'Level',0,is_legal_ml_lev) + case(f_ilu_e_) + call psb_check_def(p%dprcparm(fact_eps_),'Eps',0.0d0,is_legal_ml_eps) + end select + call psb_check_def(p%dprcparm(smooth_omega_),'omega',0.0d0,is_legal_omega) + call psb_check_def(p%iprcparm(jac_sweeps_),'Jacobi sweeps',& + & 1,is_legal_jac_sweeps) + + + p%aorig => a allocate(p%av(max_avsz),stat=info) @@ -180,41 +170,9 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) end if if (debug) write(0,*) 'Out from bldaggrmat',desc_p%matrix_data(:) - allocate(p%desc_data) - - select case(p%iprcparm(f_type_)) - - case(f_ilu_n_,f_ilu_e_) - call psb_ilu_bld(ac,desc_p,p,'F',info) - if(debug) write(0,*)me,': out of psb_ilu_bld' - if(info /= 0) then - info=4010 - ch_err='psb_ilu_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_slu_) - call psb_slu_bld(ac,desc_p,p,info) - if(debug) write(0,*)me,': out of psb_slu_bld' - if(info /= 0) then - info=4010 - ch_err='psb_slu_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_umf_) - call psb_umf_bld(ac,desc_p,p,info) - if(debug) write(0,*)me,': out of psb_umf_bld' - if(info /= 0) then - info=4010 - ch_err='psb_umf_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - end select + + call psb_baseprc_bld(ac,desc_p,p,info) ! ! We have used a separate ac because: diff --git a/src/prec/psb_dprec.f90 b/src/prec/psb_dprec.f90 index e5863e23..66eaf277 100644 --- a/src/prec/psb_dprec.f90 +++ b/src/prec/psb_dprec.f90 @@ -33,7 +33,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) +subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) use psb_serial_mod use psb_descriptor_type @@ -59,35 +59,35 @@ subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) external mpi_wtime character(len=20) :: name, ch_err - interface - subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) + interface psb_baseprc_aply + subroutine psb_dbaseprc_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_dbase_prec), 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_dbaseprcaply + 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_dbaseprc_aply end interface - interface - subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) + interface psb_mlprc_aply + subroutine psb_dmlprc_aply(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_dbase_prec), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: beta - real(kind(0.d0)),intent(inout) :: x(:), y(:) - character :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dmlprcaply + 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(inout) :: x(:), y(:) + character :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + end subroutine psb_dmlprc_aply end interface - name='psb_dprecaply' + name='psb_dprc_aply' info = 0 call psb_erractionsave(err_act) @@ -115,15 +115,15 @@ subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) write(0,*) 'Inconsistent preconditioner: neither SMTH nor BASE?' end if if (size(prec%baseprecv) >1) then - if (debug) write(0,*) 'Into mlprcaply',size(x),size(y) - call psb_dmlprcaply(prec%baseprecv,x,zero,y,desc_data,trans_,work_,info) + if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y) + call psb_mlprc_aply(prec%baseprecv,x,zero,y,desc_data,trans_,work_,info) if(info /= 0) then - call psb_errpush(4010,name,a_err='psb_dmlprcaply') + call psb_errpush(4010,name,a_err='psb_dmlprc_aply') goto 9999 end if else if (size(prec%baseprecv) == 1) then - call psb_dbaseprcaply(prec%baseprecv(1),x,zero,y,desc_data,trans_, work_,info) + call psb_baseprc_aply(prec%baseprecv(1),x,zero,y,desc_data,trans_, work_,info) else write(0,*) 'Inconsistent preconditioner: size of baseprecv???' endif @@ -144,7 +144,7 @@ subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) end if return -end subroutine psb_dprecaply +end subroutine psb_dprc_aply !!$ @@ -183,7 +183,7 @@ end subroutine psb_dprecaply !!$ !!$ -subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) +subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info) ! ! Compute Y <- beta*Y + K^-1 X ! where K is a a basic preconditioner stored in prec @@ -197,13 +197,13 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) use psb_error_mod implicit none - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbase_prec), 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 + 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 ! Local variables integer :: n_row,n_col, int_err(5) @@ -216,21 +216,21 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) external mpi_wtime character(len=20) :: name, ch_err - interface - subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) + 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_dbase_prec), intent(in) :: prec + 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_dbjacaply + end subroutine psb_dbjac_aply end interface - name='psb_dbaseprcaply' + name='psb_dbaseprc_aply' info = 0 call psb_erractionsave(err_act) @@ -280,10 +280,10 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) case(bja_) - call psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) + call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 - ch_err='psb_bjacaply' + ch_err='psb_bjac_aply' goto 9999 end if @@ -291,7 +291,7 @@ subroutine psb_dbaseprcaply(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_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) + call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info) if(info.ne.0) then info=4010 ch_err='psb_bjacaply' @@ -350,7 +350,7 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) goto 9999 end if else if (prec%iprcparm(restr_) /= psb_none_) then - write(0,*) 'Problem in PRCAPLY: Unknown value for restriction ',& + write(0,*) 'Problem in PRC_APLY: Unknown value for restriction ',& &prec%iprcparm(restr_) end if @@ -363,10 +363,10 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) end if endif - call psb_dbjacaply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) + call psb_bjac_aply(prec,tx,zero,ty,prec%desc_data,trans,aux,info) if(info.ne.0) then info=4010 - ch_err='psb_bjacaply' + ch_err='psb_bjac_aply' goto 9999 end if @@ -395,7 +395,7 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) end if case default - write(0,*) 'Problem in PRCAPLY: Unknown value for prolongation ',& + write(0,*) 'Problem in PRC_APLY: Unknown value for prolongation ',& & prec%iprcparm(prol_) end select @@ -440,7 +440,7 @@ subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) end if return -end subroutine psb_dbaseprcaply +end subroutine psb_dbaseprc_aply @@ -479,7 +479,7 @@ end subroutine psb_dbaseprcaply !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) +subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info) ! ! Compute Y <- beta*Y + K^-1 X ! where K is a a Block Jacobi preconditioner stored in prec @@ -496,7 +496,7 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) implicit none type(psb_desc_type), intent(in) :: desc_data - type(psb_dbase_prec), intent(in) :: prec + 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 @@ -514,7 +514,7 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) external mpi_wtime character(len=20) :: name, ch_err - name='psb_dbjacaply' + name='psb_bjac_aply' info = 0 call psb_erractionsave(err_act) @@ -630,7 +630,7 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) endif case default - write(0,*) 'Unknown factorization type in bjac_prcaply',prec%iprcparm(f_type_) + write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_) end select if (debugprt) write(0,*)' Y: ',y(:) @@ -739,7 +739,7 @@ subroutine psb_dbjacaply(prec,x,beta,y,desc_data,trans,work,info) end if return -end subroutine psb_dbjacaply +end subroutine psb_dbjac_aply !!$ @@ -777,7 +777,7 @@ end subroutine psb_dbjacaply !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) +subroutine psb_dmlprc_aply(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 @@ -792,13 +792,13 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) use psb_error_mod implicit none - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbase_prec), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: beta - real(kind(0.d0)),intent(inout) :: x(:), y(:) - character :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: 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(inout) :: x(:), y(:) + character :: trans + real(kind(0.d0)),target :: work(:) + integer, intent(out) :: info ! Local variables @@ -815,21 +815,21 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) external mpi_wtime character(len=20) :: name, ch_err - interface - subroutine psb_dbaseprcaply(prec,x,beta,y,desc_data,trans,work,info) + interface psb_baseprc_aply + subroutine psb_dbaseprc_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_dbase_prec), 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_dbaseprcaply + 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_dbaseprc_aply end interface - name='psb_dmlprcaply' + name='psb_dmlprc_aply' info = 0 call psb_erractionsave(err_act) @@ -843,7 +843,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) select case(baseprecv(2)%iprcparm(ml_type_)) case(no_ml_) ! Should not really get here. - write(0,*) 'Smooth preconditioner with no multilevel in MLPRCAPLY????' + write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????' case(add_ml_prec_) @@ -853,7 +853,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) 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_dbaseprcaply(baseprecv(1),x,beta,y,desc_data,trans,work,info) + 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_) @@ -912,7 +912,8 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) endif w2l=t2l - call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + & 'N',work,info) if (ismth /= no_smth_) then @@ -1015,7 +1016,8 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) t6 = mpi_wtime() w2l=t2l - call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + &'N',work,info) if(info /=0) goto 9999 if (ismth /= no_smth_) then @@ -1039,7 +1041,8 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) call psb_spmm(-one,baseprecv(2)%aorig,ty,one,tx,desc_data,info,work=work) if(info /=0) goto 9999 - call psb_dbaseprcaply(baseprecv(1),tx,one,ty,desc_data,trans,work,info) + call psb_baseprc_aply(baseprecv(1),tx,one,ty,desc_data,trans,& + & work,info) if(info /=0) goto 9999 call psb_axpby(one,ty,beta,y,desc_data,info) @@ -1073,7 +1076,8 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) call psb_axpby(one,y,zero,ty,desc_data,info) if(info /=0) goto 9999 - call psb_dbaseprcaply(baseprecv(1),x,zero,tty,desc_data,trans,work,info) + call psb_baseprc_aply(baseprecv(1),x,zero,tty,desc_data,& + & trans,work,info) if(info /=0) goto 9999 call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) @@ -1115,7 +1119,8 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) t6 = mpi_wtime() w2l=t2l - call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + & 'N',work,info) if(info /=0) goto 9999 if (ismth /= no_smth_) then @@ -1170,7 +1175,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) call psb_axpby(one,y,zero,ty,desc_data,info) if(info /=0) goto 9999 - call psb_dbaseprcaply(baseprecv(1),tx,zero,tty,desc_data,trans,work,info) + call psb_baseprc_aply(baseprecv(1),tx,zero,tty,desc_data,trans,work,info) if(info /=0) goto 9999 call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) @@ -1206,7 +1211,8 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) t6 = mpi_wtime() w2l=t2l - call psb_dbaseprcaply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,'N',work,info) + call psb_baseprc_aply(baseprecv(2),w2l,zero,t2l,baseprecv(2)%desc_data,& + & 'N',work,info) if(info /=0) goto 9999 if (ismth /= no_smth_) then @@ -1232,7 +1238,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) call psb_spmm(-one,baseprecv(2)%aorig,tty,one,tx,desc_data,info,work=work) if(info /=0) goto 9999 - call psb_dbaseprcaply(baseprecv(1),tx,one,tty,desc_data,'N',work,info) + call psb_baseprc_aply(baseprecv(1),tx,one,tty,desc_data,'N',work,info) call psb_axpby(one,tty,beta,y,desc_data,info) @@ -1247,7 +1253,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) end select case default - write(0,*) me, 'Wrong mltype into PRCAPLY ',& + write(0,*) me, 'Wrong mltype into PRC_APLY ',& & baseprecv(2)%iprcparm(ml_type_) end select @@ -1263,7 +1269,7 @@ subroutine psb_dmlprcaply(baseprecv,x,beta,y,desc_data,trans,work,info) end if return -end subroutine psb_dmlprcaply +end subroutine psb_dmlprc_aply !!$ !!$ @@ -1301,7 +1307,7 @@ end subroutine psb_dmlprcaply !!$ !!$ -subroutine psb_dprecaply1(prec,x,desc_data,info,trans) +subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) use psb_serial_mod use psb_descriptor_type @@ -1320,7 +1326,7 @@ subroutine psb_dprecaply1(prec,x,desc_data,info,trans) real(kind(1.d0)), parameter :: one=1.d0, zero=0.d0 interface - subroutine psb_dprecaply(prec,x,y,desc_data,info,trans, work) + subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) use psb_descriptor_type use psb_prec_type @@ -1332,7 +1338,7 @@ subroutine psb_dprecaply1(prec,x,desc_data,info,trans) integer, intent(out) :: info character(len=1), optional :: trans real(kind(0.d0)), optional, target :: work(:) - end subroutine psb_dprecaply + end subroutine psb_dprc_aply end interface ! Local variables @@ -1358,8 +1364,8 @@ subroutine psb_dprecaply1(prec,x,desc_data,info,trans) call psb_errpush(4010,name,a_err='Allocate') goto 9999 end if - if (debug) write(0,*) 'Prcaply1 Size(x) ',size(x), size(ww),size(w1) - call psb_dprecaply(prec,x,ww,desc_data,info,trans_,work=w1) + if (debug) write(0,*) 'Prc_aply1 Size(x) ',size(x), size(ww),size(w1) + call psb_dprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) if(info /=0) goto 9999 x(:) = ww(:) deallocate(ww,W1) @@ -1375,4 +1381,4 @@ subroutine psb_dprecaply1(prec,x,desc_data,info,trans) return end if return -end subroutine psb_dprecaply1 +end subroutine psb_dprc_aply1 diff --git a/src/prec/psb_dprecbld.f90 b/src/prec/psb_dprecbld.f90 index 4e8562db..2afa942f 100644 --- a/src/prec/psb_dprecbld.f90 +++ b/src/prec/psb_dprecbld.f90 @@ -52,47 +52,17 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) integer, intent(out) :: info character, intent(in), optional :: upd - interface psb_ilu_bld - subroutine psb_dilu_bld(a,desc_data,p,upd,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - integer, intent(out) :: info - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbase_prec), intent(inout) :: p - character, intent(in) :: upd - end subroutine psb_dilu_bld - end interface - - interface psb_slu_bld - subroutine psb_dslu_bld(a,desc_a,p,info) - use psb_serial_mod - use psb_descriptor_type - use psb_prec_type - use psb_const_mod - implicit none - - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info - end subroutine psb_dslu_bld - end interface - - interface psb_umf_bld - subroutine psb_dumf_bld(a,desc_a,p,info) - use psb_serial_mod + 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 - use psb_const_mod - implicit none - - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info - end subroutine psb_dumf_bld + type(psb_dspmat_type), target :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type),intent(inout) :: p + integer, intent(out) :: info + character, intent(in), optional :: upd + end subroutine psb_dbaseprc_bld end interface interface psb_mlprc_bld @@ -103,16 +73,16 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) use psb_const_mod implicit none - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type), intent(inout) :: p + integer, intent(out) :: info end subroutine psb_dmlprc_bld end interface ! Local scalars - Integer :: err, nnzero, n_row, n_col,I,j,k,icontxt,& - & me,mycol,nprow,npcol,mglob,lw, mtype, nrg, nzg, err_act + Integer :: err, nnzero, I,j,k,icontxt,& + & me,mycol,nprow,npcol,lw, mtype, nrg, nzg, err_act real(kind(1.d0)) :: temp, real_err(5) real(kind(1.d0)),pointer :: gd(:), work(:) integer :: int_err(5) @@ -132,25 +102,27 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) info = 0 int_err(1) = 0 icontxt = desc_a%matrix_data(psb_ctxt_) - n_row = desc_a%matrix_data(psb_n_row_) - n_col = desc_a%matrix_data(psb_n_col_) - mglob = desc_a%matrix_data(psb_m_) + if (debug) write(0,*) 'Preconditioner Blacs_gridinfo' call blacs_gridinfo(icontxt, nprow, npcol, me, mycol) if (present(upd)) then if (debug) write(0,*) 'UPD ', upd - if ((UPD.eq.'F').or.(UPD.eq.'T')) then - IUPD=UPD + if ((upd.eq.'F').or.(upd.eq.'T')) then + iupd=upd else - IUPD='F' + iupd='F' endif else - IUPD='F' + iupd='F' endif if (.not.associated(p%baseprecv)) then !! Error 1: should call precset + info=4010 + ch_err='unassociated bpv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 end if ! ! Should add check to ensure all procs have the same... @@ -158,202 +130,61 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) ! ALso should define symbolic names for the preconditioners. ! - call psb_check_def(p%baseprecv(1)%iprcparm(p_type_),'base_prec',& - & diagsc_,is_legal_base_prec) - allocate(p%baseprecv(1)%desc_data,stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - call psb_nullify_desc(p%baseprecv(1)%desc_data) - - select case(p%baseprecv(1)%iprcparm(p_type_)) - case (noprec_) - ! Do nothing. - - - case (diagsc_) - - if (debug) write(0,*) 'Precond: Diagonal scaling' - ! diagonal scaling - - call psb_realloc(n_col,p%baseprecv(1)%d,info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='psb_realloc') - goto 9999 - end if - - call psb_csrws(p%baseprecv(1)%d,a,info,trans='N') - if(info /= 0) then + if (size(p%baseprecv) >= 1) then + call init_baseprc_av(p%baseprecv(1),info) + if (info /= 0) then info=4010 - ch_err='psb_csrws' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 - end if - - if (debug) write(ilout+me,*) 'VDIAG ',n_row - do i=1,n_row - if (p%baseprecv(1)%d(i).eq.0.0d0) then - p%baseprecv(1)%d(i)=1.d0 - else - p%baseprecv(1)%d(i) = 1.d0/p%baseprecv(1)%d(i) - endif - - if (debug) write(ilout+me,*) i,desc_a%loc_to_glob(i),& - & p%baseprecv(1)%d(i) - if (p%baseprecv(1)%d(i).lt.0.d0) then - write(0,*) me,'Negative RWS? ',i,p%baseprecv(1)%d(i) - endif - end do - if (a%pl(1) /= 0) then - allocate(work(n_row),stat=info) - if (info /= 0) then - info=4000 - call psb_errpush(info,name) - goto 9999 - end if - call psb_gelp('n',a%pl,p%baseprecv(1)%d,desc_a,info) - if(info /= 0) then - info=4010 - ch_err='psb_dgelp' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(work) - endif - - if (debug) then - allocate(gd(mglob),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - call psb_gather(gd, p%baseprecv(1)%d, desc_a, info, iroot=iroot) - if(info /= 0) then - info=4010 - ch_err='psb_dgatherm' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - if (me.eq.iroot) then - write(iout+nprow,*) 'VDIAG CHECK ',mglob - do i=1,mglob - write(iout+nprow,*) i,gd(i) - enddo - endif - deallocate(gd) endif - if (debug) write(*,*) 'Preconditioner DIAG computed OK' - - - case (bja_,asm_) - - call psb_check_def(p%baseprecv(1)%iprcparm(n_ovr_),'overlap',& - & 0,is_legal_n_ovr) - call psb_check_def(p%baseprecv(1)%iprcparm(restr_),'restriction',& - & psb_halo_,is_legal_restrict) - call psb_check_def(p%baseprecv(1)%iprcparm(prol_),'prolongator',& - & psb_none_,is_legal_prolong) - call psb_check_def(p%baseprecv(1)%iprcparm(iren_),'renumbering',& - & renum_none_,is_legal_renum) - - if (debug) write(0,*)me, ': Calling PSB_ILU_BLD' - - -!!$ allocate(p%baseprecv(1)%av(bp_ilu_avsz),stat=info) - allocate(p%baseprecv(1)%av(max_avsz),stat=info) - do k=1,size(p%baseprecv(1)%av) - call psb_nullify_sp(p%baseprecv(1)%av(k)) - end do - - - select case(p%baseprecv(1)%iprcparm(f_type_)) - - case(f_ilu_n_,f_ilu_e_) - call psb_ilu_bld(a,desc_a,p%baseprecv(1),iupd,info) - if(debug) write(0,*)me,': out of psb_ilu_bld' - if(info /= 0) then - info=4010 - ch_err='psb_ilu_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_slu_) - - if(debug) write(0,*)me,': calling slu_bld' - call psb_slu_bld(a,desc_a,p%baseprecv(1),info) - if(info /= 0) then - info=4010 - ch_err='slu_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_umf_) - if(debug) write(0,*)me,': calling umf_bld' - call psb_umf_bld(a,desc_a,p%baseprecv(1),info) - if(info /= 0) then - info=4010 - ch_err='umf_bld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - case(f_none_) - write(0,*) 'Fact=None in PRECBLD Bja/ASM??' - - case default - write(0,*) 'Unknown factor type in precbld bja/asm: ',& - &p%baseprecv(1)%iprcparm(f_type_) - end select - - end select + + call psb_baseprc_bld(a,desc_a,p%baseprecv(1),info,iupd) + else + info=4010 + ch_err='size bpv' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 - if (size(p%baseprecv) >1) then + endif - if (.not.associated(p%baseprecv(2)%iprcparm)) then - info = 2222 - call psb_errpush(info,name) + 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_check_def(p%baseprecv(2)%iprcparm(ml_type_),'Multilevel type',& - & mult_ml_prec_,is_legal_ml_type) - call psb_check_def(p%baseprecv(2)%iprcparm(aggr_alg_),'aggregation',& - & loc_aggr_,is_legal_ml_aggr_kind) - call psb_check_def(p%baseprecv(2)%iprcparm(smth_kind_),'Smoother kind',& - & smth_omg_,is_legal_ml_smth_kind) - call psb_check_def(p%baseprecv(2)%iprcparm(coarse_mat_),'Coarse matrix',& - & mat_distr_,is_legal_ml_coarse_mat) - call psb_check_def(p%baseprecv(2)%iprcparm(smth_pos_),'smooth_pos',& - & pre_smooth_,is_legal_ml_smooth_pos) - call psb_check_def(p%baseprecv(2)%iprcparm(f_type_),'fact',f_ilu_n_,is_legal_ml_fact) - nullify(p%baseprecv(2)%desc_data) - select case(p%baseprecv(2)%iprcparm(f_type_)) - case(f_ilu_n_) - call psb_check_def(p%baseprecv(2)%iprcparm(ilu_fill_in_),'Level',0,is_legal_ml_lev) - case(f_ilu_e_) - call psb_check_def(p%baseprecv(2)%dprcparm(fact_eps_),'Eps',0.0d0,is_legal_ml_eps) - end select - call psb_check_def(p%baseprecv(2)%dprcparm(smooth_omega_),'omega',0.0d0,is_legal_omega) - call psb_check_def(p%baseprecv(2)%iprcparm(jac_sweeps_),'Jacobi sweeps',& - & 1,is_legal_jac_sweeps) - 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) + call init_baseprc_av(p%baseprecv(i),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(p%baseprecv(i-1)%av(ac_),p%baseprecv(i-1)%desc_data,& + & p%baseprecv(i),info) + end do + endif call psb_erractionrestore(err_act) @@ -367,5 +198,20 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) end if return +contains + + subroutine init_baseprc_av(p,info) + type(psb_dbaseprc_type), intent(inout) :: p + integer :: info + if (associated(p%av)) then + ! Have not decided what to do yet + end if + allocate(p%av(max_avsz),stat=info) + if (info /= 0) return + do k=1,size(p%av) + call psb_nullify_sp(p%av(k)) + end do + end subroutine init_baseprc_av + end subroutine psb_dprecbld diff --git a/src/prec/psb_dprecset.f90 b/src/prec/psb_dprecset.f90 index 9da2bece..a9cfd3a9 100644 --- a/src/prec/psb_dprecset.f90 +++ b/src/prec/psb_dprecset.f90 @@ -48,7 +48,7 @@ subroutine psb_dprecset(p,ptype,iv,rs,rv,info) real(kind(1.d0)), optional, intent(in) :: rv(:) integer, optional, intent(out) :: info - type(psb_dbase_prec), pointer :: bpv(:)=>null() + type(psb_dbaseprc_type), pointer :: bpv(:)=>null() character(len=len(ptype)) :: typeup integer :: isz, err diff --git a/src/prec/psb_dslu_bld.f90 b/src/prec/psb_dslu_bld.f90 index 5f649d02..9b422325 100644 --- a/src/prec/psb_dslu_bld.f90 +++ b/src/prec/psb_dslu_bld.f90 @@ -41,10 +41,10 @@ subroutine psb_dslu_bld(a,desc_a,p,info) use psb_const_mod implicit none - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type), intent(inout) :: p + integer, intent(out) :: info type(psb_dspmat_type) :: blck, atmp diff --git a/src/prec/psb_dsp_renum.f90 b/src/prec/psb_dsp_renum.f90 index eaec5598..eec817c9 100644 --- a/src/prec/psb_dsp_renum.f90 +++ b/src/prec/psb_dsp_renum.f90 @@ -45,10 +45,10 @@ subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info) implicit none ! .. array Arguments .. - type(psb_dspmat_type), intent(in) :: a,blck - type(psb_dspmat_type), intent(inout) :: atmp - type(psb_dbase_prec), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a + type(psb_dspmat_type), intent(in) :: a,blck + type(psb_dspmat_type), intent(inout) :: atmp + type(psb_dbaseprc_type), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info diff --git a/src/prec/psb_dumf_bld.f90 b/src/prec/psb_dumf_bld.f90 index 5b80586f..d9b8e98e 100644 --- a/src/prec/psb_dumf_bld.f90 +++ b/src/prec/psb_dumf_bld.f90 @@ -42,10 +42,10 @@ subroutine psb_dumf_bld(a,desc_a,p,info) use psb_const_mod implicit none - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_dbase_prec), intent(inout) :: p - integer, intent(out) :: info + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + type(psb_dbaseprc_type), intent(inout) :: p + integer, intent(out) :: info type(psb_dspmat_type) :: blck, atmp