From 54ae8fa06aaae743c030ae5788f492a2a7753ebe Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 1 Jan 2008 19:11:07 +0000 Subject: [PATCH] Reorganized AS build and application routines. --- Make.inc | 2 +- mlprec/Makefile | 8 +- mlprec/mld_das_aply.f90 | 295 +++++++++++++++++ .../{mld_dasmat_bld.f90 => mld_das_bld.f90} | 159 ++++----- mlprec/mld_dbaseprec_aply.f90 | 306 +++++++++--------- mlprec/mld_dbaseprec_bld.f90 | 31 +- mlprec/mld_dbjac_bld.f90 | 106 +++--- mlprec/mld_dilu_bld.f90 | 11 +- mlprec/mld_dprec_aply.f90 | 8 +- mlprec/mld_dslud_bld.f90 | 4 +- mlprec/mld_dsp_renum.f90 | 17 +- mlprec/mld_prec_mod.f90 | 195 ++++++----- mlprec/mld_zas_aply.f90 | 295 +++++++++++++++++ .../{mld_zasmat_bld.f90 => mld_zas_bld.f90} | 165 ++++------ mlprec/mld_zbaseprec_aply.f90 | 45 ++- mlprec/mld_zbaseprec_bld.f90 | 31 +- mlprec/mld_zbjac_bld.f90 | 97 ++++-- mlprec/mld_zilu_bld.f90 | 11 +- mlprec/mld_zprec_aply.f90 | 8 +- mlprec/mld_zslud_bld.f90 | 4 +- mlprec/mld_zsp_renum.f90 | 17 +- test/fileread/runs/drt.sh | 10 +- 22 files changed, 1214 insertions(+), 611 deletions(-) create mode 100644 mlprec/mld_das_aply.f90 rename mlprec/{mld_dasmat_bld.f90 => mld_das_bld.f90} (59%) create mode 100644 mlprec/mld_zas_aply.f90 rename mlprec/{mld_zasmat_bld.f90 => mld_zas_bld.f90} (56%) diff --git a/Make.inc b/Make.inc index 02353068..2cc6be90 100644 --- a/Make.inc +++ b/Make.inc @@ -4,7 +4,7 @@ # must be specified here with absolute pathnames # # # ########################################################## -PSBLASDIR=$(HOME)/NUMERICAL/PSBLAS2/psblas2 +PSBLASDIR=$(HOME)/NUMERICAL/PSBLAS2/psblas2-dev include $(PSBLASDIR)/Make.inc ########################################################## diff --git a/mlprec/Makefile b/mlprec/Makefile index b6ad53ef..7e9857c6 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -10,18 +10,18 @@ MODOBJS=mld_prec_type.o mld_prec_mod.o MPFOBJS=mld_daggrmat_raw_asb.o mld_daggrmat_smth_asb.o \ mld_zaggrmat_raw_asb.o mld_zaggrmat_smth_asb.o MPCOBJS=mld_slud_impl.o mld_zslud_impl.o -F90OBJS=mld_dasmat_bld.o mld_dslu_bld.o mld_dumf_bld.o mld_dilu_fct.o\ +F90OBJS=mld_das_bld.o mld_dslu_bld.o mld_dumf_bld.o mld_dilu_fct.o\ mld_dmlprec_bld.o mld_dsp_renum.o mld_dbjac_bld.o mld_dilu_bld.o \ mld_dprecbld.o mld_dprecfree.o mld_dprecset.o mld_dprecinit.o\ mld_dbaseprec_bld.o mld_ddiagsc_bld.o mld_daggrmap_bld.o \ mld_dprec_aply.o mld_dmlprec_aply.o mld_dslud_bld.o\ - mld_dbaseprec_aply.o mld_dbjac_aply.o mld_daggrmat_asb.o \ - mld_zasmat_bld.o mld_zslu_bld.o mld_zumf_bld.o mld_zilu_fct.o\ + mld_dbaseprec_aply.o mld_dbjac_aply.o mld_das_aply.o mld_daggrmat_asb.o \ + mld_zas_bld.o mld_zslu_bld.o mld_zumf_bld.o mld_zilu_fct.o\ mld_zmlprec_bld.o mld_zsp_renum.o mld_zbjac_bld.o mld_zilu_bld.o \ mld_zprecbld.o mld_zprecfree.o mld_zprecset.o mld_zprecinit.o \ mld_zbaseprec_bld.o mld_zdiagsc_bld.o mld_zaggrmap_bld.o \ mld_zprec_aply.o mld_zmlprec_aply.o mld_zslud_bld.o\ - mld_zbaseprec_aply.o mld_zbjac_aply.o mld_zaggrmat_asb.o\ + mld_zbaseprec_aply.o mld_zbjac_aply.o mld_zas_aply.o mld_zaggrmat_asb.o\ mld_diluk_fct.o mld_ziluk_fct.o mld_dilut_fct.o mld_zilut_fct.o \ $(MPFOBJS) COBJS=mld_slu_impl.o mld_umf_impl.o mld_zslu_impl.o mld_zumf_impl.o diff --git a/mlprec/mld_das_aply.f90 b/mlprec/mld_das_aply.f90 new file mode 100644 index 00000000..da5f0c80 --- /dev/null +++ b/mlprec/mld_das_aply.f90 @@ -0,0 +1,295 @@ +!!$ +!!$ +!!$ MLD2P4 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS v.2.0) +!!$ +!!$ (C) Copyright 2007 Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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. +!!$ +!!$ +! File: mld_das_aply.f90 +! +! Subroutine: mld_das_aply +! Version: real +! +! This routine applies the Additive Schwarz preconditioner by computing +! +! Y = beta*Y + alpha*op(K^(-1))*X, +! where +! - K is the base preconditioner, stored in prec, +! - op(K^(-1)) is K^(-1) or its transpose, according to the value of trans, +! - X and Y are vectors, +! - alpha and beta are scalars. +! +! +! Arguments: +! alpha - real(kind(0.d0)), input. +! The scalar alpha. +! prec - type(mld_dbaseprc_type), input. +! The base preconditioner data structure containing the local part +! of the preconditioner K. +! x - real(kind(0.d0)), dimension(:), input. +! The local part of the vector X. +! beta - real(kind(0.d0)), input. +! The scalar beta. +! y - real(kind(0.d0)), dimension(:), input/output. +! The local part of the vector Y. +! desc_data - type(psb_desc_type), input. +! The communication descriptor associated to the matrix to be +! preconditioned. +! trans - character, optional. +! If trans='N','n' then op(K^(-1)) = K^(-1); +! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). +! work - real(kind(0.d0)), dimension (:), optional, target. +! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). +! info - integer, output. +! Error code. +! +subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_das_aply + + implicit none + + ! Arguments + type(psb_desc_type),intent(in) :: desc_data + type(mld_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(in) :: x(:) + real(kind(0.d0)),intent(inout) :: y(:) + real(kind(0.d0)),intent(in) :: alpha,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), nrow_d + real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,isz, err_act + character(len=20) :: name, ch_err + + name='mld_das_aply' + info = 0 + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_data) + + call psb_info(ictxt, me, np) + + select case(trans) + case('N','n') + case('T','t','C','c') + case default + info=40 + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + select case(prec%iprcparm(mld_prec_type_)) + + case(mld_as_) + ! + ! Additive Schwarz preconditioner + ! + + if (prec%iprcparm(mld_n_ovr_)==0) then + ! + ! shortcut: this fixes performance for RAS(0) == BJA + ! + call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + if(info /= 0) then + info=4010 + ch_err='psb_bjacaply' + goto 9999 + end if + + else + + ! + ! Note: currently trans is unused + ! + + n_row = psb_cd_get_local_rows(prec%desc_data) + n_col = psb_cd_get_local_cols(prec%desc_data) + nrow_d = psb_cd_get_local_rows(desc_data) + isz=max(n_row,N_COL) + if ((6*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + aux => work(3*isz+1:) + else if ((4*isz) <= size(work)) then + aux => work(1:) + allocate(ww(isz),tx(isz),ty(isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/3*isz,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + else if ((3*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + allocate(aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + else + allocate(ww(isz),tx(isz),ty(isz),& + &aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + + endif + + tx(1:nrow_d) = x(1:nrow_d) + tx(nrow_d+1:isz) = dzero + + ! + ! Get the overlap entries of tx (tx==x) + ! + if (prec%iprcparm(mld_sub_restr_)==psb_halo_) then + call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_) + if(info /=0) then + info=4010 + ch_err='psb_halo' + goto 9999 + end if + else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then + call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + ! + ! If required, reorder tx according to the row/column permutation of the + ! local extended matrix, stored into the permutation vector prec%perm + ! + if (prec%iprcparm(mld_sub_ren_)>0) then + call psb_gelp('n',prec%perm,tx,info) + if(info /=0) then + info=4010 + ch_err='psb_gelp' + goto 9999 + end if + endif + + ! + ! Apply to tx the block-Jacobi preconditioner/solver (multiple sweeps of the + ! block-Jacobi solver can be applied at the coarsest level of a multilevel + ! preconditioner). The resulting vector is ty. + ! + call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info) + if(info /= 0) then + info=4010 + ch_err='mld_bjac_aply' + goto 9999 + end if + + ! + ! Apply to ty the inverse permutation of prec%perm + ! + if (prec%iprcparm(mld_sub_ren_)>0) then + call psb_gelp('n',prec%invperm,ty,info) + if(info /= 0) then + info=4010 + ch_err='psb_gelp' + goto 9999 + end if + endif + + select case (prec%iprcparm(mld_sub_prol_)) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(ty,prec%desc_data,info,& + ! & update=prec%iprcparm(mld_sub_prol_),work=aux) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of ty + ! + call psb_ovrl(ty,prec%desc_data,info,& + & update=prec%iprcparm(mld_sub_prol_),work=aux) + if(info /=0) then + info=4010 + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + ! + ! Compute y = beta*y + alpha*ty (ty==K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + if ((6*isz) <= size(work)) then + else if ((4*isz) <= size(work)) then + deallocate(ww,tx,ty) + else if ((3*isz) <= size(work)) then + deallocate(aux) + else + deallocate(ww,aux,tx,ty) + endif + end if + + case default + + call psb_errpush(4001,name,a_err='Invalid mld_prec_type_') + goto 9999 + + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_das_aply + diff --git a/mlprec/mld_dasmat_bld.f90 b/mlprec/mld_das_bld.f90 similarity index 59% rename from mlprec/mld_dasmat_bld.f90 rename to mlprec/mld_das_bld.f90 index ee4542cd..906378aa 100644 --- a/mlprec/mld_dasmat_bld.f90 +++ b/mlprec/mld_das_bld.f90 @@ -34,15 +34,13 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: mld_dasmat_bld.f90 +! File: mld_das_bld.f90 ! -! Subroutine: mld_dasmat_bld +! Subroutine: mld_das_bld ! Version: real ! -! This routine builds the communication descriptor associated to the extended -! matrices that form the Additive Schwarz (AS) preconditioner and retrieves -! the remote pieces needed to build the local extended matrix. If the -! preconditioner is the block-Jacobi one, the routine makes only a copy of +! This routine builds the Additive Schwarz (AS) preconditioner. +! If the preconditioner is the block-Jacobi one, the routine makes only a copy of ! the descriptor of the original matrix. ! ! @@ -71,32 +69,32 @@ ! preconditioner. Currently outfmt is set to 'CSR' by the ! calling routine mld_bjac_bld. ! -subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) +subroutine mld_das_bld(a,desc_a,p,upd,info) use psb_base_mod - use mld_prec_mod, mld_protect_name => mld_dasmat_bld - + use mld_prec_mod, mld_protect_name => mld_das_bld + Implicit None ! Arguments - integer, intent(in) :: ptype,novr - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_dspmat_type), Intent(out) :: blk - integer, intent(out) :: info - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - character(len=5), optional :: outfmt + ! Arguments + type(psb_dspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + type(mld_dbaseprc_type), intent(inout) :: p + character, intent(in) :: upd + integer, intent(out) :: info + + integer :: ptype,novr ! Local variables integer icomm Integer :: np,me,nnzero,& & ictxt, n_col,int_err(5),& - & tot_recv, n_row,nhalo, nrow_a,err_act + & tot_recv, n_row,nhalo, nrow_a,err_act, data_ integer :: debug_level, debug_unit character(len=20) :: name, ch_err - name='mld_dasmat_bld' + name='mld_as_bld' if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) @@ -106,45 +104,21 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) If (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & ' start ', upd - ictxt = psb_cd_get_context(desc_data) - icomm = psb_cd_get_mpic(desc_data) + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) Call psb_info(ictxt, me, np) tot_recv=0 - nrow_a = desc_data%matrix_data(psb_n_row_) - nnzero = Size(a%aspk) - n_col = desc_data%matrix_data(psb_n_col_) + nrow_a = psb_cd_get_local_rows(desc_a) + n_col = psb_cd_get_local_cols(desc_a) + nnzero = psb_sp_get_nnzeros(a) nhalo = n_col-nrow_a + ptype = p%iprcparm(mld_prec_type_) + novr = p%iprcparm(mld_n_ovr_) select case (ptype) - case (mld_bjac_) - ! - ! Block-Jacobi preconditioner. Copy the descriptor, just in case - ! we want to renumber the rows and columns of the matrix. - ! - call psb_sp_all(0,0,blk,1,info) - if(info /= 0) then - info=4010 - ch_err='psb_sp_all' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - blk%fida = 'COO' - blk%infoa(psb_nnz_) = 0 - If (upd == 'F') Then - call psb_cdcpy(desc_data,desc_p,info) - If(debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' done cdcpy' - if(info /= 0) then - info=4010 - ch_err='psb_cdcpy' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif case(mld_as_) ! @@ -161,23 +135,12 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) ! ! Actually, this is just block Jacobi ! - If(debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' calling allocate novr=0' - call psb_sp_all(0,0,blk,1,info) - if(info /= 0) then - info=4010 - ch_err='psb_sp_all' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - blk%fida = 'COO' - blk%infoa(psb_nnz_) = 0 + data_ = psb_no_comm_ if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Calling desccpy' if (upd == 'F') then - call psb_cdcpy(desc_data,desc_p,info) + call psb_cdcpy(desc_a,p%desc_data,info) If(debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & ' done cdcpy' @@ -191,63 +154,49 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) & write(debug_unit,*) me,' ',trim(name),& & 'Early return: P>=3 N_OVR=0' endif - return - endif + else + data_ = psb_comm_ext_ + If (upd == 'F') Then + ! + ! Build the auxiliary descriptor desc_p%matrix_data(psb_n_row_). + ! This is done by psb_cdbldext (interface to psb_cdovr), which is + ! independent of CSR, and has been placed in the tools directory + ! of PSBLAS, instead of the mlprec directory of MLD2P4, because it + ! might be used independently of the AS preconditioner, to build + ! a descriptor for an extended stencil in a PDE solver. + ! + call psb_cdbldext(a,desc_a,novr,p%desc_data,info,extype=psb_ovt_asov_) + if (info /= 0) then + info=4010 + ch_err='psb_cdbldext' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + Endif - If (upd == 'F') Then - ! - ! Build the auxiliary descriptor desc_p%matrix_data(psb_n_row_). - ! This is done by psb_cdbldext (interface to psb_cdovr), which is - ! independent of CSR, and has been placed in the tools directory - ! of PSBLAS, instead of the mlprec directory of MLD2P4, because it - ! might be used independently of the AS preconditioner, to build - ! a descriptor for an extended stencil in a PDE solver. - ! - call psb_cdbldext(a,desc_data,novr,desc_p,info,extype=psb_ovt_asov_) - if (info /= 0) then - info=4010 - ch_err='psb_cdbldext' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - Endif + End if if(debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' From cdbldext _:',desc_p%matrix_data(psb_n_row_),& - & desc_p%matrix_data(psb_n_col_) - - ! - ! Retrieve the remote pieces needed to build the local extended matrix - ! - - n_row = desc_p%matrix_data(psb_n_row_) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Before sphalo ',blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_) - Call psb_sphalo(a,desc_p,blk,info,& - & outfmt=outfmt,data=psb_comm_ext_,rowscale=.true.) - + & ' From cdbldext _:',p%desc_data%matrix_data(psb_n_row_),& + & p%desc_data%matrix_data(psb_n_col_) + + call mld_bjac_bld(a,p,upd,info,data=data_) if (info /= 0) then info=4010 - ch_err='psb_sphalo' + ch_err='mld_bjac_bld' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'After psb_sphalo ',& - & blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_) - case default + info=4001 ch_err='Invalid ptype' call psb_errpush(info,name,a_err=ch_err) goto 9999 - + End select if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),'Done' @@ -263,5 +212,5 @@ subroutine mld_dasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) end if Return -End Subroutine mld_dasmat_bld +End Subroutine mld_das_bld diff --git a/mlprec/mld_dbaseprec_aply.f90 b/mlprec/mld_dbaseprec_aply.f90 index 01993841..ab3fca23 100644 --- a/mlprec/mld_dbaseprec_aply.f90 +++ b/mlprec/mld_dbaseprec_aply.f90 @@ -83,7 +83,7 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) implicit none -! Arguments + ! Arguments type(psb_desc_type),intent(in) :: desc_data type(mld_dbaseprc_type), intent(in) :: prec real(kind(0.d0)),intent(in) :: x(:) @@ -96,9 +96,8 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! Local variables integer :: n_row,n_col, int_err(5), nrow_d real(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) - character ::diagl, diagu - integer :: ictxt,np,me,isz, err_act - character(len=20) :: name, ch_err + integer :: ictxt,np,me,isz, err_act + character(len=20) :: name, ch_err name='mld_dbaseprec_aply' info = 0 @@ -108,9 +107,6 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) call psb_info(ictxt, me, np) - diagl='U' - diagu='U' - select case(trans) case('N','n') case('T','t','C','c') @@ -124,17 +120,17 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) select case(prec%iprcparm(mld_prec_type_)) case(mld_noprec_) - ! - ! No preconditioner - ! + ! + ! No preconditioner + ! call psb_geaxpby(alpha,x,beta,y,desc_data,info) case(mld_diag_) - ! - ! Diagonal preconditioner - ! - + ! + ! Diagonal preconditioner + ! + if (size(work) >= size(x)) then ww => work else @@ -158,9 +154,9 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) end if case(mld_bjac_) - ! - ! Block-Jacobi preconditioner - ! + ! + ! Block-Jacobi preconditioner + ! call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info /= 0) then @@ -170,160 +166,170 @@ subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) end if case(mld_as_) - ! - ! Additive Schwarz preconditioner - ! - - if (prec%iprcparm(mld_n_ovr_)==0) then - ! - ! shortcut: this fixes performance for RAS(0) == BJA - ! - call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - if(info /= 0) then - info=4010 - ch_err='psb_bjacaply' - goto 9999 - end if - - else + ! + ! Additive Schwarz preconditioner + ! + call mld_as_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + if(info /= 0) then + info=4010 + ch_err='mld_as_aply' + goto 9999 + end if - ! - ! Note: currently trans is unused - ! - - n_row = psb_cd_get_local_rows(prec%desc_data) - n_col = psb_cd_get_local_cols(prec%desc_data) - nrow_d = psb_cd_get_local_rows(desc_data) - isz=max(n_row,N_COL) - if ((6*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - aux => work(3*isz+1:) - else if ((4*isz) <= size(work)) then - aux => work(1:) - allocate(ww(isz),tx(isz),ty(isz),stat=info) - if (info /= 0) then - call psb_errpush(4025,name,i_err=(/3*isz,0,0,0,0/),& - & a_err='real(kind(1.d0))') - goto 9999 - end if - else if ((3*isz) <= size(work)) then - ww => work(1:isz) - tx => work(isz+1:2*isz) - ty => work(2*isz+1:3*isz) - allocate(aux(4*isz),stat=info) - if (info /= 0) then - call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& - & a_err='real(kind(1.d0))') - goto 9999 - end if - else - allocate(ww(isz),tx(isz),ty(isz),& - &aux(4*isz),stat=info) - if (info /= 0) then - call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& - & a_err='real(kind(1.d0))') - goto 9999 + if (.false.) then + if (prec%iprcparm(mld_n_ovr_)==0) then + ! + ! shortcut: this fixes performance for RAS(0) == BJA + ! + call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + if(info /= 0) then + info=4010 + ch_err='psb_bjacaply' + goto 9999 end if - endif + else - tx(1:nrow_d) = x(1:nrow_d) - tx(nrow_d+1:isz) = dzero + ! + ! Note: currently trans is unused + ! - ! - ! Get the overlap entries of tx (tx==x) - ! + n_row = psb_cd_get_local_rows(prec%desc_data) + n_col = psb_cd_get_local_cols(prec%desc_data) + nrow_d = psb_cd_get_local_rows(desc_data) + isz=max(n_row,N_COL) + if ((6*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + aux => work(3*isz+1:) + else if ((4*isz) <= size(work)) then + aux => work(1:) + allocate(ww(isz),tx(isz),ty(isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/3*isz,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + else if ((3*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + allocate(aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + else + allocate(ww(isz),tx(isz),ty(isz),& + &aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='real(kind(1.d0))') + goto 9999 + end if + + endif + + tx(1:nrow_d) = x(1:nrow_d) + tx(nrow_d+1:isz) = dzero + + ! + ! Get the overlap entries of tx (tx==x) + ! if (prec%iprcparm(mld_sub_restr_)==psb_halo_) then - call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_) - if(info /=0) then - info=4010 - ch_err='psb_halo' + call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_) + if(info /=0) then + info=4010 + ch_err='psb_halo' + goto 9999 + end if + else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then + call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_') goto 9999 end if - else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then - call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_') - goto 9999 - end if - ! - ! If required, reorder tx according to the row/column permutation of the - ! local extended matrix, stored into the permutation vector prec%perm - ! - if (prec%iprcparm(mld_sub_ren_)>0) then - call psb_gelp('n',prec%perm,tx,info) - if(info /=0) then - info=4010 - ch_err='psb_gelp' - goto 9999 - end if - endif - - ! - ! Apply to tx the block-Jacobi preconditioner/solver (multiple sweeps of the - ! block-Jacobi solver can be applied at the coarsest level of a multilevel - ! preconditioner). The resulting vector is ty. - ! - call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info) - if(info /= 0) then - info=4010 - ch_err='mld_bjac_aply' - goto 9999 - end if + ! + ! If required, reorder tx according to the row/column permutation of the + ! local extended matrix, stored into the permutation vector prec%perm + ! + if (prec%iprcparm(mld_sub_ren_)>0) then + call psb_gelp('n',prec%perm,tx,info) + if(info /=0) then + info=4010 + ch_err='psb_gelp' + goto 9999 + end if + endif - ! - ! Apply to ty the inverse permutation of prec%perm - ! - if (prec%iprcparm(mld_sub_ren_)>0) then - call psb_gelp('n',prec%invperm,ty,info) + ! + ! Apply to tx the block-Jacobi preconditioner/solver (multiple sweeps of the + ! block-Jacobi solver can be applied at the coarsest level of a multilevel + ! preconditioner). The resulting vector is ty. + ! + call mld_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info) if(info /= 0) then info=4010 - ch_err='psb_gelp' + ch_err='mld_bjac_aply' goto 9999 end if - endif - select case (prec%iprcparm(mld_sub_prol_)) - - case(psb_none_) - ! - ! Would work anyway, but since it is supposed to do nothing ... - ! call f90_psovrl(ty,prec%desc_data,update=prec%a_restrict) ! + ! Apply to ty the inverse permutation of prec%perm + ! + if (prec%iprcparm(mld_sub_ren_)>0) then + call psb_gelp('n',prec%invperm,ty,info) + if(info /= 0) then + info=4010 + ch_err='psb_gelp' + goto 9999 + end if + endif + + select case (prec%iprcparm(mld_sub_prol_)) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(ty,prec%desc_data,info,& + ! & update=prec%iprcparm(mld_sub_prol_),work=aux) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of ty + ! + call psb_ovrl(ty,prec%desc_data,info,& + & update=prec%iprcparm(mld_sub_prol_),work=aux) + if(info /=0) then + info=4010 + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_') + goto 9999 + end select - case(psb_sum_,psb_avg_) ! - ! Update the overlap of ty + ! Compute y = beta*y + alpha*ty (ty==K^(-1)*tx) ! - call psb_ovrl(ty,prec%desc_data,info,& - & update=prec%iprcparm(mld_sub_prol_),work=aux) - if(info /=0) then - info=4010 - ch_err='psb_ovrl' - goto 9999 - end if + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + - case default - call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_') - goto 9999 - end select - - ! - ! Compute y = beta*y + alpha*ty (ty==K^(-1)*tx) - ! - call psb_geaxpby(alpha,ty,beta,y,desc_data,info) - - - if ((6*isz) <= size(work)) then - else if ((4*isz) <= size(work)) then - deallocate(ww,tx,ty) - else if ((3*isz) <= size(work)) then - deallocate(aux) - else - deallocate(ww,aux,tx,ty) - endif + if ((6*isz) <= size(work)) then + else if ((4*isz) <= size(work)) then + deallocate(ww,tx,ty) + else if ((3*isz) <= size(work)) then + deallocate(aux) + else + deallocate(ww,aux,tx,ty) + endif + end if end if + case default call psb_errpush(4001,name,a_err='Invalid mld_prec_type_') diff --git a/mlprec/mld_dbaseprec_bld.f90 b/mlprec/mld_dbaseprec_bld.f90 index 1638b86e..247f022f 100644 --- a/mlprec/mld_dbaseprec_bld.f90 +++ b/mlprec/mld_dbaseprec_bld.f90 @@ -149,7 +149,31 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd) goto 9999 end if - case(mld_bjac_,mld_as_) + case(mld_bjac_) + + call mld_check_def(p%iprcparm(mld_sub_ren_),'renumbering',& + & mld_renum_none_,is_legal_renum) + call mld_check_def(p%iprcparm(mld_sub_solve_),'fact',& + & mld_ilu_n_,is_legal_ml_fact) + + 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 + + ! Build the local part of the base preconditioner + call mld_bjac_bld(a,p,iupd,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='mld_bjac_bld') + goto 9999 + end if + + + case(mld_as_) ! Block Jacobi and additive Schwarz preconditioners/smoothers call mld_check_def(p%iprcparm(mld_n_ovr_),'overlap',& @@ -165,7 +189,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd) ! Set parameters for using SuperLU_dist on the local submatrices if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then - p%iprcparm(mld_n_ovr_) = 0 + p%iprcparm(mld_n_ovr_) = 0 p%iprcparm(mld_smooth_sweeps_) = 1 end if @@ -174,7 +198,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd) & ': Calling mld_bjac_bld' ! Build the local part of the base preconditioner - call mld_bjac_bld(a,desc_a,p,iupd,info) + call mld_as_bld(a,desc_a,p,iupd,info) if(info /= 0) then info=4010 call psb_errpush(info,name,a_err='mld_bjac_bld') @@ -182,6 +206,7 @@ subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd) end if case default + info=4001 ch_err='Unknown mld_prec_type_' call psb_errpush(info,name,a_err=ch_err) diff --git a/mlprec/mld_dbjac_bld.f90 b/mlprec/mld_dbjac_bld.f90 index 80f2d53e..68473688 100644 --- a/mlprec/mld_dbjac_bld.f90 +++ b/mlprec/mld_dbjac_bld.f90 @@ -84,27 +84,32 @@ ! a - type(psb_dspmat_type), input. ! The sparse matrix structure containing the local part of the ! matrix to be preconditioned or factorized. -! desc_a - type(psb_desc_type), input. -! The communication descriptor associated to a. ! p - type(mld_dbaseprec_type), input/output. ! The 'base preconditioner' data structure containing the local ! part of the preconditioner or solver at the current level. +! ! info - integer, output. ! Error code. +! data - integer Which index list in desc_a should be used to retrieve +! rows, default psb_comm_halo_ +! psb_no_comm_ don't retrieve remote rows +! psb_comm_halo_ use halo_index +! psb_comm_ext_ use ext_index +! psb_comm_ovrl_ DISABLED for this routine. ! -subroutine mld_dbjac_bld(a,desc_a,p,upd,info) +subroutine mld_dbjac_bld(a,p,upd,info,data) use psb_base_mod use mld_prec_mod, mld_protect_name => mld_dbjac_bld implicit none - -! Arguments + + ! Arguments type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in) :: desc_a type(mld_dbaseprc_type), intent(inout) :: p integer, intent(out) :: info character, intent(in) :: upd + integer, intent(in), optional :: data ! Local Variables integer :: i, k, m @@ -112,9 +117,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) character :: trans, unitd type(psb_dspmat_type) :: blck, atmp integer :: debug_level, debug_unit - integer :: err_act, n_row, nrow_a,n_col + integer :: err_act, n_row, nrow_a,n_col, data_ integer :: ictxt,np,me - character(len=20) :: name + character(len=20) :: name, ch_err character(len=5), parameter :: coofmt='COO', csrfmt='CSR' if(psb_get_errstatus().ne.0) return @@ -123,7 +128,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(p%desc_data) call psb_info(ictxt, me, np) m = a%m @@ -151,25 +156,50 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),': Start',& & p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_) - - ! - ! Build the communication descriptor for the Additive Schwarz - ! preconditioner and retrieves the remote pieces of the local - ! extended matrix needed by that preconditioner. If the - ! preconditioner is the block-Jacobi one, only a copy of the - ! descriptor of the original matrix is made. - ! - call mld_asmat_bld(p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info,outfmt=csrfmt) - - if (info/=0) then - call psb_errpush(4010,name,a_err='mld_asmat_bld') - goto 9999 + if (present(data)) then + data_ = data + else + data_ = psb_no_comm_ end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ': out of mld_asmat_bld' + select case (data_) + case(psb_no_comm_) + If (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' calling allocate novr=0' + call psb_sp_all(0,0,blck,1,info) + if(info /= 0) then + info=4010 + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + blck%fida = 'COO' + blck%infoa(psb_nnz_) = 0 + + case(psb_comm_halo_,psb_comm_ext_) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Before sphalo ',blck%fida,blck%m,psb_nnz_,blck%infoa(psb_nnz_) + Call psb_sphalo(a,p%desc_data,blck,info,data=data_,rowscale=.true.) + + if (info /= 0) then + info=4010 + ch_err='psb_sphalo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'After psb_sphalo ',& + & blck%fida,blck%m,psb_nnz_,blck%infoa(psb_nnz_) + case default + info=35 + call psb_errpush(info,name,i_err=(/5,data,0,0,0/)) + goto 9999 + end select ! ! Treat separately the case the local matrix has to be reordered @@ -177,9 +207,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) ! select case(p%iprcparm(mld_sub_ren_)) - ! - ! A reordering of the local matrix is required. - ! + ! + ! A reordering of the local matrix is required. + ! case (1:) ! @@ -187,7 +217,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) ! according to the value of p%iprcparm(sub_ren_). The reordered ! matrix is stored into atmp, using the COO format. ! - call mld_sp_renum(a,desc_a,blck,p,atmp,info) + call mld_sp_renum(a,blck,p,atmp,info) if (info/=0) then call psb_errpush(4010,name,a_err='mld_sp_renum') goto 9999 @@ -233,7 +263,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) ! ILU(k)/MILU(k)/ILU(k,t) factorization. ! call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_) - if (info == 0) call mld_ilu_bld(atmp,p%desc_data,p,upd,info) + if (info == 0) call mld_ilu_bld(atmp,p,upd,info) if (info/=0) then call psb_errpush(4010,name,a_err='mld_ilu_bld') goto 9999 @@ -282,9 +312,9 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) goto 9999 end if - ! - ! No reordering of the local matrix is required - ! + ! + ! No reordering of the local matrix is required + ! case(0) ! @@ -321,7 +351,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) call psb_errpush(4010,name,a_err='clip & psb_spcnv csr 4') goto 9999 end if - + k = psb_sp_get_nnzeros(p%av(mld_ap_nd_)) call psb_sum(ictxt,k) @@ -338,7 +368,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) ! ! Compute the incomplete LU factorization. ! - call mld_ilu_bld(a,desc_a,p,upd,info,blck=blck) + call mld_ilu_bld(a,p,upd,info,blck=blck) if (info/=0) then call psb_errpush(4010,name,a_err='mld_ilu_bld') goto 9999 @@ -362,7 +392,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) if (info == 0) call psb_sp_clip(atmp,p%av(mld_ap_nd_),info,& & jmin=atmp%m+1,rscale=.false.,cscale=.false.) - + if (info == 0) call psb_spcnv(p%av(mld_ap_nd_),info,& & afmt='csr',dupl=psb_dupl_add_) if (info /= 0) then @@ -463,7 +493,7 @@ subroutine mld_dbjac_bld(a,desc_a,p,upd,info) ! ! Compute the LU factorization. ! - if (info == 0) call psb_ipcoo2csc(atmp,info,clshr=.true.) + if (info == 0) call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_) if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/mld_dilu_bld.f90 b/mlprec/mld_dilu_bld.f90 index a3380bc8..34eb2c8a 100644 --- a/mlprec/mld_dilu_bld.f90 +++ b/mlprec/mld_dilu_bld.f90 @@ -70,8 +70,6 @@ ! only the 'original' local part of the matrix to be factorized, ! i.e. the rows of the matrix held by the calling process according ! to the initial data distribution. -! desc_a - type(psb_desc_type), input. -! The communication descriptor associated to a. ! p - type(mld_dbaseprc_type), input/output. ! The 'base preconditioner' data structure. In input, p%iprcparm ! contains information on the type of factorization to be computed. @@ -88,7 +86,7 @@ ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see mld_bjac_bld), then blck does not contain any row. ! -subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck) +subroutine mld_dilu_bld(a,p,upd,info,blck) use psb_base_mod use mld_prec_mod, mld_protect_name => mld_dilu_bld @@ -98,7 +96,6 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck) ! Arguments type(psb_dspmat_type), intent(in), target :: a type(mld_dbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd integer, intent(out) :: info type(psb_dspmat_type), intent(in), optional :: blck @@ -116,7 +113,7 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck) call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(p%desc_data) call psb_info(ictxt, me, np) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' @@ -148,14 +145,14 @@ subroutine mld_dilu_bld(a,desc_a,p,upd,info,blck) end if endif - nrow_a = psb_cd_get_local_rows(desc_a) + nrow_a = psb_sp_get_nrows(a) nztota = psb_sp_get_nnzeros(a) if (present(blck)) then nztota = nztota + psb_sp_get_nnzeros(blck) end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ': out get_nnzeros',nztota,a%m,a%k + & ': out get_nnzeros',nztota,a%m,a%k,nrow_a n_row = p%desc_data%matrix_data(psb_n_row_) p%av(mld_l_pr_)%m = n_row diff --git a/mlprec/mld_dprec_aply.f90 b/mlprec/mld_dprec_aply.f90 index a46dd174..dc03b155 100644 --- a/mlprec/mld_dprec_aply.f90 +++ b/mlprec/mld_dprec_aply.f90 @@ -205,7 +205,6 @@ subroutine mld_dprec_aply1(prec,x,desc_data,info,trans) character(len=1), optional :: trans ! Local variables - character :: trans_ integer :: ictxt,np,me, err_act real(kind(1.d0)), pointer :: WW(:), w1(:) character(len=20) :: name @@ -217,11 +216,6 @@ subroutine mld_dprec_aply1(prec,x,desc_data,info,trans) ictxt = psb_cd_get_context(desc_data) call psb_info(ictxt, me, np) - if (present(trans)) then - trans_=trans - else - trans_='N' - end if allocate(ww(size(x)),w1(size(x)),stat=info) if (info /= 0) then @@ -231,7 +225,7 @@ subroutine mld_dprec_aply1(prec,x,desc_data,info,trans) goto 9999 end if - call mld_precaply(prec,x,ww,desc_data,info,trans_,work=w1) + call mld_precaply(prec,x,ww,desc_data,info,trans=trans,work=w1) if (info /= 0) then call psb_errpush(4010,name,a_err='mld_precaply') goto 9999 diff --git a/mlprec/mld_dslud_bld.f90 b/mlprec/mld_dslud_bld.f90 index 6e393fc5..0919039b 100644 --- a/mlprec/mld_dslud_bld.f90 +++ b/mlprec/mld_dslud_bld.f90 @@ -34,9 +34,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: mld_dslu_bld.f90 +! File: mld_dslud_bld.f90 ! -! Subroutine: mld_dslu_bld +! Subroutine: mld_dsludist_bld ! Version: real ! ! This routine computes the LU factorization of of a distributed matrix, diff --git a/mlprec/mld_dsp_renum.f90 b/mlprec/mld_dsp_renum.f90 index d299a1d3..95371bc7 100644 --- a/mlprec/mld_dsp_renum.f90 +++ b/mlprec/mld_dsp_renum.f90 @@ -61,8 +61,6 @@ ! part of the matrix to be reordered, i.e. the rows of the matrix ! held by the calling process according to the initial data ! distribution. -! desc_a - type(psb_desc_type), input. -! The communication descriptor associated to a. ! blck - type(psb_dspmat_type), input. ! The sparse matrix structure containing the remote rows of the ! matrix to be reordered, that have been retrieved by mld_asmat_bld @@ -81,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info) +subroutine mld_dsp_renum(a,blck,p,atmp,info) use psb_base_mod use mld_prec_mod, mld_protect_name => mld_dsp_renum @@ -92,7 +90,6 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info) type(psb_dspmat_type), intent(in) :: a,blck type(psb_dspmat_type), intent(inout) :: atmp type(mld_dbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info ! Local variables @@ -107,7 +104,7 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info) name='mld_dsp_renum' call psb_erractionsave(err_act) - ictxt=psb_cd_get_context(desc_a) + ictxt=psb_cd_get_context(p%desc_data) call psb_info(ictxt, me, np) ! @@ -127,12 +124,6 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info) if (p%iprcparm(mld_sub_ren_)==mld_renum_glb_) then - ! - ! Compute the row and column reordering coherent with the - ! global indices - ! - - mglob = psb_cd_get_global_rows(desc_a) ! ! Remember: we have switched IA1=COLS and IA2=ROWS. ! Now identify the set of distinct local column indices. @@ -144,9 +135,7 @@ subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info) goto 9999 end if - do k=1,nnr - itmp2(k) = p%desc_data%loc_to_glob(k) - enddo + call psb_loc_to_glob(itmp2(1:nnr),p%desc_data,info,iact='I') ! ! Compute reordering. We want new(i) = old(perm(i)). ! diff --git a/mlprec/mld_prec_mod.f90 b/mlprec/mld_prec_mod.f90 index d741dc3a..4e6cca06 100644 --- a/mlprec/mld_prec_mod.f90 +++ b/mlprec/mld_prec_mod.f90 @@ -45,29 +45,6 @@ module mld_prec_mod use mld_prec_type - interface mld_precbld - subroutine mld_dprecbld(a,desc_a,prec,info,upd) - use psb_base_mod - use mld_prec_type - implicit none - type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dprec_type), intent(inout) :: prec - integer, intent(out) :: info - character, intent(in),optional :: upd - end subroutine mld_dprecbld - subroutine mld_zprecbld(a,desc_a,prec,info,upd) - use psb_base_mod - use mld_prec_type - implicit none - type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_zprec_type), intent(inout) :: prec - integer, intent(out) :: info - character, intent(in),optional :: upd - end subroutine mld_zprecbld - end interface - interface mld_precinit subroutine mld_dprecinit(p,ptype,info,nlev) use psb_base_mod @@ -202,48 +179,35 @@ module mld_prec_mod end subroutine mld_zprec_aply1 end interface - interface mld_baseprc_bld - subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd) - use psb_base_mod - use mld_prec_type - type(psb_dspmat_type), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dbaseprc_type),intent(inout) :: p - integer, intent(out) :: info - character, intent(in), optional :: upd - end subroutine mld_dbaseprc_bld - subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd) - use psb_base_mod - use mld_prec_type - type(psb_zspmat_type), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_zbaseprc_type),intent(inout) :: p - integer, intent(out) :: info - character, intent(in), optional :: upd - end subroutine mld_zbaseprc_bld - end interface - - interface mld_mlprec_bld - subroutine mld_dmlprec_bld(a,desc_a,p,info) + interface mld_baseprec_aply + subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_base_mod use mld_prec_type - type(psb_dspmat_type), intent(inout), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_dbaseprc_type), intent(inout), target :: p - integer, intent(out) :: info - end subroutine mld_dmlprec_bld - subroutine mld_zmlprec_bld(a,desc_a,p,info) + type(psb_desc_type),intent(in) :: desc_data + type(mld_dbaseprc_type), intent(in) :: prec + real(kind(0.d0)),intent(in) :: x(:) + real(kind(0.d0)),intent(inout) :: 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 mld_dbaseprec_aply + subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_base_mod use mld_prec_type - type(psb_zspmat_type), intent(inout), target :: a - type(psb_desc_type), intent(in), target :: desc_a - type(mld_zbaseprc_type), intent(inout),target :: p - integer, intent(out) :: info - end subroutine mld_zmlprec_bld + type(psb_desc_type),intent(in) :: desc_data + type(mld_zbaseprc_type), intent(in) :: prec + complex(kind(1.d0)),intent(in) :: x(:) + complex(kind(1.d0)),intent(inout) :: 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 mld_zbaseprec_aply end interface - interface mld_baseprec_aply - subroutine mld_dbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + interface mld_as_aply + subroutine mld_das_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_base_mod use mld_prec_type type(psb_desc_type),intent(in) :: desc_data @@ -254,8 +218,8 @@ module mld_prec_mod character(len=1) :: trans real(kind(0.d0)),target :: work(:) integer, intent(out) :: info - end subroutine mld_dbaseprec_aply - subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + end subroutine mld_das_aply + subroutine mld_zas_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) use psb_base_mod use mld_prec_type type(psb_desc_type),intent(in) :: desc_data @@ -266,7 +230,7 @@ module mld_prec_mod character(len=1) :: trans complex(kind(1.d0)),target :: work(:) integer, intent(out) :: info - end subroutine mld_zbaseprec_aply + end subroutine mld_zas_aply end interface interface mld_mlprec_aply @@ -323,6 +287,91 @@ module mld_prec_mod end subroutine mld_zbjac_aply end interface + + interface mld_precbld + subroutine mld_dprecbld(a,desc_a,prec,info,upd) + use psb_base_mod + use mld_prec_type + implicit none + type(psb_dspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dprec_type), intent(inout) :: prec + integer, intent(out) :: info + character, intent(in),optional :: upd + end subroutine mld_dprecbld + subroutine mld_zprecbld(a,desc_a,prec,info,upd) + use psb_base_mod + use mld_prec_type + implicit none + type(psb_zspmat_type), intent(in), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_zprec_type), intent(inout) :: prec + integer, intent(out) :: info + character, intent(in),optional :: upd + end subroutine mld_zprecbld + end interface + + interface mld_baseprc_bld + subroutine mld_dbaseprc_bld(a,desc_a,p,info,upd) + use psb_base_mod + use mld_prec_type + type(psb_dspmat_type), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dbaseprc_type),intent(inout) :: p + integer, intent(out) :: info + character, intent(in), optional :: upd + end subroutine mld_dbaseprc_bld + subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd) + use psb_base_mod + use mld_prec_type + type(psb_zspmat_type), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_zbaseprc_type),intent(inout) :: p + integer, intent(out) :: info + character, intent(in), optional :: upd + end subroutine mld_zbaseprc_bld + end interface + + interface mld_as_bld + subroutine mld_das_bld(a,desc_a,p,upd,info) + use psb_base_mod + use mld_prec_type + type(psb_dspmat_type), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dbaseprc_type),intent(inout) :: p + character, intent(in) :: upd + integer, intent(out) :: info + end subroutine mld_das_bld + subroutine mld_zas_bld(a,desc_a,p,upd,info) + use psb_base_mod + use mld_prec_type + type(psb_zspmat_type), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_zbaseprc_type),intent(inout) :: p + character, intent(in) :: upd + integer, intent(out) :: info + end subroutine mld_zas_bld + end interface + + interface mld_mlprec_bld + subroutine mld_dmlprec_bld(a,desc_a,p,info) + use psb_base_mod + use mld_prec_type + type(psb_dspmat_type), intent(inout), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_dbaseprc_type), intent(inout), target :: p + integer, intent(out) :: info + end subroutine mld_dmlprec_bld + subroutine mld_zmlprec_bld(a,desc_a,p,info) + use psb_base_mod + use mld_prec_type + type(psb_zspmat_type), intent(inout), target :: a + type(psb_desc_type), intent(in), target :: desc_a + type(mld_zbaseprc_type), intent(inout),target :: p + integer, intent(out) :: info + end subroutine mld_zmlprec_bld + end interface + interface mld_diag_bld subroutine mld_ddiag_bld(a,desc_data,p,upd,info) use psb_base_mod @@ -345,43 +394,41 @@ module mld_prec_mod end interface interface mld_bjac_bld - subroutine mld_dbjac_bld(a,desc_data,p,upd,info) + subroutine mld_dbjac_bld(a,p,upd,info,data) use psb_base_mod use mld_prec_type - integer, intent(out) :: info type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data type(mld_dbaseprc_type), intent(inout) :: p character, intent(in) :: upd + integer, intent(out) :: info + integer, intent(in), optional :: data end subroutine mld_dbjac_bld - subroutine mld_zbjac_bld(a,desc_data,p,upd,info) + subroutine mld_zbjac_bld(a,p,upd,info,data) use psb_base_mod use mld_prec_type integer, intent(out) :: info type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data type(mld_zbaseprc_type), intent(inout) :: p character, intent(in) :: upd + integer, intent(in), optional :: data end subroutine mld_zbjac_bld end interface interface mld_ilu_bld - subroutine mld_dilu_bld(a,desc_data,p,upd,info,blck) + subroutine mld_dilu_bld(a,p,upd,info,blck) use psb_base_mod use mld_prec_type integer, intent(out) :: info type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data type(mld_dbaseprc_type), intent(inout) :: p character, intent(in) :: upd type(psb_dspmat_type), intent(in), optional :: blck end subroutine mld_dilu_bld - subroutine mld_zilu_bld(a,desc_data,p,upd,info,blck) + subroutine mld_zilu_bld(a,p,upd,info,blck) use psb_base_mod use mld_prec_type integer, intent(out) :: info type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data type(mld_zbaseprc_type), intent(inout) :: p character, intent(in) :: upd type(psb_zspmat_type), intent(in), optional :: blck @@ -538,22 +585,20 @@ module mld_prec_mod end interface interface mld_sp_renum - subroutine mld_dsp_renum(a,desc_a,blck,p,atmp,info) + subroutine mld_dsp_renum(a,blck,p,atmp,info) use psb_base_mod use mld_prec_type type(psb_dspmat_type), intent(in) :: a,blck type(psb_dspmat_type), intent(inout) :: atmp type(mld_dbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info end subroutine mld_dsp_renum - subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info) + subroutine mld_zsp_renum(a,blck,p,atmp,info) use psb_base_mod use mld_prec_type type(psb_zspmat_type), intent(in) :: a,blck type(psb_zspmat_type), intent(inout) :: atmp type(mld_zbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info end subroutine mld_zsp_renum end interface diff --git a/mlprec/mld_zas_aply.f90 b/mlprec/mld_zas_aply.f90 new file mode 100644 index 00000000..bc4c481c --- /dev/null +++ b/mlprec/mld_zas_aply.f90 @@ -0,0 +1,295 @@ +!!$ +!!$ +!!$ MLD2P4 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS v.2.0) +!!$ +!!$ (C) Copyright 2007 Alfredo Buttari University of Rome Tor Vergata +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ +!!$ 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 MLD2P4 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 MLD2P4 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. +!!$ +!!$ +! File: mld_zas_aply.f90 +! +! Subroutine: mld_zas_aply +! Version: real +! +! This routine applies the Additive Schwarz preconditioner by computing +! +! Y = beta*Y + alpha*op(K^(-1))*X, +! where +! - K is the base preconditioner, stored in prec, +! - op(K^(-1)) is K^(-1) or its transpose, according to the value of trans, +! - X and Y are vectors, +! - alpha and beta are scalars. +! +! +! Arguments: +! alpha - real(kind(0.d0)), input. +! The scalar alpha. +! prec - type(mld_dbaseprc_type), input. +! The base preconditioner data structure containing the local part +! of the preconditioner K. +! x - real(kind(0.d0)), dimension(:), input. +! The local part of the vector X. +! beta - real(kind(0.d0)), input. +! The scalar beta. +! y - real(kind(0.d0)), dimension(:), input/output. +! The local part of the vector Y. +! desc_data - type(psb_desc_type), input. +! The communication descriptor associated to the matrix to be +! preconditioned. +! trans - character, optional. +! If trans='N','n' then op(K^(-1)) = K^(-1); +! if trans='T','t' then op(K^(-1)) = K^(-T) (transpose of K^(-1)). +! work - real(kind(0.d0)), dimension (:), optional, target. +! Workspace. Its size must be at least 4*psb_cd_get_local_cols(desc_data). +! info - integer, output. +! Error code. +! +subroutine mld_zas_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_prec_mod, mld_protect_name => mld_zas_aply + + implicit none + + ! Arguments + type(psb_desc_type),intent(in) :: desc_data + type(mld_zbaseprc_type), intent(in) :: prec + complex(kind(0.d0)),intent(in) :: x(:) + complex(kind(0.d0)),intent(inout) :: y(:) + complex(kind(0.d0)),intent(in) :: alpha,beta + character(len=1) :: trans + complex(kind(0.d0)),target :: work(:) + integer, intent(out) :: info + + ! Local variables + integer :: n_row,n_col, int_err(5), nrow_d + complex(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) + integer :: ictxt,np,me,isz, err_act + character(len=20) :: name, ch_err + + name='mld_zas_aply' + info = 0 + call psb_erractionsave(err_act) + + ictxt = psb_cd_get_context(desc_data) + + call psb_info(ictxt, me, np) + + select case(trans) + case('N','n') + case('T','t','C','c') + case default + info=40 + int_err(1)=6 + ch_err(2:2)=trans + goto 9999 + end select + + select case(prec%iprcparm(mld_prec_type_)) + + case(mld_as_) + ! + ! Additive Schwarz preconditioner + ! + + if (prec%iprcparm(mld_n_ovr_)==0) then + ! + ! shortcut: this fixes performance for RAS(0) == BJA + ! + call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) + if(info /= 0) then + info=4010 + ch_err='psb_bjacaply' + goto 9999 + end if + + else + + ! + ! Note: currently trans is unused + ! + + n_row = psb_cd_get_local_rows(prec%desc_data) + n_col = psb_cd_get_local_cols(prec%desc_data) + nrow_d = psb_cd_get_local_rows(desc_data) + isz=max(n_row,N_COL) + if ((6*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + aux => work(3*isz+1:) + else if ((4*isz) <= size(work)) then + aux => work(1:) + allocate(ww(isz),tx(isz),ty(isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/3*isz,0,0,0,0/),& + & a_err='complex(kind(1.d0))') + goto 9999 + end if + else if ((3*isz) <= size(work)) then + ww => work(1:isz) + tx => work(isz+1:2*isz) + ty => work(2*isz+1:3*isz) + allocate(aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='complex(kind(1.d0))') + goto 9999 + end if + else + allocate(ww(isz),tx(isz),ty(isz),& + &aux(4*isz),stat=info) + if (info /= 0) then + call psb_errpush(4025,name,i_err=(/4*isz,0,0,0,0/),& + & a_err='complex(kind(1.d0))') + goto 9999 + end if + + endif + + tx(1:nrow_d) = x(1:nrow_d) + tx(nrow_d+1:isz) = dzero + + ! + ! Get the overlap entries of tx (tx==x) + ! + if (prec%iprcparm(mld_sub_restr_)==psb_halo_) then + call psb_halo(tx,prec%desc_data,info,work=aux,data=psb_comm_ext_) + if(info /=0) then + info=4010 + ch_err='psb_halo' + goto 9999 + end if + else if (prec%iprcparm(mld_sub_restr_) /= psb_none_) then + call psb_errpush(4001,name,a_err='Invalid mld_sub_restr_') + goto 9999 + end if + + ! + ! If required, reorder tx according to the row/column permutation of the + ! local extended matrix, stored into the permutation vector prec%perm + ! + if (prec%iprcparm(mld_sub_ren_)>0) then + call psb_gelp('n',prec%perm,tx,info) + if(info /=0) then + info=4010 + ch_err='psb_gelp' + goto 9999 + end if + endif + + ! + ! Apply to tx the block-Jacobi preconditioner/solver (multiple sweeps of the + ! block-Jacobi solver can be applied at the coarsest level of a multilevel + ! preconditioner). The resulting vector is ty. + ! + call mld_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,trans,aux,info) + if(info /= 0) then + info=4010 + ch_err='mld_bjac_aply' + goto 9999 + end if + + ! + ! Apply to ty the inverse permutation of prec%perm + ! + if (prec%iprcparm(mld_sub_ren_)>0) then + call psb_gelp('n',prec%invperm,ty,info) + if(info /= 0) then + info=4010 + ch_err='psb_gelp' + goto 9999 + end if + endif + + select case (prec%iprcparm(mld_sub_prol_)) + + case(psb_none_) + ! + ! Would work anyway, but since it is supposed to do nothing ... + ! call psb_ovrl(ty,prec%desc_data,info,& + ! & update=prec%iprcparm(mld_sub_prol_),work=aux) + + + case(psb_sum_,psb_avg_) + ! + ! Update the overlap of ty + ! + call psb_ovrl(ty,prec%desc_data,info,& + & update=prec%iprcparm(mld_sub_prol_),work=aux) + if(info /=0) then + info=4010 + ch_err='psb_ovrl' + goto 9999 + end if + + case default + call psb_errpush(4001,name,a_err='Invalid mld_sub_prol_') + goto 9999 + end select + + ! + ! Compute y = beta*y + alpha*ty (ty==K^(-1)*tx) + ! + call psb_geaxpby(alpha,ty,beta,y,desc_data,info) + + + if ((6*isz) <= size(work)) then + else if ((4*isz) <= size(work)) then + deallocate(ww,tx,ty) + else if ((3*isz) <= size(work)) then + deallocate(aux) + else + deallocate(ww,aux,tx,ty) + endif + end if + + case default + + call psb_errpush(4001,name,a_err='Invalid mld_prec_type_') + goto 9999 + + end select + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_errpush(info,name,i_err=int_err,a_err=ch_err) + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine mld_zas_aply + diff --git a/mlprec/mld_zasmat_bld.f90 b/mlprec/mld_zas_bld.f90 similarity index 56% rename from mlprec/mld_zasmat_bld.f90 rename to mlprec/mld_zas_bld.f90 index 4b02036b..6521646b 100644 --- a/mlprec/mld_zasmat_bld.f90 +++ b/mlprec/mld_zas_bld.f90 @@ -34,15 +34,13 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: mld_zasmat_bld.f90 +! File: mld_zas_bld.f90 ! -! Subroutine: mld_zasmat_bld +! Subroutine: mld_zas_bld ! Version: complex ! -! This routine builds the communication descriptor associated to the extended -! matrices that form the Additive Schwarz (AS) preconditioner and retrieves -! the remote pieces needed to build the local extended matrix. If the -! preconditioner is the block-Jacobi one, the routine makes only a copy of +! This routine builds the Additive Schwarz (AS) preconditioner. +! If the preconditioner is the block-Jacobi one, the routine makes only a copy of ! the descriptor of the original matrix. ! ! @@ -55,13 +53,9 @@ ! a - type(psb_zspmat_type), input. ! The sparse matrix structure containing the local part of the ! matrix to be preconditioned. -! blk - type(psb_zspmat_type), output. -! The sparse matrix structure containing the remote rows that -! extend the local matrix according to novr. If novr = 0 then -! blk does not contain any row. -! desc_data - type(psb_desc_type), input. +! desc_a - type(psb_desc_type), input. ! The communication descriptor of the sparse matrix a. -! desc_p - type(psb_desc_type), output. +! desc_p - type(psb_desc_type), output. ! The communication descriptor associated to the extended ! matrices that form the AS preconditioner. ! info - integer, output. @@ -71,32 +65,32 @@ ! preconditioner. Currently outfmt is set to 'CSR' by the ! calling routine mld_bjac_bld. ! -subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) +subroutine mld_zas_bld(a,desc_a,p,upd,info) use psb_base_mod - use mld_prec_mod, mld_protect_name => mld_zasmat_bld + use mld_prec_mod, mld_protect_name => mld_zas_bld Implicit None ! Arguments - integer, intent(in) :: ptype,novr - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_zspmat_type), Intent(out) :: blk - integer, intent(out) :: info - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - character(len=5), optional :: outfmt + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + type(mld_zbaseprc_type), intent(inout) :: p + character, intent(in) :: upd + integer, intent(out) :: info + + integer :: ptype,novr ! Local variables integer icomm Integer :: np,me,nnzero,& & ictxt, n_col,int_err(5),& - & tot_recv, n_row,nhalo, nrow_a,err_act + & tot_recv, n_row,nhalo, nrow_a,err_act, data_ integer :: debug_level, debug_unit character(len=20) :: name, ch_err - name='mld_zasmat_bld' + name='mld_as_bld' if(psb_get_errstatus() /= 0) return info=0 call psb_erractionsave(err_act) @@ -106,45 +100,21 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) If (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & ' start ', upd - ictxt = psb_cd_get_context(desc_data) - icomm = psb_cd_get_mpic(desc_data) + ictxt = psb_cd_get_context(desc_a) + icomm = psb_cd_get_mpic(desc_a) Call psb_info(ictxt, me, np) tot_recv=0 - nrow_a = desc_data%matrix_data(psb_n_row_) - nnzero = Size(a%aspk) - n_col = desc_data%matrix_data(psb_n_col_) + nrow_a = psb_cd_get_local_rows(desc_a) + n_col = psb_cd_get_local_cols(desc_a) + nnzero = psb_sp_get_nnzeros(a) nhalo = n_col-nrow_a + ptype = p%iprcparm(mld_prec_type_) + novr = p%iprcparm(mld_n_ovr_) select case (ptype) - case (mld_bjac_) - ! - ! Block-Jacobi preconditioner. Copy the descriptor, just in case - ! we want to renumber the rows and columns of the matrix. - ! - call psb_sp_all(0,0,blk,1,info) - if(info /= 0) then - info=4010 - ch_err='psb_sp_all' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - blk%fida = 'COO' - blk%infoa(psb_nnz_) = 0 - If (upd == 'F') Then - call psb_cdcpy(desc_data,desc_p,info) - If(debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' done cdcpy' - if(info /= 0) then - info=4010 - ch_err='psb_cdcpy' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif case(mld_as_) ! @@ -161,23 +131,12 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) ! ! Actually, this is just block Jacobi ! - If(debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ' calling allocate novr=0' - call psb_sp_all(0,0,blk,1,info) - if(info /= 0) then - info=4010 - ch_err='psb_sp_all' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - blk%fida = 'COO' - blk%infoa(psb_nnz_) = 0 + data_ = psb_no_comm_ if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & 'Calling desccpy' if (upd == 'F') then - call psb_cdcpy(desc_data,desc_p,info) + call psb_cdcpy(desc_a,p%desc_data,info) If(debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& & ' done cdcpy' @@ -191,63 +150,49 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) & write(debug_unit,*) me,' ',trim(name),& & 'Early return: P>=3 N_OVR=0' endif - return - endif + else + data_ = psb_comm_ext_ + If (upd == 'F') Then + ! + ! Build the auxiliary descriptor desc_p%matrix_data(psb_n_row_). + ! This is done by psb_cdbldext (interface to psb_cdovr), which is + ! independent of CSR, and has been placed in the tools directory + ! of PSBLAS, instead of the mlprec directory of MLD2P4, because it + ! might be used independently of the AS preconditioner, to build + ! a descriptor for an extended stencil in a PDE solver. + ! + call psb_cdbldext(a,desc_a,novr,p%desc_data,info,extype=psb_ovt_asov_) + if (info /= 0) then + info=4010 + ch_err='psb_cdbldext' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + Endif - If (upd == 'F') Then - ! - ! Build the auxiliary descriptor desc_p%matrix_data(psb_n_row_). - ! This is done by psb_cdbldext (interface to psb_cdovr), which is - ! independent of CSR, and has been placed in the tools directory - ! of PSBLAS, instead of the mlprec directory of MLD2P4, because it - ! might be used independently of the AS preconditioner, to build - ! a descriptor for an extended stencil in a PDE solver. - ! - call psb_cdbldext(a,desc_data,novr,desc_p,info,extype=psb_ovt_asov_) - if (info /= 0) then - info=4010 - ch_err='psb_cdbldext' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - Endif + End if if(debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ' From cdbldext _:',desc_p%matrix_data(psb_n_row_),& - & desc_p%matrix_data(psb_n_col_) - - ! - ! Retrieve the remote pieces needed to build the local extended matrix - ! - - n_row = desc_p%matrix_data(psb_n_row_) - - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'Before sphalo ',blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_) - Call psb_sphalo(a,desc_p,blk,info,& - & outfmt=outfmt,data=psb_comm_ext_,rowscale=.true.) - + & ' From cdbldext _:',p%desc_data%matrix_data(psb_n_row_),& + & p%desc_data%matrix_data(psb_n_col_) + + call mld_bjac_bld(a,p,upd,info,data=data_) if (info /= 0) then info=4010 - ch_err='psb_sphalo' + ch_err='mld_bjac_bld' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & 'After psb_sphalo ',& - & blk%fida,blk%m,psb_nnz_,blk%infoa(psb_nnz_) - case default + info=4001 ch_err='Invalid ptype' call psb_errpush(info,name,a_err=ch_err) goto 9999 - + End select if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),'Done' @@ -263,5 +208,5 @@ subroutine mld_zasmat_bld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) end if Return -End Subroutine mld_zasmat_bld +End Subroutine mld_zas_bld diff --git a/mlprec/mld_zbaseprec_aply.f90 b/mlprec/mld_zbaseprec_aply.f90 index eab5cfc8..9d234a27 100644 --- a/mlprec/mld_zbaseprec_aply.f90 +++ b/mlprec/mld_zbaseprec_aply.f90 @@ -96,9 +96,8 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! Local variables integer :: n_row,n_col, int_err(5), nrow_d complex(kind(1.d0)), pointer :: ww(:), aux(:), tx(:),ty(:) - character ::diagl, diagu - integer :: ictxt,np,me,isz, err_act - character(len=20) :: name, ch_err + integer :: ictxt,np,me,isz, err_act + character(len=20) :: name, ch_err name='mld_zbaseprec_aply' info = 0 @@ -108,9 +107,6 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) call psb_info(ictxt, me, np) - diagl='U' - diagu='U' - select case(trans) case('N','n') case('T','t','C','c') @@ -124,17 +120,17 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) select case(prec%iprcparm(mld_prec_type_)) case(mld_noprec_) - ! - ! No preconditioner - ! + ! + ! No preconditioner + ! call psb_geaxpby(alpha,x,beta,y,desc_data,info) case(mld_diag_) - ! - ! Diagonal preconditioner - ! - + ! + ! Diagonal preconditioner + ! + if (size(work) >= size(x)) then ww => work else @@ -158,21 +154,21 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) end if case(mld_bjac_) - ! - ! Block-Jacobi preconditioner - ! + ! + ! Block-Jacobi preconditioner + ! call mld_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) if(info /= 0) then - info=4010 - ch_err='mld_bjac_aply' - goto 9999 + info=4010 + ch_err='mld_bjac_aply' + goto 9999 end if case(mld_as_) - ! - ! Additive Schwarz preconditioner - ! + ! + ! Additive Schwarz preconditioner + ! if (prec%iprcparm(mld_n_ovr_)==0) then ! @@ -289,8 +285,9 @@ subroutine mld_zbaseprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) case(psb_none_) ! ! Would work anyway, but since it is supposed to do nothing ... - ! call f90_psovrl(ty,prec%desc_data,update=prec%a_restrict) - ! + ! call psb_ovrl(ty,prec%desc_data,info,& + ! & update=prec%iprcparm(mld_sub_prol_),work=aux) + case(psb_sum_,psb_avg_) ! diff --git a/mlprec/mld_zbaseprec_bld.f90 b/mlprec/mld_zbaseprec_bld.f90 index ba4fd521..0350bede 100644 --- a/mlprec/mld_zbaseprec_bld.f90 +++ b/mlprec/mld_zbaseprec_bld.f90 @@ -149,7 +149,31 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd) goto 9999 end if - case(mld_bjac_,mld_as_) + case(mld_bjac_) + + call mld_check_def(p%iprcparm(mld_sub_ren_),'renumbering',& + & mld_renum_none_,is_legal_renum) + call mld_check_def(p%iprcparm(mld_sub_solve_),'fact',& + & mld_ilu_n_,is_legal_ml_fact) + + 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 + + ! Build the local part of the base preconditioner + call mld_bjac_bld(a,p,iupd,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='mld_bjac_bld') + goto 9999 + end if + + + case(mld_as_) ! Block Jacobi and additive Schwarz preconditioners/smoothers call mld_check_def(p%iprcparm(mld_n_ovr_),'overlap',& @@ -165,7 +189,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd) ! Set parameters for using SuperLU_dist on the local submatrices if (p%iprcparm(mld_sub_solve_)==mld_sludist_) then - p%iprcparm(mld_n_ovr_) = 0 + p%iprcparm(mld_n_ovr_) = 0 p%iprcparm(mld_smooth_sweeps_) = 1 end if @@ -174,7 +198,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd) & ': Calling mld_bjac_bld' ! Build the local part of the base preconditioner - call mld_bjac_bld(a,desc_a,p,iupd,info) + call mld_as_bld(a,desc_a,p,iupd,info) if(info /= 0) then info=4010 call psb_errpush(info,name,a_err='mld_bjac_bld') @@ -182,6 +206,7 @@ subroutine mld_zbaseprc_bld(a,desc_a,p,info,upd) end if case default + info=4001 ch_err='Unknown mld_prec_type_' call psb_errpush(info,name,a_err=ch_err) diff --git a/mlprec/mld_zbjac_bld.f90 b/mlprec/mld_zbjac_bld.f90 index c00317e9..41d04047 100644 --- a/mlprec/mld_zbjac_bld.f90 +++ b/mlprec/mld_zbjac_bld.f90 @@ -91,8 +91,14 @@ ! part of the preconditioner or solver at the current level. ! info - integer, output. ! Error code. +! data - integer Which index list in desc_a should be used to retrieve +! rows, default psb_comm_halo_ +! psb_no_comm_ don't retrieve remote rows +! psb_comm_halo_ use halo_index +! psb_comm_ext_ use ext_index +! psb_comm_ovrl_ DISABLED for this routine. ! -subroutine mld_zbjac_bld(a,desc_a,p,upd,info) +subroutine mld_zbjac_bld(a,p,upd,info,data) use psb_base_mod use mld_prec_mod, mld_protect_name => mld_zbjac_bld @@ -101,10 +107,10 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) ! Arguments type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in) :: desc_a type(mld_zbaseprc_type), intent(inout) :: p integer, intent(out) :: info character, intent(in) :: upd + integer, intent(in), optional :: data ! Local Variables integer :: i, k, m @@ -112,9 +118,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) character :: trans, unitd type(psb_zspmat_type) :: blck, atmp integer :: debug_level, debug_unit - integer :: err_act, n_row, nrow_a,n_col + integer :: err_act, n_row, nrow_a,n_col, data_ integer :: ictxt,np,me - character(len=20) :: name + character(len=20) :: name, ch_err character(len=5), parameter :: coofmt='COO', csrfmt='CSR' if(psb_get_errstatus().ne.0) return @@ -123,7 +129,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(p%desc_data) call psb_info(ictxt, me, np) m = a%m @@ -151,25 +157,50 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),': Start',& & p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_) - - ! - ! Build the communication descriptor for the Additive Schwarz - ! preconditioner and retrieves the remote pieces of the local - ! extended matrix needed by that preconditioner. If the - ! preconditioner is the block-Jacobi one, only a copy of the - ! descriptor of the original matrix is made. - ! - call mld_asmat_bld(p%iprcparm(mld_prec_type_),p%iprcparm(mld_n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info,outfmt=csrfmt) - - if (info/=0) then - call psb_errpush(4010,name,a_err='mld_asmat_bld') - goto 9999 + if (present(data)) then + data_ = data + else + data_ = psb_no_comm_ end if - if (debug_level >= psb_debug_outer_) & - & write(debug_unit,*) me,' ',trim(name),& - & ': out of mld_asmat_bld' + select case (data_) + case(psb_no_comm_) + If (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & ' calling allocate novr=0' + call psb_sp_all(0,0,blck,1,info) + if(info /= 0) then + info=4010 + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + blck%fida = 'COO' + blck%infoa(psb_nnz_) = 0 + + case(psb_comm_halo_,psb_comm_ext_) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'Before sphalo ',blck%fida,blck%m,psb_nnz_,blck%infoa(psb_nnz_) + Call psb_sphalo(a,p%desc_data,blck,info,data=data_,rowscale=.true.) + + if (info /= 0) then + info=4010 + ch_err='psb_sphalo' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'After psb_sphalo ',& + & blck%fida,blck%m,psb_nnz_,blck%infoa(psb_nnz_) + case default + info=35 + call psb_errpush(info,name,i_err=(/5,data,0,0,0/)) + goto 9999 + end select ! ! Treat separately the case the local matrix has to be reordered @@ -177,9 +208,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) ! select case(p%iprcparm(mld_sub_ren_)) - ! - ! A reordering of the local matrix is required. - ! + ! + ! A reordering of the local matrix is required. + ! case (1:) ! @@ -187,7 +218,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) ! according to the value of p%iprcparm(sub_ren_). The reordered ! matrix is stored into atmp, using the COO format. ! - call mld_sp_renum(a,desc_a,blck,p,atmp,info) + call mld_sp_renum(a,blck,p,atmp,info) if (info/=0) then call psb_errpush(4010,name,a_err='mld_sp_renum') goto 9999 @@ -233,7 +264,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) ! ILU(k)/MILU(k)/ILU(k,t) factorization. ! call psb_spcnv(atmp,info,afmt='csr',dupl=psb_dupl_add_) - if (info == 0) call mld_ilu_bld(atmp,p%desc_data,p,upd,info) + if (info == 0) call mld_ilu_bld(atmp,p,upd,info) if (info/=0) then call psb_errpush(4010,name,a_err='mld_ilu_bld') goto 9999 @@ -282,9 +313,9 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) goto 9999 end if - ! - ! No reordering of the local matrix is required - ! + ! + ! No reordering of the local matrix is required + ! case(0) ! @@ -321,7 +352,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) call psb_errpush(4010,name,a_err='clip & psb_spcnv csr 4') goto 9999 end if - + k = psb_sp_get_nnzeros(p%av(mld_ap_nd_)) call psb_sum(ictxt,k) @@ -338,7 +369,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) ! ! Compute the incomplete LU factorization. ! - call mld_ilu_bld(a,desc_a,p,upd,info,blck=blck) + call mld_ilu_bld(a,p,upd,info,blck=blck) if (info/=0) then call psb_errpush(4010,name,a_err='mld_ilu_bld') goto 9999 @@ -463,7 +494,7 @@ subroutine mld_zbjac_bld(a,desc_a,p,upd,info) ! ! Compute the LU factorization. ! - if (info == 0) call psb_ipcoo2csc(atmp,info,clshr=.true.) + if (info == 0) call psb_spcnv(atmp,info,afmt='csc',dupl=psb_dupl_add_) if (info == 0) call mld_umf_bld(atmp,p%desc_data,p,info) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& diff --git a/mlprec/mld_zilu_bld.f90 b/mlprec/mld_zilu_bld.f90 index 0a5290a6..679a7cc7 100644 --- a/mlprec/mld_zilu_bld.f90 +++ b/mlprec/mld_zilu_bld.f90 @@ -70,8 +70,6 @@ ! only the 'original' local part of the matrix to be factorized, ! i.e. the rows of the matrix held by the calling process according ! to the initial data distribution. -! desc_a - type(psb_desc_type), input. -! The communication descriptor associated to a. ! p - type(mld_zbaseprc_type), input/output. ! The 'base preconditioner' data structure. In input, p%iprcparm ! contains information on the type of factorization to be computed. @@ -88,7 +86,7 @@ ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see mld_bjac_bld), then blck does not contain any row. ! -subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck) +subroutine mld_zilu_bld(a,p,upd,info,blck) use psb_base_mod use mld_prec_mod, mld_protect_name => mld_zilu_bld @@ -98,7 +96,6 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck) ! Arguments type(psb_zspmat_type), intent(in), target :: a type(mld_zbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd integer, intent(out) :: info type(psb_zspmat_type), intent(in), optional :: blck @@ -116,7 +113,7 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck) call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - ictxt = psb_cd_get_context(desc_a) + ictxt = psb_cd_get_context(p%desc_data) call psb_info(ictxt, me, np) if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' @@ -148,14 +145,14 @@ subroutine mld_zilu_bld(a,desc_a,p,upd,info,blck) end if endif - nrow_a = psb_cd_get_local_rows(desc_a) + nrow_a = psb_sp_get_nrows(a) nztota = psb_sp_get_nnzeros(a) if (present(blck)) then nztota = nztota + psb_sp_get_nnzeros(blck) end if if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),& - & ': out get_nnzeros',nztota,a%m,a%k + & ': out get_nnzeros',nztota,a%m,a%k,nrow_a n_row = p%desc_data%matrix_data(psb_n_row_) p%av(mld_l_pr_)%m = n_row diff --git a/mlprec/mld_zprec_aply.f90 b/mlprec/mld_zprec_aply.f90 index cca301a9..aed207d6 100644 --- a/mlprec/mld_zprec_aply.f90 +++ b/mlprec/mld_zprec_aply.f90 @@ -206,7 +206,6 @@ subroutine mld_zprec_aply1(prec,x,desc_data,info,trans) character(len=1), optional :: trans ! Local variables - character :: trans_ integer :: ictxt,np,me, err_act complex(kind(1.d0)), pointer :: WW(:), w1(:) character(len=20) :: name @@ -218,11 +217,6 @@ subroutine mld_zprec_aply1(prec,x,desc_data,info,trans) ictxt = psb_cd_get_context(desc_data) call psb_info(ictxt, me, np) - if (present(trans)) then - trans_=trans - else - trans_='N' - end if allocate(ww(size(x)),w1(size(x)),stat=info) if (info /= 0) then @@ -232,7 +226,7 @@ subroutine mld_zprec_aply1(prec,x,desc_data,info,trans) goto 9999 end if - call mld_precaply(prec,x,ww,desc_data,info,trans_,work=w1) + call mld_precaply(prec,x,ww,desc_data,info,trans=trans,work=w1) if (info /= 0) then call psb_errpush(4010,name,a_err='mld_precaply') goto 9999 diff --git a/mlprec/mld_zslud_bld.f90 b/mlprec/mld_zslud_bld.f90 index 435b0177..afea847e 100644 --- a/mlprec/mld_zslud_bld.f90 +++ b/mlprec/mld_zslud_bld.f90 @@ -34,9 +34,9 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: mld_zslu_bld.f90 +! File: mld_zslud_bld.f90 ! -! Subroutine: mld_zslu_bld +! Subroutine: mld_zsludist_bld ! Version: real ! ! This routine computes the LU factorization of of a distributed matrix, diff --git a/mlprec/mld_zsp_renum.f90 b/mlprec/mld_zsp_renum.f90 index 5d7a79a9..7dff9275 100644 --- a/mlprec/mld_zsp_renum.f90 +++ b/mlprec/mld_zsp_renum.f90 @@ -61,8 +61,6 @@ ! part of the matrix to be reordered, i.e. the rows of the matrix ! held by the calling process according to the initial data ! distribution. -! desc_a - type(psb_desc_type), input. -! The communication descriptor associated to a. ! blck - type(psb_zspmat_type), input. ! The sparse matrix structure containing the remote rows of the ! matrix to be reordered, that have been retrieved by mld_asmat_bld @@ -81,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info) +subroutine mld_zsp_renum(a,blck,p,atmp,info) use psb_base_mod use mld_prec_mod, mld_protect_name => mld_zsp_renum @@ -92,7 +90,6 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info) type(psb_zspmat_type), intent(in) :: a,blck type(psb_zspmat_type), intent(inout) :: atmp type(mld_zbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a integer, intent(out) :: info ! Local variables @@ -107,7 +104,7 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info) name='mld_zsp_renum' call psb_erractionsave(err_act) - ictxt=psb_cd_get_context(desc_a) + ictxt=psb_cd_get_context(p%desc_data) call psb_info(ictxt, me, np) ! @@ -127,12 +124,6 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info) if (p%iprcparm(mld_sub_ren_)==mld_renum_glb_) then - ! - ! Compute the row and column reordering coherent with the - ! global indices - ! - - mglob = psb_cd_get_global_rows(desc_a) ! ! Remember: we have switched IA1=COLS and IA2=ROWS. ! Now identify the set of distinct local column indices. @@ -144,9 +135,7 @@ subroutine mld_zsp_renum(a,desc_a,blck,p,atmp,info) goto 9999 end if - do k=1,nnr - itmp2(k) = p%desc_data%loc_to_glob(k) - enddo + call psb_loc_to_glob(itmp2(1:nnr),p%desc_data,info,iact='I') ! ! Compute reordering. We want new(i) = old(perm(i)). ! diff --git a/test/fileread/runs/drt.sh b/test/fileread/runs/drt.sh index 59f49ba6..0a33f865 100755 --- a/test/fileread/runs/drt.sh +++ b/test/fileread/runs/drt.sh @@ -16,7 +16,7 @@ do # 3rd batch: 14bis,64bis 4 sweeps echo "mpirun -np $np -machinefile locm df_bench >>log.part$part.ren$renum.${np}p" #/usr/local/mpich-gcc42/bin/mpirun -np $np -machinefile locm df_bench <>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <